A Robust Algorithm for Real-Time Phasor and Frequency Estimation under Diverse System Conditions

: This paper presents a comprehensive approach for performing phasor and frequency estimation from voltage and/or current signals of the modern power system. Undesirable components, such as decaying DC, if present in the input signal, are ﬁrst attenuated using a complex-gain ﬁlter. The initial estimates of phasor and frequency are obtained next using the discrete Fourier transform and an improved estimation of signal parameters via rotational invariance technique, respectively. Finally, the accuracy of phasor and frequency estimates are increased based on the identiﬁed system condition. Simulations performed to evaluate the proposed approach conﬁrm that it can do fast and accurate estimation of phasor and frequency under diverse operating conditions, making it ideal for wide-area monitoring, protection, and control applications in power systems.

The proposed research satisfies this need and goes beyond the state-of-the-art in the following ways. First, stray components (such as DDC), if present in the input signal, are removed/attenuated by a filter that is designed based on the concept of complex frequencies.
The initial values of phasor and frequency are then estimated in parallel from the filtered signal using DFT and an improved version of estimation of signal parameters via 1.
It efficiently suppresses undesirable components present in the input signal, and can therefore be applied to both voltage and current signals.

2.
It provides fast, accurate, and robust phasor and frequency estimates under diverse operating conditions (faults, frequency ramps, power swings, harmonics). 3. It detects abnormal system conditions quickly (in less than one-cycle).

4.
It satisfies the requirements of relevant IEEE and IEC Standards (this paper focuses on P-class PMUs which have more stringent speed requirements) [29,30].

Proposed Filter for DDC Rejection
Considering the most general setting, let the input signal be described by where A 0 denotes a DC offset; K denotes the total number of sinusoidal components; A k and θ k denote the magnitude and phase angle of the k-th sinusoid. The angular frequency of the k-th sinusoid is denoted by ω k and is equal to 2π f k , where f k is the frequency of the k-th sinusoid; the nominal frequency and angular frequency are denoted by f nom and ω nom , respectively. D 0 and τ are the magnitude and time constant of the DDC component and are called DDC parameters. The above signal can be written in discretized form as A k cos(ω k n∆t + θ k ) + D 0 E n (2) in which ∆t is the sampling interval and is equal to 1/ f s (where f s is the sampling frequency), n denotes the sample number, and E ∆ = e −∆t/τ , called the damping factor, depends on the DDC time constant [4]. It is difficult to predetermine the DDC parameters as they depend on fault characteristics, such as, fault inception time/angle, fault location, and fault resistance [4][5][6]31]. However, their presence (usually in the current signals) can cause up to 15% error in the estimated phasor [31]. This can result in a maximum overreach of 115% if the first zone of a distance relay is set at 100% of the line impedance. Thus, DDC suppression is important for safe and reliable operation of the smart grid.
The goals of the proposed filter design are: (i) removing/attenuating the DDC for different values of its time constant, (ii) preventing variation in parameters of the fundamental component during the removal process, and (iii) attaining fast response time (or equivalently low-order characteristic). The frequencies present in (1) are shown in Figure 1, where the frequencies of the DC offset, DDC, and sinusoidal components are located at the origin, negative side of the real axis, and imaginary axis of the s-plane, respectively. The DDC time constant (τ) varies due to fault characteristics (such as fault impedance and fault location) between 0.5 and 5 cycles [31,32]; hence, a corresponding lower bound (τ l ) and upper bound (τ u ) are specified. Since the DDC frequency is equal to −1/τ in the s-plane, the zeros of the proposed filter should lie between [−1/τ l , −1/τ u ], which is called DDC frequency band (see Figure 1), to filter out the DDC. However, the filter's zeros will cause phase shifting of the fundamental component. To compensate for this phase-shift, filter's poles can be added (as was done in [6]). However, doing so will result in: (i) neutralization of the effect of DDC attenuation by the filter's zeros, (ii) additional error in phasor estimation due to the vicinity of the added poles to the fundamental frequency, and (iii) oscillatory behavior in the estimated magnitude and phase angle. To solve this problem, in this paper, we propose the design of a complex gain filter. The transfer function of the proposed filter is given by where the total number of zeros (or equivalently the filter order) is denoted by R which plays a key role in the response time of the proposed filter. The filter gain is expressed by K F which is complex (i.e., K F ∆ = |K F |e jγ F ). z m s are inserted in the DDC frequency band with the aim of reaching minimum multiplication of their distance to any frequency in this band. This can be formulated as the following optimization problem: Note that this optimization problem ensures that the proposed filter suppresses DDC components regardless of their time-constant values, without any re-tuning.
Next, to compensate for the phase shift due to the added zeros, the phase-angle of the filter gain γ F is adjusted as follows (see Figure 1): Finally, the magnitude of the filter gain |K F | is determined with the aim of reaching unit total filter gain at nominal frequency f nom . This is given by Since the proposed filter is designed for ω nom , its performance may deteriorate at offnominal frequencies. Under such circumstances, the filter's parameters can be readjusted as described in Section 6.2. As the number of filter's zeros increases, the filter can suppress DDC components more effectively, at the cost of increased response time. This trade-off between the proposed-filter's capacity for DDC attenuation and its response time is investigated in Section 7.1.

