DFT-Based Identification of Oscillation Modes from PMU Data Using an Exponential Window Function

The characteristics of oscillation modes, such as interarea, regional, and subsynchronous modes, can vary during a power system fault, which can cause switching and control actions in the power system. Transient data of the modal response due to such a fault can be acquired through phasor measurement units (PMUs). When the transient data have a long duration, it is desirable to perform modal identification separately on each segment of the transient data, so as to reflect the varying characteristics of oscillation modes. A conventional discrete Fourier transform (DFT)-based method for parametric modal identification cannot be efficiently applied to such a segment dataset. In this paper, a DFT-based method with an exponential window function is proposed to identify oscillation modes from each segment of transient data in PMUs. This window function allows the application of the least squares method (LSM) for modal identification in the same manner as the conventional method. The accuracy of identification of the proposed method is compared with those of the conventional method and a Prony method through synthetic data of transient responses. Its feasibility is also verified by identifying real-world oscillation modes from transient data in PMUs.


Introduction
Power systems can have various oscillation modes, such as interarea, regional, and subsynchronous modes. An interarea mode has an oscillation frequency of less than 1 Hz and oscillates in opposite directions between faraway regions in the power system [1]. A regional mode has an oscillation frequency between 1 and 2 Hz and oscillates among generators in a region of the power system. A subsynchronous mode has an oscillation frequency above 5 Hz and oscillates because of interactions among the mechanical turbine generator system, the electrical power transmission system of series compensation, and electronic devices in high-voltage DC (HVDC) transmission systems and flexible AC transmission systems (FACTS) [2,3]. These oscillations can have a severe effect on the stability and transfer capacity of the power system. It is essential to monitor and control the oscillation modes for efficient operation of the power system.
Waveforms of modal oscillations can be observed well in a transient response of power systems owing to a sudden disturbance caused by a fault in the power system or a large change in the power load [4][5][6]. The Fourier transform magnitude of this transient response has a distinguishable peak at the oscillation frequency of each mode, owing to low modal damping. Nowadays, the modal oscillation frequencies and damping factors can be estimated from the transient data, which can be obtained from measurements of the active power, phasor angle, and frequency in phasor measurement units (PMUs). Most PMUs have data measurement rates below 60 per second. Such PMU data are not suitable for identifying a subsynchronous oscillation mode with an oscillation frequency of up to 60 Hz. Some PMUs have data measurement rates above 60 per second: their data enable subsynchronous modes with low-frequency oscillations to be identified effectively [7]. Currently, PMUs are installed on transmission grids. Their installation is now expanding to distributed grids and microgrids for monitoring and protection [8,9]. PMU data on distribution grids and microgrids generally have signal-to-noise ratios (SNRs) lower than transmission grids.
Parametric identification of oscillation modes can be performed on transient data in the time domain or frequency domain. A variety of time-domain-based methods, such as Prony [10], Matrix pencil [11], and Hankel total least squares (HTLS) methods, have been reported so far [12]. They fit transient data to their autoregressive (AR) models of oscillation modes directly. These methods are very sensitive to noise, which inevitably means that many fake modes of noise are added to the AR model. The fake mode number can be assigned through singular value decomposition (SVD) [10], which requires a little computation time. A discrete Fourier transform (DFT)-based method of parametric modal identification was designed to curve-fit DFT coefficients of transient data into a transfer function of oscillation modes in the frequency domain [13,14]. Such curve-fitting is performed on small frequency ranges around each modal peak in the DFT magnitude, which can lead to a short computation time compared to that involved in the time-domain methods. However, this method can be applied for identification of stable modes whose oscillations converge to zeros in the transient data.
Oscillation phenomena in power systems can sometimes last for a long duration. In such transient data, oscillation modes can be unstable or typically exhibit non-stationary behaviors due to the nature of switching events and other control actions in the power system [15,16]. This means that the oscillation frequency and damping factors may vary slightly from time to time in transient data. For this reason, it is necessary to split such transient data into many segments and identify oscillation modes from each segment sequentially. A smaller segment is more desirable as transient data are less stationary. A conventional DFT-based method cannot be applied directly to the segment where modal oscillations do not converge to zero. For the application of the conventional DFT-based method, a window function is required to force the modal oscillations to converge to zero in the segment.
In this paper, a DFT-based method with an exponential window function is proposed to identify oscillation modes from each segment of transient data in PMUs. This exponential function enables the conventional DFT-based method to be applied to each segment of transient data. The identification accuracy of the proposed methods is compared with conventional DFT-based and Prony methods through modal identification using synthetic data of transient responses. Its feasibility is demonstrated by identifying oscillation modes from PMU transient data of power systems in the United States and Central America.