Phasor Estimation at Off-Nominal Frequencies
DFT is one of the most widely used techniques for phasor estimation due to its simplicity in implementation, low computational burden, harmonic filtering, and noise alleviation [1]. By applying DFT on the filtered signal, denoted by subscript f , we get where r is the number of the first sample in the data string, N denotes the number of samples per cycle, andX 1 refers to estimated fundamental phasor at nominal frequency. Note that the DC offset and harmonic components, present in (2), are filtered out by the DFT, while the DDC has already been attenuated by the proposed filter (see Section 2). However, when the frequency deviates from its nominal value f nom , the phasor estimated using DFT can be erroneous. This requires additional correction, as illustrated below.
Consider an off-nominal frequency condition in which the angular frequency is different from its nominal value ω nom ; i.e., let the actual angular frequency be ω 1 . The DFT of the filtered input signal in this situation can be written aŝ whereX 1 is the estimated phasor by DFT at actual frequency ω 1 . After performing trigonometric simplifications on (7) and (8), we get X 1 = PX 1 e jr(ω 1 −ω nom )∆t + QX * 1 e −jr(ω 1 +ω nom )∆t (9) where P and Q are correction coefficients that are independent of 'r' and dependent on N and frequency deviation ∆ f . Variations of P and Q for different values of N and ∆ f are shown in Tables A1 and A2 in Appendix A.1. In a typical power system, the variations in the frequency rarely exceed 0.5 Hz in either direction from its nominal value [33]. As shown in Tables A1 and Appendix A.1, Q is negligibly small in comparison to P in the range of −0.5 Hz and +0.5 Hz. The impact of Q on the accuracy of the estimated phasor by the proposed algorithm is investigated in more detail in Appendix A.2. From the conducted investigation, it is realized that the impact of Q on the accuracy of the estimated phasor by the proposed algorithm is small (in comparison to the total vector error (TVE) limit of 1% [29,30]). Consequently, (9) is simplified and rewritten as follows: SinceṔ depends on the actual frequency, it is necessary to estimate the frequency precisely and rapidly for correcting the phasor estimate. The frequency estimation algorithm described in Section 4 can not only improve accuracy of the estimated phasor at off-nominal frequency conditions (see Section 6.2), but also estimate the frequency quickly and accurately for other power system applications.