Models for Modal Identification in Time and Frequency Domains
An oscillation mode in the power system is triggered by a sudden disturbance caused by a power event such as a line fault, generator trip, or huge load change. This disturbance can be modeled as an impulse function at the beginning of the power event. A modal oscillation owing to the disturbance can be expressed as a second order differential equation. When there are K oscillation modes in the power system, their transient response after the event can be given as the following sum y(t) of sinusoidal functions with the magnitudes of exponential functions: where q k (t) is the k-th mode with a damping factor of ξ k and an oscillation frequency of ω k = 2π f k , f 1 < f 2 · · · < f K . Here, A k and θ k are the amplitude and phase angle of the k-th mode, respectively.
The natural frequency and damping ratio of the k-th mode can be calculated to be k = ω k 1 − ζ 2 k and ζ k = ξ k / k , respectively. The initial time of a power event is denoted as t = 0 and t e is the duration of the transient response. The modal parameters ξ k , f k , A k , and θ k can vary slightly from time to time, according to control actions of the power system after the event.
The Laplace transform of y(t) is given as the following for Y(s): where E k (s) has a transcendental function of e −(s+ξ k )t e . Here, Y(s) becomes a rational function if E k (s) is ignorable, which is satisfied under the conditions of stable oscillation modes with e −ξ k t e ≈ 0, t e > t. The Fourier transform of y(t) is obtained as Y( jω), where j = √ −1. The modal parameters ξ k , f k , A k , and θ k can be estimated in the frequency domain by curve-fitting Y( jω) into DFT coefficients of transient data in PMUs.
For modal identification in the time domain, the continuous-time signal y(t) should be converted into a discrete-time signal, which is represented with an autoregressive (AR) model. The impulse invariance transform [17] is applied to this conversion with a sampling interval of T s . Then, the AR model is given as where y[n] = y(nT s ) and N is the data number for t e = NT s . The Prony method can estimate the modal parameters ξ k , f k , A k , and θ k by fitting the AR model into transient data in PMUs. Actually, PMU data contain measurement noise [18] and a DC component [19,20]. The measurement noise can generally be assumed to be white noise. The DC component is an instantaneous drop or rise that reflects any imbalances between the supply and demand of active power in the power system. The DC component as well as the measurement noise degrades the identification accuracy of the time-domain-based method more severely than that of the frequency-domain-based method. In this paper, a second-order Butterworth digital bandpass filter is introduced for pre-processing to reject the DC component from transient data in PMUs. Its lower cutoff frequency is designed to be 0.01 Hz, which is far away from the lowest frequency of interarea oscillations. Its upper cutoff frequency is set at ( f s /2 − 0.01) Hz, where f s = 1/T s . This bandwidth leads to ignoring the filtering effects on y[n] and the measurement noise. Then, a pre-filtered signal from transient data in PMUs can be given as the following for y w [n]: where w[n] is the measurement noise. The magnitude of Y( jω) has sharp peaks around ω k due to low modal damping. This indicates that the modal number K can be determined by counting the number of dominant peaks in the DFT magnitude of y w [n]. The shape of these peaks can be more apparent as the DFT frequency resolution is smaller; the resolution is defined as the frequency interval between two adjacent DFT coefficients. For this reason, zero data values are padded properly to y w [n]. Then, the DFT coefficient of y w [n] was calculated as the following for Y w [m]: where M is the total data number, including the zero-padding data, and the index m corresponds to m/(MT s ) Hz in the frequency domain. Frequencies at the peaks are measured as f k , k = 1, 2, · · · , K, and are rearranged as ( f 1 < f 2 < · · · < f K ) according to their size. They can be used effectively for the discrimination of fake modes in the time-domain-based method and for assignment of curve-fitting frequency ranges per oscillation mode in the frequency-domain-based method.