Proposed Frequency Estimation Algorithm
ESPRIT is a well-known subspace-based technique for frequency estimation due to its strong immunity against noise [28]. However, due to its high computational complexity [25], ESPRIT is not suitable for power system applications that require fast estimation of frequency. Hence, an improved version of the ESPRIT, called iESPRIT, is proposed in this paper. The input signal for frequency estimation is the output of the filter that was proposed in Section 2, and can be expressed as (13) where w f [n] denotes an additive noise component with the power of σ n 2 . Note that the filter was designed in such a way that it kept the parameters of the fundamental component unchanged. Equation (13) can be equivalently written in the form of Y ∆ = 2K + 1 complex exponentials distorted with noise, called the harmonic model, as follows: (14) in For simplicity, we denote the Y α k s and ω k s for k = −K, · · · , 0, · · · , K with α 1 , α 2 , · · · , α Y and ω 1 , ω 2 , · · · , ω Y , respectively, and let ω k where ω k is thenormalized angular frequency of the k-th exponential component. Then, by using M consecutive samples of x f [n], an M-length data vector x f (n) is defined as follows: .
where the Y columns of matrix H are length-M signal frequency vectors described by for k = 1, · · · , Y. The matrix H is simplified using (16) as In addition, matrix Φ is the diagonal matrix of phase shifts and can be expressed as Since the frequencies of the complex exponentials ω k s completely describe the matrix Φ, the estimated frequencies can be obtained by finding Φ.
Let us consider two overlapping sub-matrices within the length-M window vector. In this regard, Then, two new (M − 1) × Y matrices which are called H 1 and H 2 are defined as follows: Using (15) and (19), we have: As seen in (19), the matrices H 1 and H 2 are related as Now, a G × M matrix, X f , called the data matrix, can be obtained from an L-length of input data in which L ∆ = G + M − 1 and G is the total number of data vectors used to construct X f . In conventional ESPRIT, Φ is calculated from the singular value decomposition (SVD) of X f , as shown below [28] X where L G×G and U M×M are the left and right singular vectors (both of them are unitary), and Σ G×M is a diagonal matrix containing singular values ordered in descending magnitude. However, SVD of X f will have a complexity of O(GM 2 + MG 2 + M 3 ) [34], which makes ESPRIT too slow for frequency estimation in P-class PMUs [1,29,30]. To overcome this problem, the proposed iESPRIT estimates an autocorrelation matrix from X f as shown belowR whereR x can be expressed in terms of its eigenvalues and eigenvectors using eigendecomposition (ED) asR where Λ is an M × M diagonal matrix containing the eigenvalues of theR x in descending order. The columns of Q are also the corresponding eigenvectors. It is proved in Appendix B that the eigenvectors ofR x (columns of Q) are equal to the columns of U which is used in ESPRIT (see (23)). Note that the complexity of ED ofR x is O(M 3 ) [34], which is much lower than the complexity of SVD [35]. The Y largest diagonal elements of Λ (i.e., λ 1 , λ 2 , · · · , λ Y ) are equal to M|α k | 2 + σ n 2 (for k = 1, 2, · · · , Y), while the remaining M − Y + 1 diagonal elements of Λ (i.e., λ Y+1 , · · · , λ M ) have similar values as σ n 2 [28]. As a result, the total number of complex exponentials Y can be obtained using the criterion of the relative slope variation between theR x 's eigenvalues as shown in Figure 2. Note that the fundamental frequency corresponds to the largest singular value in Λ. In order to find Φ, the matrix Q is split into two sub-matrices for separation of the signal and noise subspaces where Q S is an M × Y sub-matrix that corresponds to the Y largest magnitudes of the eigenvalues (see (25)). Since the elements of the matrix H defined in (17) completely describe the frequencies of exponential components existing in the input signal, H similar to Q S must also lie in the same subspace, called the signal subspace. As a result, an invertible transformation W exists that maps Q S into H as Note that it is not necessary to calculate the transformation W, and it is only formulated as a mapping between the matrices H and Q S within the signal subspace. The matrix Q S can also be split into two smaller M − 1-dimensional subspaces as was done with the matrix H where both Q 1 and Q 2 are (M − 1) × Y sub-matrices. Since H 1 and H 2 correspond to the same subspaces as Q 1 and Q 2 , (27) must also hold for these subspaces, as shown below: As seen in (21), the rotational matrix Φ exists that relates (rotates) H 1 to H 2 . In the same way, a rotational matrix (denoted by Ψ in this paper) must also exist that relates (rotates) Q 1 to Q 2 as Note that the matrices Q 1 and Q 2 are known from the ED ofR x . Therefore, from (30), the rotational matrix Ψ can be estimated using LS aŝ Finally, using (21), (29), and (30), we get Now, as the diagonal entries ofΦ (φ k = e jω k for k = 1, 2, · · · , Y) are simply the eigenvalues ofΨ, the estimated frequencies can be obtained aŝ Recall ω k = 2π f k / f s ; therefore, the estimated frequencies arê