Conventional DFT-Based Method
A conventional DFT-based method can be applied to Y( jω) only if Y(s) is a rational function, which is satisfied under the condition of E k (s) ≈ 0 in Equation (3). In this condition, Y(s) can be represented with 4K coefficients of a k and b k , k = 1, 2, · · · , 2K as Each Q k (s) is obtained from a k and b k by factorizing Y(s) through its characteristic root α k of the k-th oscillation mode as where β k is a complex number corresponding to α k and * represents the conjugate of a complex number. The k-th modal parameters of ω k = 2π f k , ξ k , A k , and θ k in Equation (1) can be calculated from α k and β k as follows: where |·| and ∠(·) represent the magnitude and angle of a complex number, respectively. The K oscillation modes in Y(s) of Equation (7) are identified through estimation of the coefficients a k and b k , which is accomplished by curve-fitting Y w [m] of Equation (6) into Y( jω) of Equation (7) in the frequency domain. To apply the LSM for estimation of a k and b k , the complex values of Equation (7) by s = jω can be divided into the two following linear equations for a k and b k in the real and imaginary parts: where Re{·} and Im{·} represent the real and imaginary parts of a complex number, respectively.
To reduce the noise effect as well as the computation burden for modal identification, it is desirable to perform curve-fitting on Y w [m] around f k , where most of the modal oscillation power is concentrated. Thus, Y w [m] of the following frequency ranges can be applied for the curve-fitting: where f c is a design parameter. Oscillation frequencies of interarea and regional modes can be distributed at an interval of up to 0.1 Hz in the frequency domain. Subsynchronus frequencies can be typically distributed farther than 1 Hz from each other. For these reasons, the curve-fitting frequency range for interarea modes and regional modes is set to f c = 0.05 Hz, and that for synchronous modes is set to f c = 0.5 Hz in this paper. Because a mode has four unknown parameters, namely ω k , ξ k , A k , and θ k , the sample number M c of the angular frequency ω in Equations (10) and (11) per mode should be greater than 2. The number M c was chosen experimentally to be 5 by comparing the identification accuracy obtained during the simulations discussed in Section 3.2. Then, M is calculated from M c = 2M f c T s + 1. When angular frequency samples per mode in Equation (12) are given as ω k,l , 0 ≤ l < M C , the curve-fitting of Y( jω) of Equations (10) and (11) represents the following matrix form [13]: The LSM is applied to Equation (13) by replacing Y jω k,l in B with Y w [m] in Equation (6), where the index m corresponds to ω k,l in the frequency domain. Estimates of a k and b k are obtained as Here, α k and β k in Equation (8) are obtained from the estimated a k and b k . The k-th modal parameters are estimated from α k and β k asf k ,ξ k ,Â k , andθ k through Equation (9).

Prony Method with SVD
In a time-domain-based method, y w [n] in Equation (4) is expressed with the coefficient number of 2K for an AR model of oscillation modes. A Prony method estimates the AR model coefficients using the LSM, whose estimation error is extremely dependent on the noise w[n] in y w [n]. The estimation error can be reduced significantly by adding fake modes of the noise to the AR model. When the number of characteristic roots in y w [n] is assigned to P 2K, including fake modes, y w [n] can be represented with unknown coefficients c p , p = 1, 2, · · · , P, employing the following AR model: In this paper, P is set to the largest integer smaller than N/3 [21]. The coefficients c p are obtained from the following matrix equation [10,13]: . . .
To reduce the noise effect in modal identification, SVD can be performed on the Hankel matrix H, whose singular values (SVs) can be partitioned into the following form: where U 1 U 2 and V 1 V 2 are the left and the right singular vectors, respectively. Here, S 1 and S 2 are the diagonal matrices with SVs, σ 1 ≥ σ 2 · · · ≥ σ P ; ρ is a design parameter that functions as the rank threshold for choice of P, based on which the matrix dimensions of U 1 , V 1 , and S 1 are given as 2P × P, P × P, and P × P, respectively [22]. In this paper, ρ is assigned to 0.99 [10]. By truncating S 2 of the smallest SVs, the coefficients c p are estimated by applying the LSM to U 1 S 1 V T 1 C = G. Characteristic roots of the AR model of Equation (17) are obtained from estimates of c p . Their oscillation frequencies and damping factors can be calculated as ω p = 2π f p and ξ p , respectively: where λ p is the characteristic root.
The k-th modal characteristic root γ k can be chosen to be λ p , which minimizes f k − f p , p = 1, 2, · · · , P. The k-th modal oscillation frequency and damping factor are estimated with the chosen p aŝ f k = f p andξ k = ξ p , respectively. Here, y w [n] can be represented with complex magnitudes µ k and µ * k of the k-th mode, corresponding to γ k and γ * k as The magnitude µ k is calculated by applying the LSM to the above equation. Estimates of A k and θ k in Equation (1) can be obtained aŝ

Exponential Window Function for Modal Identification
A transient response of oscillation modes with a long duration can be nonstationary, owing to time-dependent control and nonlinear dynamics in the power system. To make the modal identification reflect a nonstationary phenomenon, the identification needs to be performed separately for sequential segments of the transient data. Unfortunately, the conventional DFT-based method has difficulty identifying the modes in segments where modal oscillations do not reach zero exponentially. This difficulty originates from Laplace transform of such a segment, which contains the transcendental function E k (s) 0 in Equation (3).
To overcome the limited use of the conventional DFT-based method, the proposed method introduces a window function into the Laplace transform. When a portion of y(t) is given as x(τ), Energies 2019, 12, 4357 7 of 16 0 ≤ τ < T h for modal identification, the following exponential window function h(τ) is introduced to the Laplace transform of x(τ): where ξ 0 > 0 is a design parameter. The Laplace transform of x(τ) = h(τ)x(τ) is given as the following X(s): When h(T h ) = e −ξ o T h ≈ 0, E k (s) is negligible, which makes X(s) approximate to a rational function with coefficients a k and b k , as follows: The closer that h(T h ) gets to zero, the more X(s) approaches a rational function. However, the modal oscillation intensity in x(τ) decreases, which can degrade the identification accuracy of modes. For this reason, in this paper, the approximation of X(s) to the rational function is achieved by assigning h(T h ) to 1 % of h(0) = 1, which leads to ξ 0 T h = 4.61. The assignment of h(T h ) needs to be studied further.
The function h(τ) acts as a kind of weight function for x(τ). The effect of this weight was investigated by dividing h(τ) into the first part of 0 ≤ τ < T h /2 and the second part of T h /2 ≤ τ < T h . The power ratio of the second part to the whole h(τ) is calculated as This power ratio implies that the second part rarely contributes to estimation of a k and b k in Equation (27). The effectual identification interval due to h(τ) can be regarded as the first part. The second part is used accessorily in modal identification. For this reason, in the proposed method, modal identification is done for each segment of the interval Γ = T h /2 and uses that segment and the next segment. These two segments are given as x(τ). On the other hand, the conventional DFT-based and Prony methods use one segment.
A shorter segment can reflect the nonstationary phenomenon of oscillation modes more effectively in modal identification. The oscillation mode can be identified properly from a segment whose interval is at least longer than two cycles of its oscillation frequency in the conventional DFT-based and Prony methods. In this paper, the interval Γ for identification of interarea and regional modes is set to be about two cycles of the lowest oscillation frequency f 1 , and for identification of subsynchronus modes is set to be one second.
The transient data y w [n] in PMUs are segmented at every interval of Γ. In the proposed method, a segment and its next segment are combined into x[l], 0 ≤ l < L = 2Γ/T s , corresponding to , was calculated by padding (M − L)-zero data. The coefficients a k and b k in Equation (27) can be estimated by the same curve-fitting as in Equation (10) through Equation (12), where X( jω) is substituted with Y( jω). A modal oscillation signal in Equation (1) and the exponential window function can be unified as follows: The rational function of X(s) can be factorized through its characteristic roots α k , corresponding to q k (τ) as where β k is a complex number corresponding to α k . Then, α k and β k are obtained from a k and b k . The modal parameters of ω k = 2π f k , ξ k , A k , and θ k are estimated from α k and β k as follows:

Simulations through Synthetic Data
The identification accuracy of the proposed method was compared with those of the conventional DFT-based and Prony methods in simulations, where the synthetic signal of a transient response was generated as follows: Instant oscillation frequencies and damping factors of the two modes were calculated as where f k (t) and ξ k (t) are nonstationary as time-varying parameters. Synthetic data for modal identification were obtained as y w [n] = y[n] + w[n] with a sampling interval of T s = 1/30 s and an additive white Gaussian noise of w[n]. The segment interval in y w [n] was assigned as Γ = 10 s, which corresponds to two cycles of the first oscillation frequency of 0.2 Hz. The variance of w[n] in each segment was adjusted to meet a given SNR of y w [n]. The duration of y w [n], 0 ≤ n < 3600, in the simulation was 2 min. Here, y w [n] was split into 12 segments of Γ = 10 s in the simulation, where the identification time index was denoted as the end time of each segment. Each of the two combined segments was given as x[l], 0 ≤ l < 600 in the proposed method.
The two modes were identified from each segment by the three methods. By using the estimated modal parametersω k ,ξ k ,Â k , andθ k , fitting values of x[l] in each segment were calculated as the followingx[l]:x [l] = K k=1Â k e −ξ k nT s sin ω k lT s +θ k , 0 ≤ l < 300 The identification accuracies of the three methods can be compared through the following error square ratio (ESR): Two SNRs of 50 and 10 dB were used for the identification of modes in the simulation. Figure 1a,b show y w [n] of 50 dB and its DFT magnitude, respectively. The variation of the modal oscillation frequencies and damping ratios made y w [n] nonstationary, which resulted in distributing its modal spectrum between 0.2 and 0.5 Hz. The two modes were identified 11 times from 12 segments of y w [n].     The ESRs of the three methods were calculated in the simulation and are shown in Figure 3a,b for 50 and 10 dB SNRs, respectively. Estimates of the modal oscillation frequencies and damping factors are plotted in Figures 4 and 5 for 50 and 10 dB SNRs, respectively. The ESR of the proposed method improved significantly compared to that of the conventional DFT-based method. Through simulations on a variety of SNRs, it was confirmed that ESRs of the proposed method were similar to or higher than those of the Prony method when the SNR was lower than 10 dB. This results from performing the curve-fitting in the DFT-based method on a small frequency region around each modal peak frequency, where the SNR was relatively higher than other frequency regions. In the case of the SNR higher than       To compare the identification accuracy of the three methods through simulations, another synthetic oscillation signal of a synchronous mode was generated as the following for y(t): This synchronous mode has an oscillation frequency of f 1 = 5.5 Hz and a damping factor of ξ 1 = 0. Synthetic data for modal identification were obtained as y w [n] = y[n] + w[n], with a sampling interval of T s = 1/120. The white Gaussian noise w[n] was generated to give a 10 dB SNR to y w [n]. For modal identification, y w [n] was split into 20 segments at intervals of Γ = 1 s.
The subsynchronous mode was identified 19 times sequentially from the 20 segments by the three methods, ESRs of which are shown in Figure 6a. Estimates of the modal oscillation frequency and damping factor are plotted in Figure 6b, c. As mentioned in the previous simulation of a 10 dB SNR in Figure 3b, the ESR of the proposed method was also shown to be similar to that of the Prony method in Figure 6a.  To compare the identification accuracy of the three methods through simulations, another synthetic oscillation signal of a synchronous mode was generated as the following for ( ): This synchronous mode has an oscillation frequency of 1 = 5.   Figure 3b, the ESR of the proposed method was also shown to be similar to that of the Prony method in Figure 6a.