Fast Detection of System Operating Conditions
The input-signal parameters (magnitude, phase angle, frequency) may abruptly change/oscillate due to power system events, such as faults. Termed abnormal conditions in this paper, these sudden changes can quickly be detected usingR x . Let the estimated autocorrelation matrix corresponding to the previous data string be denoted byR x old . The new autocorrelation matrixR x new is estimated using (24) when the most recent data string is obtained. We can now create an M × M matrix ∆R x as shown below: The Frobenius norm of ∆R x , ||∆R x || F , is given by A threshold ζ can now be defined for ||∆R x || F to determine the system's current operating condition as follows: 1.
If ||∆R x || F ≤ ζ, it implies that the system operates under normal conditions and the signal parameters have varied slightly. In such conditions, the estimated phasors from previous and current data strings can be used for accurate estimation of frequency (see Section 6.1).

2.
If ||∆R x || F > ζ, it implies that the system is faced with an abnormal condition resulting in the signal parameters differing significantly from the previous data string. In such conditions, the proposed iESPRIT must be used for precise estimation of frequency, and subsequent correction of the estimated phasor (see Section 6.2).
The impact of ζ on the accuracy and computational burden of the proposed approach is investigated in detail in Section 8.6.

Accuracy Enhancement of Estimated Phasor and Frequency under Detected Conditions
The main goal of this section is to improve the accuracy of the estimated phasor (Section 3) and frequency (Section 4) based on the system operating condition (normal/abnormal) identified in Section 5.

Implementation of the Phasor Estimation Algorithm for Frequency Estimation under Normal Conditions
Let the estimated fundamental phasor obtained from the previous N-length data string with sample numbers {r, r + 1, · · · , r + N − 1} be denoted byX 1old , which is calculated as Similarly,X 1new is obtained from the current N-length data string with sample numbers {r + 1, r + 2, · · · , r + N} at ω 1new = 2π f 1new aŝ It is proved in Appendix C that, under normal conditions, the parameterṔ (defined in (12)) can be obtained as follows.Ṕ Now, asṔ only depends on N, r, and ∆ω 1 , for a pre-determined value of N,Ṕ can be calculated offline for different values of r and ∆ω 1 in the form of a look-up table. Then, in real-time, the ∆ω 1 can be obtained with high accuracy from its correspondingṔ using the look-up table. Finally, the estimated frequency can be updated during normal conditions aŝ The steps of the proposed frequency estimation algorithm are summarized in Algorithm 1.

Implementation of the Proposed Frequency Estimation Algorithm for Phasor Estimation under Abnormal Conditions
The first valid estimated frequency realized using Algorithm 1 is obtained after receiving the first L-length of input data; let this frequency be denoted byf 1,L . It will be shown in Section 7.4 that the proposed iESPRIT is able to estimate the frequency precisely even when L is lower than one-cycle-length of input data (i.e., L < N), which is a significant achievement especially for protection applications. Therefore, it is possible to correct the estimated phasors by usingṔ which is dependent on the accurate frequencyω 1 (see (9), (10) and (12)). To do so, we calculate the median (Mdn) of a set F that contains the previous and current estimated frequencies (i.e., F = {f 1,L ,f 1,L+1 , · · · ,f 1,N }) as shown in (41).
Employing the median brings additional robustness against noise. The value of the actual angular frequencyω 1 , obtained from (41), can be used to correct the estimated phasors during abnormal conditions. Moreover,ω 1 can also be used to improve the performance of the proposed filter design (see Section 2) by replacing ω nom withω 1 in (5) and (6). In this regard, the filter's parameters can be calculated offline, for the pre-chosen filter order (see Section 7.1), for different values of ω 1 using (5) and (6), and stored in the form of a look-up table. In real time, the filter's parameters can be adjusted quickly usinĝ ω 1 and the look-up table.
Finally, note that the proposed frequency estimator itself provides strong immunity against noise and can be used for various power system applications that require accurate and fast frequency measurements, such as under-frequency relaying, power system stabilization/restoration, protection against loss of synchronism, and rapid synthetic inertia control.
Flowchart of the proposed approach for both phasor and frequency estimation is shown in Figure 3 in which normal and abnormal system operating conditions have been distinguished with green dashed and red solid arrows, respectively.