Interarea Mode of Central American Power System
The Central American power system was operated on two islands on July 28, 2012 [5]. The larger island included the power systems of Mexico, Guatemala, El Salvador, and Honduras. The smaller island included the power systems of Nicaragua, Costa Rica, and Panama. Undamped active power oscillations between Guatemala and El Salvador started to occur right after the synchronization of the two networks. The Guatemalan wide-area protection scheme (WAPS) issued a trip command to open the interconnection between Guatemala and El Salvador. Data of an active power flow from the Guatemala Este substation to the Moyuta substation in Guatemala was acquired through a PMU with a measurement rate of 30 per second and are plotted in Figure 7a, where modal oscillations start at 9.3 s and subside after the interconnection opening time of 54.3 s. These active power data were filtered with the Butterworth bandpass filter for rejection of a DC component. The filtered data [ ], 0 ≤ < 1350, are shown in Figure 7b, where 0 s in the time axis is 9.3 s in Figure 7a. The DFT magnitude of [ ] is shown in Figure 7c, wherein there is a peak of an interarea mode at around 0.2 Hz.
For modal identification, [ ] from 0 to 40 s was split into four segments at intervals of =

Interarea Mode of Central American Power System
The Central American power system was operated on two islands on July 28, 2012 [5]. The larger island included the power systems of Mexico, Guatemala, El Salvador, and Honduras. The smaller island included the power systems of Nicaragua, Costa Rica, and Panama.  Figure 7b, where 0 s in the time axis is 9.3 s in Figure 7a. The DFT magnitude of y w [n] is shown in Figure 7c, wherein there is a peak of an interarea mode at around 0.2 Hz.
For modal identification, y w [n] from 0 to 40 s was split into four segments at intervals of Γ = 10 s. The interarea mode was identified three times sequentially using the four segments of y w [n] with the three methods. Figure 8 shows identification results of the first segment, where x[l] and its estimated signals are plotted in Figure 8a, and their DFT magnitudes are plotted in Figure 8c. The x[l] of the first two combined segments and its estimated signals are plotted in Figure 8b and their DFT are plotted in Figure 8d. Estimated oscillation frequencies and damping factors of the interarea mode are plotted in Figure 9, which shows that the estimates in the proposed method are similar to those in the Prony method, but the estimates of damping factor in the conventional DFT-based method are significantly different from those in the two other methods, as shown in Figure 9b. This implies that the Prony and proposed methods yielded reliable identification results.

Regional Mode of New England Power System
The New England power system is in the northeastern part of the Eastern Interconnection in USA. Oscillations with a significant MW magnitude were observed in multiple locations of the New England power system on July 20, 2017. A regional mode of 1.13 Hz was engaged in the oscillations, whose data were acquired from Line 2 at substation 2 through a PMU with a measurement rate of 30 per second [6]. PMU frequency data obtained at a 3-minute interval for the oscillations are plotted in Figure 10a, for which oscillation data for modal identification were extracted from 70 to 110 s. These frequency data were filtered with the Butterworth bandpass filter for rejection of a DC component. The filtered oscillation data [ ], 0 ≤ < 1200, are plotted in Figure 10b, where 0 s in the time axis corresponds to 70 s in Figure 10a. The DFT magnitude of [ ] is plotted in Figure 10c, which shows a peak for the regional mode at around 1.13 Hz.
For modal identification, [ ] was split into 16 segments at intervals of = 2.5 s. The regional mode was identified 15 times sequentially using the 16 segments of [ ] in the three methods. [ ] and its estimated signals of the first segment are plotted in Figure 11a, and their DFT magnitudes are plotted in Figure 11c. ̅ [ ] of the first two combined segments and its estimated signals are plotted in Figure 11b and their DFT are plotted in Figure 11d. The estimated oscillation frequencies of the three methods are shown to be similar in Figure 12a. The estimated damping factors in the Prony and proposed methods are similar to each other, but they are significantly different from those in the conventional DFT-based method. This implies that the Prony and proposed methods yielded reliable identification results.