Choosing Optimal Value for the Filter Order R
To determine the filter order R of the complex gain filter (see Section 2), 100 synthetic test signals containing both fundamental and DDC components are generated. In the test signals, the D 0 and τ (DDC parameters) are randomly selected from a range of [0, 1 p.u.] and [0.5, 5] (cycles), respectively [5]. Seven different values {3, 4, 6, 8, 10, 15, 20} are considered for R. The input signals are passed through the proposed filter and DFT is then performed for estimation of the fundamental phasor. For each filter order, the maximum values of TVE [29,30] and settling time t S [4], obtained over the 100 test signals are given in Table 1. As can be seen from Table 1, by increasing R, the DDC is suppressed more effectively resulting in reduction of the maximum TVE at the cost of increasing the maximum t S ; i.e., there is a trade-off in choosing the optimal value of the filter order. Table 1 indicates that an acceptable value for the order of the proposed filter would be six since this provides a reasonable balance between DDC rejection (significantly lower maximum TVE than 1%) and response time.

Choosing Optimal Value for N
The impact of N (or equivalently f s ) on the accuracy of the estimated phasor by the proposed algorithm is investigated in this section. In this regard, a synthetic test signal containing a fundamental sinusoidal component that is distorted by DDC, harmonics, and noise is used as an input to the proposed algorithm. Six different values {16, 32, 128, 512, 1024, 16,384} are considered for N. For each value, the fundamental phasor is estimated by the proposed algorithm, and the maximum TVE is calculated, which is given in Table 2. As can be seen in Table 2, the phasor estimation accuracy of the proposed algorithm converges to its limit value and is negligibly affected by increasing the number of samples per cycle, N (or equivalently the sampling frequency, f s ). On the other hand, increasing the sampling frequency causes (i) difficulties in hardware implementation of the phasor and frequency estimation algorithms in numerical relays, and (ii) increased computational burden, which is a major concern for the protection application. Based on Table 2, an appropriate value of N was found to be 32.

Choosing Optimal Value for M
The value of M plays a key role in computational complexity and robustness of the proposed iESPRIT against noise. As such, seven different values {3, 4, 5, 6, 8, 10, 15} are considered for M. A synthetic test signal containing a sinusoidal component (with known frequency) which is distorted by noise is used as an input to the proposed iESPRIT. Note that the value of L (length of data string in (22)) is kept constant at N. The maximum absolute frequency estimation error (denoted by Max. |FE|) [30] and the execution time taken for estimation are presented in Table 3 for each of the seven values of M. Note that all of the simulations done in this paper have been carried out on a system with Core 2 Duo CPU (2.67 GHz) processor and 4 GB RAM. It can be concluded from Table 3 that more robustness against noise and higher precision are achieved by choosing greater values for M, at the cost of increasing the execution time due to increasing complexity of calculations. Therefore, there is a tradeoff between the precision of estimation and computational burden. From Table 3, an appropriate value of M is found to be four.

Choosing Optimal Value for L
The value for L (length of data string) plays a key role in convergence speed of the proposed iESPRIT. In this regard, six different values {N/4, N/2, 3N/4, N, 3N/2, 2N} are considered for L. Likewise, as we did in Section 7.3, a synthetic signal is used while M is kept constant at four.
The Max |FE| and settling time t S are given in Table 4 for the six values of L. Considering optimal accuracy and convergence speed of the proposed iESPRIT, an appropriate value of L is found to be N/2. Note that the empirically obtained values of the parameters identified in Section 7 are kept constant for all the simulations performed in Section 8.

Simulation Results
In the following sub-sections, several static and dynamic test signals that comply with the IEEE and IEC Standards [29,30], as well as power system signals, are used to investigate the performance of the proposed approach. The above-mentioned standards only consider changes in the input signal parameters (i.e., magnitude, phase angle, and frequency). However, in this paper, the performance of the proposed algorithm is also evaluated in the presence of stray components, such as DDC and noise. Four recently published phasor and/or frequency estimation methods, applicable to P-class PMUs, have been selected for comparison purposes: (i) A phasor estimation method presented in [5], which is based on DFT and down sampling (called "DS-DFT"), (ii) a phasor and frequency estimation method presented in [13], which is based on LS and amplitude modulation modeling of the input signal (called "AM"), (iii) a phasor and frequency estimation method presented in [13], which is based on LS and phase modulation modeling of the input signal (called "PM"), and finally (iv) a subspace-based frequency estimation method inspired by the multiple signal classification (MUSIC) technique presented in [24] (called "A 2 MUSIC").
The nominal frequency f nom , sampling frequency f s , and ζ are set to 60 Hz, 1920 Hz, and 10 −4 , respectively.