Regional Mode of New England Power System
The New England power system is in the northeastern part of the Eastern Interconnection in USA. Oscillations with a significant MW magnitude were observed in multiple locations of the New England power system on July 20, 2017. A regional mode of 1.13 Hz was engaged in the oscillations, whose data were acquired from Line 2 at substation 2 through a PMU with a measurement rate of 30 per second [6]. PMU frequency data obtained at a 3-min interval for the oscillations are plotted in Figure 10a, for which oscillation data for modal identification were extracted from 70 to 110 s. These frequency data were filtered with the Butterworth bandpass filter for rejection of a DC component. The filtered oscillation data y w [n], 0 ≤ n < 1200, are plotted in Figure 10b, where 0 s in the time axis corresponds to 70 s in Figure 10a. The DFT magnitude of y w [n] is plotted in Figure 10c, which shows a peak for the regional mode at around 1.13 Hz.  For modal identification, y w [n] was split into 16 segments at intervals of Γ = 2.5 s. The regional mode was identified 15 times sequentially using the 16  and its estimated signals of the first segment are plotted in Figure 11a, and their DFT magnitudes are plotted in Figure 11c. x[l] of the first two combined segments and its estimated signals are plotted in Figure 11b and their DFT are plotted in Figure 11d. The estimated oscillation frequencies of the three methods are shown to be similar in Figure 12a. The estimated damping factors in the Prony and proposed methods are similar to each other, but they are significantly different from those in the conventional DFT-based method. This implies that the Prony and proposed methods yielded reliable identification results.   . Estimated modal parameters of the New England power system: (a) oscillation frequency of the regional mode; (b) damping factor of the regional mode.

Conclusions
In this paper, a DFT-based method with an exponential window function has been proposed to identify oscillation modes from segments of modal transient data which involve a long duration with a nonstationary trend. The proposed method overcomes the limited use of the conventional DFTbased method, which cannot be applied efficiently to such segment data. Through simulations on synthetic oscillation data, it was shown that the identification accuracy of the proposed method was significantly improved compared to the conventional DFT-based method. In the simulations, when the SNR was below 10 dB, the identification accuracy of the proposed method was not lower than Figure 12. Estimated modal parameters of the New England power system: (a) oscillation frequency of the regional mode; (b) damping factor of the regional mode.

Conclusions
In this paper, a DFT-based method with an exponential window function has been proposed to identify oscillation modes from segments of modal transient data which involve a long duration with a nonstationary trend. The proposed method overcomes the limited use of the conventional DFT-based method, which cannot be applied efficiently to such segment data. Through simulations on synthetic oscillation data, it was shown that the identification accuracy of the proposed method was significantly improved compared to the conventional DFT-based method. In the simulations, when the SNR was below 10 dB, the identification accuracy of the proposed method was not lower than that of the Prony method. The proposed method can be applied effectively to low-SNR PMU data, which can be acquired on noisy transmission grids, distribution grids, and microgrids. The feasibility of the proposed method was verified by identifying real-world oscillation modes from PMU data, where the proposed and Prony methods yielded similar estimates of modal parameters.
Author Contributions: J.K.H. proposed the core idea, programed the identification methods, and performed the simulations. S.S. and J.S.C. contributed to the collection and analyses of the PMU data. H.-C.K. contributed to the writing of this manuscript and revised the paper.