Performance Evaluation in the Presence of DC Offset, DDC, Harmonics, and Noise in Both Nominal and Off-Nominal Frequency Conditions
The following test signal is used to investigate the performance of the proposed approach for phasor and frequency estimation in the presence of DC offset, DDC, harmonics, and noise, in both nominal and off-nominal frequency conditions where the total number of sinusoids K is set equal to 50 [29,30]. Five-hundred test signals are generated by selecting different values for A 0 , f 1 (varied from 58 to 62 Hz [29]), D 0 (varied from 0 to 1 p.u.), τ (varied from 0.5 to 5 cycles), and signal-to-noise ratio (SNR) (varied from 10 to 80 dB) [24]. The fundamental phasor and frequency are estimated by the proposed approach and other understudied algorithms, and the mean, root mean square (RMS), and maximum values of maximum of TVE and |FE| are calculated. From the results shown in Table 5, it is clear that the proposed approach gives higher accuracy for both phasor (Max. TVE) and frequency (Max. |FE|) estimation in comparison to the other methods.

Performance Evaluation for Dynamic Frequency Ramp Test
In this test, for the same test signals generated in Section 8.1, f 1 is ramped to f 1 + 5 (Hz) from 0 to 5 (s) at a rate of 1 Hz/s [3,29,30]. The mean, RMS, and maximum values of the two indices, namely, Max. TVE and Max. |FE|, for the 500 dynamic test signals are given in Table 6.
As seen in Table 6, the proposed approach gives better performance for both phasor and frequency estimation under frequency ramps in comparison to the other methods. Apart from the high estimation accuracy, the proposed approach also has faster response time under dynamic conditions which will be explored in more details in Section 8.4.

Modulation Test
Transmission line faults/outages and bulk load switching often cause rotor-angle oscillations which may lead to power swings and subsequently mal-operation of the protective relays [32]. Generally, power swings are mathematically modeled as a modulated sinusoidal signal [13,29,30,32] x 2 (t) = X m (1 + k x cos (ω m t)) × cos (ω nom t + k a cos (ω m t − π)) where k x and k a denote amplitude and phase modulation coefficients, respectively, both of which are set to 0.1 [29,30]. In this test, the frequency of the input signal remains constant at 60 Hz, but the magnitude and phase-angle oscillate with an angular frequency of ω m = 2π f m where f m changes from 0.1 Hz to 2 Hz in steps of 0.1 Hz as recommended by the Standards [29,30] for P-class PMUs. The maximum values of TVE and |FE| are given in Table 7 for different algorithms. As can be observed from Table 7, the proposed approach gives better frequency estimation accuracy compared to the other algorithms. Since the algorithms presented in [13] are based on amplitude modulation (AM) and phase modulation (PM) modeling of the input signals, they have the highest phasor estimation accuracy (their Max. TVE is lowest). However, the proposed phasor estimation approach also satisfies the modulation test requirements for P-class PMUs specified in the Standards [29,30].

Step Change Test
Abrupt changes in the impedance and network configuration during faults may cause current/voltage transients which are mathematically expressed as [13,29,30,32] x 3 (t) = X m (1 + k x u(t)) × cos (ω nom t + k a u(t)) where k x and k a are the magnitudes of the unit step function u(t) and set to 0.1 and π/18 rad, respectively, as recommended in the Standards [29,30]. Two well-known indices, namely, response time and maximum over/undershoot, which are defined in [30], are used to evaluate the performance of the different algorithms from the viewpoints of speed and accuracy. The results are given in Table 8 in which the response time and maximum over/undershoot are denoted by Rsp. and Max O/U.S., respectively. Based on the results obtained in Table 8, it is clear that the proposed approach has superior performance in terms of speed and accuracy, and its performance indices can comply with the requirements for P-class PMUs mentioned in the Standards [29,30].

Performance Evaluation Using Power System Signals
Power system signals of the IEEE 9-bus test system [13] are used to evaluate the performance of the proposed approach for both phasor and frequency estimation when severe oscillations in magnitude, phase angle, and frequency of the input signals take place. Two different scenarios are considered [13]: Scenario 1: Outage of line 7-8 at t = 0 s. Voltage of bus 6 is used as the input signal. Scenario 2: G 2 (connected to bus 2) trips at t = 50 ms. Current of line 6-4 is used as the input signal.
The maximum values of TVE and |FE| are given in Table 9 for the two scenarios. As can be observed from Table 9, the proposed approach provides higher accuracy compared to the other methods. It is also worth mentioning that the abnormal events in both the scenarios were detected by the proposed approach in less than one cycle.

Computational Burden and Sensitivity Analysis
The threshold ζ has a significant impact on the accuracy and the computational complexity of the proposed approach (see Section 5). Hence, a sensitivity analysis is performed to investigate the effect of ζ on the accuracy of the estimated phasors and frequencies and the computational burden of the proposed approach. In this regard, Max. TVE, Max. |FE|, and average execution time (taken for each estimation) are given in Table 10 for different values of ζ for each of the test signals used in this section. As can be observed from Table 10, higher accuracy in estimation of phasor and frequency is achieved by choosing smaller values for ζ, at the cost of increasing the execution time (or equivalently computational burden/complexity) i.e., there is a trade-off in choosing the optimal value of the threshold ζ. Table 10 indicates that an acceptable value for ζ would be 10 −4 since this provides a reasonable balance between the estimation accuracy and computational burden. The mean execution time of the proposed approach (over all of the test signals in this section) for each synchrophasor and frequency estimation by the adoption of 10 −4 for ζ was around 0.34214 (ms), which makes it suitable for realtime applications.

Noise Propagation through the Proposed Algorithm
In this section, the following test signal is used to investigate the propagation of noise through the proposed algorithm: Cramér Rao bounds (CRBs) that give lower bounds on the variances of estimated signal parameters (i.e., magnitude, phase angle, and frequency) are obtained for the abovementioned test signal (45) as follows: [36].
In order to investigate the propagation of noise through the proposed algorithm, 100 test signals are generated by choosing different values for f 1 between 58 and 62 (Hz) [29,30]. The variance of the noise in the sinusoidal signal σ n 2 is set to be 0.01 2 [36]. For each of the generated test signals, the signal parameters are estimated by the proposed algorithm. The mean square error (MSE) [37] as well as the CRBs are given in Table 11. As can be seen from Table 11, the values of MSE for all of the estimated parameters are very close to their corresponding CRBs.

Conclusions
This paper proposes a fast, accurate, robust, and computationally-efficient approach for estimating phasor and frequency that can be applied to both voltage and current signals without making simplifying assumptions regarding the components present in the signal. The first step of the proposed methodology involves the removal of the DDC that may be present in the input signal using optimal filtering. Care is taken to ensure that the parameters of the fundamental component are left unaltered. In the next step, initial phasor and frequency estimates are obtained in parallel using DFT and the proposed iESPRIT, respectively. Finally, to further improve accuracy, the initial estimates of phasor and frequency are updated, based on the identified system condition.
The simulation results confirm that the proposed approach has high accuracy, fast response time, low computational complexity, and is robust against grid disturbances and off-nominal frequency conditions. Furthermore, the proposed approach fulfills the requirements recommended by the IEEE and IEC standards for P-class PMUs in both steady-state and dynamic conditions. This makes the proposed approach suitable for a variety of power system monitoring, protection, and control applications that need fast and accurate estimates of voltage and current phasors and system frequency.
To provide more robustness against stray components such as harmonics and interharmonics, particularly for measurement applications, the combination of more recent versions of DFT such as interpolated DFT (IpDFT) with iESPRIT will be explored in the future.  It is clear from (A7), (23), and the (25) that the squared magnitudes of the data matrix singular values are equal to the eigenvalues of theR x scaled by a factor of G (i.e., Σ 2 = GΛ), and the columns of the matrix U are exactly equal to the columns of Q (i.e., Q = U).

Appendix C. Proof ofṔ Calculation Using (39)
The DFT of the previous data string in (37) is split aŝ Similarly, (38) decomposes as follows for the current string: According to (7), (8), (10) and (12), the following equation exists that relates J to Z: By substituting (A8) and (A9) into (A10), we havê Now, note the following assumptions (which are not violated under normal operation of power systems): • Deviation of frequency from its nominal value is very small (hence, e −jNω 1new ∆t 1).

•
The input signal is approximately periodic (hence, x[r] x[r + N]).
Based on these assumptions, (A11) can be simplified to obtain (39).