Computationally Efficient Synchrophasor Estimation: Delayed In-Quadrature Interpolated DFT

This article proposes a synchrophasor estimation (SE) algorithm that leverages the use of a delayed in-quadrature complex signal to mitigate the self-interference of the fundamental tone. The estimator, which uses a three-point IpDFT combined with a three-cycle Hanning window, incorporates a new detection mechanism to iteratively estimate and remove the effects caused by interfering tones within the out-of-band interference (OOBI) range. The main feature of the method is its ability to detect interfering tones with an amplitude lower than that adopted by the IEC/IEEE Std. 60255-118, this detection being notably challenging. Furthermore, it simultaneously satisfies all the accuracy requirements for the P and M phasor measurement unit (PMU) performance classes while offering a reduction in the total computational cost compared with other state-of-the-art techniques.


I. INTRODUCTION
S INCE their introduction in the 1980s [1], Phasor Measure- ment Units (PMUs) have evolved and are currently present in the power grids of all developed countries [2], used successfully in many applications, such as wide-area monitoring protection and control, model validation and state estimation [3].To ensure the inter-operability and compatibility of PMUs from different manufacturers, an industry standard, whose latest installment corresponds to the IEC/IEEE Std.60255-118 [4], has been in continuous development since 1992 [5].The standard [4] defines a set of tests representative of the operating conditions of the power system, and the corresponding performance requirements that all PMUs must satisfy based on their class and reporting rate.Two different performance classes are defined in [4]: class P, for applications requiring a fast response, and class M, for those requiring high accuracy and rejection of interharmonics or subharmonics.
Over the years many synchrophasor estimation (SE) techniques have been explored in the literature such as wavelets [6], Prony [7], Taylor-Fourier [8]- [10], Shanks [11], Kalman filtering [12]- [14] and adaptive filters [15], among others.Most commercial PMUs employ DFT-based SE algorithms [16], which can produce accurate estimates by simply processing a few DFT bins [17].However, two main limitations compromise the performance of DFT-based techniques: aliasing and spectral leakage [18].Although the first can be easily solved by increasing the sampling rate or adopting an anti-aliasing filter, spectral leakage, which includes long-and short-range interferences, requires more elaborate solutions [18].Long-range leakage refers to the mutual interference caused by all the tones that make up the signal spectrum, while short-range leakage, also known as scalloping or picket fence effect, is the error caused by the displacement of the maximum bin.Long-range leakage can be mitigated by windowing [19], while interpolation of DFT bins [20] addresses short-range leakage.Combining both windowing and interpolation, as proposed in [21], represents the foundation of modern interpolated DFT (IpDFT) techniques.To satisfy the underlying assumption that the signal parameters are fixed within the observation window, and to comply with the time and latency requirements set by [4], intervals of just a few nominal cycles are generally selected [22]- [24].In turn, this causes a coarse frequency resolution which, in the case of a real-valued power system signal whose main tone nominal values fall close to DC, results in the proximity between the positive and negative images.Their mutual interaction, also known as self-interference, constitutes the main source of error in the estimation process [24].
Different approaches have been explored to deal with negative frequency infiltration [17], [22]- [28].[25] introduced a multi-point weighted IpDFT that reduces the effects of long-range spectral leakage.[17] employed a method based on a novel three-point weighted IpDFT combined with the use of maximum sidelobe decay (MSD) windows, which showed a great rejection of interference from the negative image.[27] developed new cosine windows, called Maximum Image interference Rejection with Rapid Sidelobe Decay rate (MIR-RSD) windows, which showed a high rejection of selfinterference and long-range leakage from other narrowband disturbances.[22] presented an algorithm (e-IpDFT) that iteratively approximated and compensated for the effects of selfinterference by taking advantage of the symmetry of the DFT spectrum with respect to the DC component.The method was further extended in [23] to further compensate for an additional generic interfering tone.The resulting method, the so-called i-IpDFT, was proven to meet the requirements for the P and M classes.Both e-IpDFT and i-IpDFT were proved to be compu-tationally efficient and experimentally deployed and validated on industry-grade field-programmable gate arrays (FPGAs) [29], [30].The same idea of estimating and compensating for the effects of self-interference and long-range leakage from additional tones presented in [23] is also used in [28].While [23] uses an iterative formulation, [28] derives analytical expressions to achieve faster convergence.However, [28] does not account for subharmonic or interharmonic components.
A different approach to mitigate negative image infiltration is to generate a complex signal from the real one.HT-IpDFT [24] employs a Hilbert filter to approximate the analytical signal and suppress the negative spectrum.Although the method reduces the computational complexity of the iterative process compared to [23], it does not meet the combined requirements of classes P and M for harmonic distortion (1%) and phase step tests in [4].Another suitable complex signal to cancel the negative image is the one defined by the signal itself and its imaginary in-quadrature signal.Many quadrature signal generation (QSG) techniques exist and have been used in other applications, such as single-phase phase-locked loops (PLLs) [31].In [26] an in-quadrature signal is obtained by applying the Clarke transform to a three-phase signal generated by buffering and shifting a noise-filtered version of the measured single-phase signal.The complex signal is further enhanced with a weighted least-squares Taylor-Fourier (WLS-TF) filter to match the magnitudes of both in-quadrature components and is used with a Blackman window filter.
This paper proposes a SE algorithm that takes advantage of the use of a delayed in-quadrature complex signal generated in the time domain to mitigate the self-interference of the fundamental tone.The estimator, named Time-Delay IpDFT (TD-IpDFT), uses a three-point IpDFT combined with a Hanning window and incorporates a new detection mechanism to iteratively estimate and remove the effects caused by interfering tones equal to or greater than 5% within the out-ofband interference (OOBI) range.It simultaneously satisfies all accuracy requirements for the P and M classes, while offering a reduction in the total computational cost compared to the i-IpDFT [23], [30].
The paper is structured as follows.Section II reviews the fundamentals of IpDFT-based SE techniques.Section III describes the theoretical basis of how a complex signal consisting of in-quadrature components allows the suppression of the negative image spectrum.Section IV presents the structure, formulation, and adjustment of the parameters of the TD-IpDFT SE algorithm.Section V presents the complete characterization of the algorithm for P-and M-class PMUs indicated by [4].Section VI discusses and compares the results with those obtained by other state-of-the-art SE methods.Finally, VII summarizes the main findings and gives closing remarks.

II. SE IPDFT TECHNIQUES: FUNDAMENTALS
We refer to a formulation that considers a Hanning window function, which offers a good compromise between sidelobe decay and mainlobe width [21], and a three-point DFT interpolation scheme, as it reduces the effects of long-range leakage [25].This configuration, already adopted in [23], [24], is also used for the proposed TD-IpDFT.
A. Three-point (3p) IpDFT based on the Hanning window Consider a discrete single-tone steady-state signal x(n): where T s represents the sampling period and A 0 , f 0 and ϕ 0 the fundamental amplitude, frequency and initial phase, respectively.If a window of N consecutive samples is selected, the DFT spectrum of the signal X(k) is given by (2), where w(n) is the discrete windowing function and According to the convolution theorem, a conjugate spectrum is obtained, where each spectral component X(k) is the result of the combined contributions of a positive and negative image of the fundamental tone, each being the scaled and rotated versions of the DFT of the window function W (k) shifted at ±f 0 /∆ f1 : where V 0 = A 0 e jϕ0 .The distance between positive and negative images depends on the frequency resolution ∆ f , which is the reciprocal of the length of the selected window T = N T s .For the Hanning window, W (k) corresponds to: where W R (k) is the DFT of the rectangular window, also known as the Dirichlet kernel: Given a set of estimates of the signal parameters ( f , Â, φ) and knowing the windowing function adopted, the spectral contributions of positive and negative images can be estimated by: For the general case of incoherent sampling2 , the fundamental tone of the signal falls between subsequent bins.Thus, the signal frequency can be expressed as where k m is the index of the highest bin and δ a fractional correction term that, for the Hanning window (δ H ), can be analytically calculated by interpolating the three highest DFT bins as [25]: where With δ H determined, the amplitude and phase of the fundamental tone are given by:

III. NEGATIVE FUNDAMENTAL IMAGE SUPPRESSION
A. Delayed in-quadrature technique Consider a complex discrete single-tone steady-state signal ȳ(n), with components y(n) and y(n − d θ ), θ being the resulting phase shift that delay d θ (d θ ∈ N) introduces at the corresponding frequency f 0 : , n ∈ N (10) where ω 0 = 2πf 0 is the angular frequency3 and θ = ω 0 d θ T s .If (10) is expressed in terms of complex exponentials and the positive and negative frequency components are grouped, the following expression is obtained: where σ + and σ -are the complex positive and negative delay gains and correspond to: σ -= 1 + e j(π/2+θ) (12b) The complex signal ȳ(n) can be seen as a 'filtered' version of y(n), where the amplitude and phase of each positive and negative frequency component are altered respectively by (12a) and (12b) based on the relative angular shift θ.Let ... be the round-to-the-nearest integer function, rounding the equidistant values away from zero to the integer with a larger magnitude.Fig 1 shows the frequency response of a 'filtered' delayed in-quadrature complex signal ȳ(n) generated considering a frequency f from a signal y(n) with normalized frequency f 0 [pu] = f 0 /f .The magnitude and phase response for each positive and negative frequency component are given, respectively, by |σ + |, ∠σ + and |σ -|, ∠σ -with d θ = f s /(4f ) .The frequency response is 4 pu periodic and is centered on the frequency considered for in-quadrature generation f (1 pu).If there is no discrepancy between the signal frequency (f 0 ) and that used for the generation of the in-quadrature component, i.e. (f 0 = f → f 0 [pu] = 1), a perfect in-quadrature signal is obtained, resulting in cancelation of the negative image (σ -= 0) and doubling of the positive one (σ + = 2).Otherwise, if the signal frequency falls in its vicinity (f 0 f → f 0 [pu] 1), still a significant mitigation of the negative image occurs, but the effects can only be quantified by (12) once the signal frequency has been estimated.

B. Adopted delayed structure
A delayed in-quadrature complex signal generation technique is adopted.Conceptually, an in-quadrature component, given a sampled single-tone steady-state signal (x(n)) as defined in (1), can be generated by delaying it by f s /(4f 0 ) .Since f 0 is not known in advance, a two-step delayed inquadrature complex signal generation method based on the IpDFT, the so-called Time-Delay QSG (TD-QSG), is proposed.The method is described in Algorithm 1, where the functions DFT and IpDFT refer respectively to (2) and ( 7)- (9).
First, an initial approximation of the in-quadrature complex signal (x o (n)) is generated delaying the sampled signal (x(n)) by d 0 samples.These correspond to the rounded delay given by the nominal frequency f n (lines 1-2).The spectrum of the complex signal is then calculated and windowed in the frequency domain considering a three-cycle Hanning window 4(lines 3-4).A three-point IpDFT is then applied to just obtain a refined estimate of the signal frequency ( f0 ) (line 5), which is subsequently used to obtain a refined delay d f and inquadrature complex signal (x f (n)) (lines 6-7).

IV. TD-IPDFT TECHNIQUE
This section describes the TD-IpDFT algorithm by providing: (i) a description of its structure and formulation (Section IV-A); (ii) a novel trigger mechanism for OOBI identification and removal; (iii) a sensitivity analysis to choose its most suitable parameters (Section IV-B), and (iv) a preliminary analysis of its computational complexity (Section IV-C).The proposed method is an improvement over the i-IpDFT proposed in [23], [30] and its aim is to match the performance of the i-IpDFT while presenting a lower computational cost by using a delayed in-quadrature complex signal to mitigate the effects of self-interference.

A. Structure and Formulation
The TD-IpDFT algorithm incorporating the OOBI compensation routine is summarized using the diagram in Fig. 2 and the pseudocode in Algorithm 2, where the functions TD-QSG, TD-SR, TD-APc and wf refer, respectively, to Algorithm 1, Algorithm 3, Algorithm 4 and ( 5)-( 6), whereas Fig. 3 shows the spectra at different steps of the process.First, the delayed inquadrature complex signal (x f ) is generated using Algorithm 1 (line 1), followed by the calculation of its DFT spectrum and the application of the Hanning window in the frequency domain (lines 2-3).The IpDFT is then applied to obtain a first estimate of the signal parameters (line 4).Subsequently, an interference compensation loop is initiated (line 6) after the initialization of some auxiliary variables (line 5).These include the initial estimates of the positive ( X0 i+ ) and negative ( X0 i-) spectrums of a potential interference tone, the previous energy of the residual spectrum (R 0 e ), along with the residual energy exit (τ Re ) and interference (τ i ) trigger flags.
Within the loop, provided that the residual energy exit flag has not been triggered in the previous iteration, the contribution of the fundamental tone is first estimated ( Xq−1 0 ) by means of Algorithm 3 (line 8).Algorithm 3, named Time-Delay Spectral Reconstruction (TD-SR), calculates the spectrum of any tone 'α' in a delayed in-quadrature complex signal (x f ) given estimates of its frequency ( fα ) and those of the uncorrected amplitude ( Âα+ ) and phase ( φα+ ) of its positive image.This is done by calculating the positive (σ +α ) and negative (σ -α ) complex delay gains (Algorithm 3: lines 1-2), which can then be used to obtain the spectral contributions of both positive and negative images (Algorithm 3: lines 3-6).Together with the previously estimated negative spectrum of the interference tone ( Xq−1 i-), Xq−1 0 is removed from the original spectrum (line 9) and used to determine whether such narrowband components may be present.The analysis of such bins, which should correspond to a positive image of a spurious tone, is performed during the first iteration.Among them, only those where an interference tone is expected are considered, i.e., for a three-cycle window k Once the highest magnitude bin k c has been determined (13), its total energy and that of its two closest neighbors are aggregated into E c (14) (line 11).The calculated E c is then compared with three heuristically defined threshold levels (λ l o , λ u o , λ i ), which allow for the identification of a potential interference (see Section IV-B).These thresholds measure the relation of E c to the total original spectral energy E o (15a) (λ l o , λ u o ) and to the entire set of interference bins For an efficient implementation, a routine stop criterion is introduced inspired by the one proposed in [24], where the incremental ratio between consecutive values of E i /E o was evaluated.Instead, the evolution of the total residual spectral energy R e is used.This is defined as the total remaining energy after subtracting both the estimated fundamental and the interfering tones inferred by the previous iteration, from the original spectrum (16).The adopted metric is more computationally efficient than the one used in [24] since it substitutes divisions with subtractions.If the variation falls below a predefined threshold level (λ Re ), τ Re is set to zero to stop the iterative process (see Section IV-B) (lines [21][22][23][24][25]. The second IpDFT is then applied to estimate the parameters of the positive image of the interfering component (line 22), which are used to approximate the complete spectrum of the interfering tone ( Xq i (k))(line 23).Finally, the last IpDFT can be applied to X f H (k)− Xq i (k), obtaining an improved estimate of the fundamental { f q 0 , Âq 0+ , φq 0+ } (line 24).The process is looped Q times or until τ Re is triggered.The final results, the initial estimate, if no interferences are found, or the latest, once τ Re is triggered or the maximum number of iterations is reached, are then corrected to account for the amplitude and phase alterations due to the use of the delayed inquadrature complex signal (x f ).If no interferences are found, the amplitude and phase corrections can be applied directly (line 26), since the fundmanetal positive complex delay gain (σ q−1 +0 ) is already obtained when Algorithm 3 is applied to estimate the fundamental spectrum (line 8).Otherwise, this is done by means of Algorithm 4 (lines 30 and 35), named Time-Delay Amplitude and Phase correction (TD-APc), which corrects the estimated amplitude and phase of the positive image of any tone 'α' in xf given fα , Âα+ , and φα+ .Finally, fundamental frequency estimations at two successive reporting times (m and m−1) are used to calculate the Rate-Of-Change-Of-Frequency (ROCOF) at the reporting time m with a firstorder backward approximation of a first-order derivative.
F r denotes the reporting rate.

B. Novel Trigger Mechanism for OOBI Compenstation
A novel trigger mechanism for the OOBI that builds on and expands the normalized spectral energy method proposed in [23] is adopted.In [23] an energy threshold capable of identifying interfering components with distortion equal to or greater than 10% was adjusted based on a characterization of the energy content shown by each test condition specified in [4].A first refinement is proposed on the basis of an improved selection of the DFT bins to calculate the spectral energy.Although all bins in the interference spectrum were considered in [23], only k c , and its two closest neighbors, are now selected.The ratio between E c and E o is called spectral energy ratio.As shown in [23], the amplitude and phase steps are the main limiting signal dynamics when setting an adequate energy threshold.However, their energy is spread across the Fig. 3. Spectra at different steps of Algorithm 2 for a discrete signal x(n) affected by additive white Gaussian noise with a signal-to-noise ratio equal to 80dB The rectangular window is used to represent all spectra, so the effects of long-range spectral leakage are noticeable.Similarly, a 35% harmonic distortion is selected for the same reason.whole spectrum, whereas, in the case of an interfering tone, it will be concentrated around its frequency.This is leveraged by introducing an additional metric to enhance the distinction between interfering tones and other spurious contributions of different origin.The so-called spectral energy concentration ratio allows to measure the spread of spectral energy around k c by taking the relation between E c and E i .
Setting appropriate values for λ l o , λ u o and λ i allows for the identification and correction of OOBIs below the limit of 10% set by [4].To account for more realistic conditions and derive more robust thresholds, multi-dynamic signals, beyond the scope of [4], have been simulated.Fig. 4 shows the variability of E c /E o and E c /E i using a boxplot representation for each test condition in [4].All regular tests have been combined with an amplitude modulation (except for the amplitude modulation itself and the OOBI) and evaluated for SNR equal to 60 and 80 dB.Amplitude modulations were selected because they present the second highest energy content, just after the steps.For each test, a 10% amplitude modulation and the highest modulating frequency within the range [0.1−5] Hz, compliant with [4], has been used.Disregarding the amplitude modulated steps (PS* and AS*), Fig. 4 (top) reveals how a minimum of 5% OOBI can be distinguished from the other signal types by selecting a λ l o = 7.4 • 10 −4 .Lower levels of OOBI fall within the phase and amplitude modulated range (PM*).Additionally, the 10% amplitude modulated harmonic distortion (HD 10%*) also exhibits some outliers that exceed this limit.These correspond to the second harmonic which can also be corrected with the iterative process.Moreover, setting a λ u o = 2.4 • 10 −3 allows distinguishing OOBIs with an amplitude greater than or equal to 9% from amplitude modulated steps (PS* and AS*).
Regarding the value of λ i , only an examination of the operating conditions that fall within the band defined by both λ u o and λ l o is necessary.Namely, the amplitude modulated steps.Therefore, a value of 0.765 has been chosen, which allows to account for the spectral energy concentration ratio of all OOBI distortion levels from 5% to 9%.Both OOBI cases are plotted (while the latter represents the most critical case) since the ratio decreases as the distortion level increases.A small overlap can still be seen in Fig. 4 (bottom) regarding the PS* as well as some outliers for the AS*.However, the simulation results according to the standard [4] did not show any impact on their respective tests.All the cases (top and bottom) reveal that noise levels have no significant impact on the magnitude of neither E c /E o nor E c /E i .Only for the latter in the signal frequency (SF*), amplitude modulation (AM*) and harmonic distortion (HD*) tests can the effects of noise be noticed.
Finally, the adjustment of λ Re described in Section IV-A, is analyzed.This is done by simulating 200 consecutive runs of a 10% distortion OOBI test for 23 different potential threshold values.These have been selected within the 10 −8 to 10 −11 range based on a previous coarse estimation.Each run consists of 43 windows of analysis for each interfering tone.The analysis is done accounting for SNRs equal to 60 and 80 dB and allowing the iterative process to run, if unstopped, for a sizeable number of iterations.Fig. 5 shows the variability of the error in estimating the correction term δ and in the total number of iterations executed for a maximum Q equal to 200 (Fig. 5(a)) and 37 (Fig. 5(b)) .The error estimates of δ consider the global maximum error among all simulated interfering tones.Despite a noticeable performance trend in the error estimates of δ for the 80 dB case, a more uniform behavior is shown for a 60 dB noise.As shown in Fig. 5(a) for the 60dB case, although no single run reached the maximum iteration limit (Q) of 200, all considered threshold values λ Re , with the exception of λ Re = 10 −8 and λ Re = 5 • 10 −9 , presented runs where a large number of executed iterations were required.Although the adoption of λ Re = 10 −8 or λ Re = 5•10 −9 would limit computational cost by reducing the maximum number of required iterations, it would also compromise accuracy, as both values result in notably higher error estimates of δ at a lower noise level.Fig. 5(b) examines the impact of limiting the maximum number of iterations (Q) to 37 on the error estimates of δ 5 .The results show that there are no significant changes in the error values obtained compared to Fig. 5(a).Thus, a λ Re = 9.5 • 10 −10 is selected as a trade-off between accuracy and computational performance.This value ensures the lowest maximum error in the estimates of δ (1.3683 • 10 −4 (60 dB)) when considering an absolute maximum of 37 iterations. 5A value of Q = 37 has been selected as it corresponds to the maximum number of iterations required to attain the lowest maximum error in the estimates of δ within the 80dB case in Fig. 5(a).

C. Computational Complexity
To evaluate the feasibility of the implementation of the TD-IpDFT in an embedded device, and compare its performance with other state-of-the-art techniques, its computational complexity is here analyzed.Table I summarizes the total number of arithmetic operations required by the TD-IpDFT (both in case an interfering tone is detected or not) versus those used by its constituent functions, e.g., IpDFT.As in [23] 1 For comparison with other existing methods the operations to obtain the delayed in-quadrature complex spectrum from the buffered windowed DFT bins of the real signal, as shown in Fig. 6, have been included.
and [24], the difference between simple operations (+| − |×), complex operations (÷| exp | sin |∠), and function calls (such as calls to predefined subroutines or algorithms, e.g.IpDFT) is drawn considering a field programmable gate array-based device (FPGA).Likewise, the total cost is expressed in terms of the total number of DFT bins calculated, K, and the maximum number of executions of the iterative process Q.
The results show a total of 675 simple operations and 164 complex operations are required if no interference is detected and 174 + 1065Q simple and 34 + 287Q complex operations considering the OOBI iterative compensation for 8 DFT bins.For a maximum of 37 iterations the total number of operations results in 39579 simple and 10653 complex operations.No explicit method has been specified for the calculation of the DFT bins, nor for generating the delayed in-quadrature complex signal (x f (n)).As indicated in [23], if a small number of DFT bins are needed, recursive computation methods are generally more efficient.Among them, the mSDFT technique [32] is guaranteed stable without sacrificing accuracy and has been used by previous implementations of IpDFT algorithms in embedded FPGA-based PMU devices [29][30].Consequently, Fig. 6 presents a preliminary implementation of the windowed delayed in-quadrature complex signal spectrum.For efficiency, as suggested in [32], a common buffer is used for all mSDFT modules.Moreover, the delayed in-quadrature complex spectrum is obtained by buffering the windowed DFT bins of the real signal X o H (k). As indicated in Algorithm 1 the method relies on an IpDFT to refine the final delay based on an initial nominal quadrature and thus select the appropriate past DFT bins from the buffer D f l .The size of the buffer D f l depends on the lowest expected signal frequency, sampling rate, and the required number of DFT bins K.For a minimum frequency of 45 Hz, 50 kHz sampling rate, and 8 bins, its size will be roughly 1.5 times that of the one required by the mSDFT.

V. PERFORMANCE ASSESSMENT
In this section, a complete evaluation of the TD-IpDFT algorithm is performed by comparing its performance with the static and dynamic accuracy limits indicated in [4] for the performance classes P and M. The assessment is carried out in a MATLAB simulation environment in terms of total vector error (TVE), frequency error (FE), ROCOF error (RFE) and response (R t ) and delay times (D t ) for the step tests.As in [23], [24], a nominal frequency of 50Hz and a reporting rate of 50 frames per second (fps) were selected to limit the number of tests.Similarly, the Hanning window function and a threepoint interpolation scheme are adopted as they, respectively, represent a good trade-off between the width of the main lobe and the levels of the side lobes [21], and reduce the effect of long-range spectral leakage [25].Each test is performed using a three nominal cycle window and considering two levels of additive white Gaussian noise.Noise levels with signal-tonoise ratios (SNR) equal to 60 and 80 dB have been selected, as previously used in [23], [24], as they allow to consider the uncertainty of the measurement and simulate more realistic conditions.The results are presented by means of stacked graphs that summarize the resulting performance against the maximum permissible limits set by [4] (Figs. 7 -14 ).Finally, the maximum values resulting from each test case are reported together again with the accuracy limits indicated in [4] in Tables III, IV, and V.All simulations are carried out according to the parameters given in Table II.

A. Static Tests
Three static tests are defined in [4] to evaluate the performance of the algorithm under steady-state conditions: the signal frequency range test, the harmonic distorsion test and the out-of-band interference (OOBI) test.The results of all static tests are presented in Figs.(7 -9) and Table III.
The results of the signal frequency range test (Fig. 7) show how the accuracy of the method does not depend on the fundamental tone frequency but rather on the total noise level, resulting in errors of one higher order of magnitude for the lower SNR.As reported in Table III, maximum TVE values of 0.009% (60 dB) and 0.001% (80 dB) are obtained, well below the required 1% limit.Likewise, similar trends also apply for 1 Maximum limit values taken from [4] for the 10% case.Fig. 7. Signal frequency range test [4].
the FE and RFE maximum errors, where respective values of 1.15 mHz (60 dB) and 0.13 mHz (80 dB) for the frequency and 0.099 Hz/s (60 dB) and 0.010 Hz/s (80 dB) for the ROCOF are obtained, in accordance with the limit requirements of 5 mHz and 0.1 Hz/s.As was the case in [23], [24], it is worth noting that, under the higher noise conditions, spurious RFE values can marginally exceed the most stringent limit of class M.
The results of the harmonic distortion test are shown in Fig. 8 for a THD of 10% and in Table III for both THDs of 1% and 10%.Again, no significant performance difference is detected on the basis of the order of the interference tone, but rather on the total noise level.The 1% performance limit is far from the maximum TVE values reported for the 60 dB case of 0.024% (1% THD) and 0.022% (10% THD) in Table III.Specifically, the maximum values of FE and RFE were 1.09 mHz (1% THD) and 0.95 mHz (10% THD) for the frequency and 0.094 Hz/s (1% THD) and 0.089 Hz/s (10% THD) for the ROCOF.All are within the most demanding 5 mHz and 0.4 Hz/s limit requirements for devices of class P.
Finally, regarding the OOBI test, Fig. 9 shows the maximum value of TVE, FE and RFE obtained for each interference tone among the three simulated fundamental frequency values of 47.5, 50 and 52.5 Hz considering a total interharmonic distortion of 10% as required by [4].The closer the interference is to the fundamental tone, the more difficult it is to detect and remove it.This trend is evident in Fig. 9 for the lower noise level, as the maximum registered values increase as they approach the fundamental tone frequency.These maximum values, together with those for a 5% distortion, disaggregated by each fundamental frequency, are presented in Table III.All values are well within the performance requirements of class M. It is important to note that no performance requirement is provided in [4] for interfering signals with a magnitude other than 10% of the fundamental.Thus, the same reference values have been considered for the 5% case.As shown, the newly proposed trigger mechanism for OOBI is capable of correcting interferences with amplitudes equal to or greater than 5%, showing a very similar performance to the 10% case.This represents an improvement over the detection mechanism used in [23], [24], which was only able to detect and correct interferences equal to or greater than 10%.

B. Dynamic Tests
Two dynamic tests are defined in [4] to evaluate the performance the algorithm under time-varying conditions: the measurement test and the frequency ramp test.The results of all dynamic tests are presented in Figs. 10 -12 and Table IV.
The most demanding requirements set by [4] correspond to class P, with TVE, FE, and RFE values, respectively, of 3%, 60 mHz, and 2.3 Hz/s for both amplitude and phase modulations.In the case of frequency ramps, class M sets the most stringent limits with TVE, FE, and RFE thresholds of 1%, 10 mHz and 0.2 Hz/s, respectively.As shown in Figs. 10 and 11, both measurement bandwidth tests show how the algorithm exhibits similar TVEs.The total error increases with the modulation frequency and becomes the predominant source of error, masking the effect of noise for modulating frequencies above 1 Hz.However, the trends for the FE and RFE differ.On the one hand, in the phase modulation test (Fig. 11), the TD-IpDFT method reveals a pattern like the one of the TVE where higher modulating frequencies result  in higher errors, and noise as the main source of error is overridden.However, in the amplitude modulation test (Fig. 10), the algorithm can provide accurate frequency and ROCOF estimates without being affected by the modulation frequency.All maximum errors fall well within the limits set in [4].The Fig. 12. Frequency ramp test.Worst-case TVE, FE, and RFE as functions of the ramp rate.[4] worst-case results of the frequency ramp test are shown in Fig. 12 for ramps with different positive and negative rates.The maximum values are presented in Table IV.The TD-IpDFT is shown to provide accurate frequency and ROCOF estimates, in line with those of the signal frequency test and unaffected by the ramp rate, and maximum TVEs of 0.044% (60 dB) and 0.038% (80 dB) which depend on the magnitude of the ramp rate but not on its sign, fully meeting performance requirements.

C. Step Tests
Instantaneous changes in the amplitude (±10%) or phase (±π/18) of the signal are defined in [4] to evaluate the performance of the algorithm in a transient event.The results of both steps are shown in Figs. 13   function of their respective response times, which means that each time axis is centered at the moment when the accuracy limit is first exceeded.Likewise, the estimated phase and amplitude are shown as functions of their respective delay time, with their time axis centered at the instant each step occurs.The results show how the proposed algorithm meets all the requirements, with all estimates within the limits.It is important to note that no significant impact is found by the noise level during the transient, besides the pre-and poststeady-state accuracy already shown by the static tests.

VI. COMPARISON WITH SOTA IPDFT
This section is intended to discuss the results obtained with the TD-IpDFT compared to other IpDFT-based SE methods in the existing literature, namely the i-IpDFT [23] and the HT- IpDFT [24], since all these SE methods have been designed to meet the accuracy requirements for the P and M performance classes.
Focusing on the OOBI, which is the most computationally demanding test, the HT-IpDFT [24] has been shown to offer a lower computational cost compared to the i-IpDFT [23] due to the cancelation of negative frequency components offered by the adoption of the analytic signal.However, the HT-IpDFT does not meet the combined requirements of classes P and M for harmonic distortion (1%) and phase step tests in [4].The TD-IpDFT satisfies all the accuracy requirements for the P and M classes, and also offers a reduction in the total computational cost compared to the i-IpDFT [23], [30].This is because it mitigates the effects of self-interference by using the in-quadrature complex signal for estimate, while the i-IpDFT, which is based on the e-IpDFT [22], must iteratively for them on each successive iteration 6 .Taking into account the same parameterization for both methods, that is, the same type and length of the window and the same maximum number of iterations and DFT bins, the i-IpDFT requires a total of 844 + 1751Q simple operations and 218 + 474Q complex operations based on the updated tuning presented in [30], whereas the TD-IpDFT can perform the estimation with just 174 + 1065Q and 34 + 287Q operations respectively.This represents approximately, for simple and 6 To draw a meaningful comparison a simplified version of Algorithm 2, aimed at single-tone signals, can be defined consisting of lines 1-4 and line 30.This simplified TD-IpDFT can be shown to be more computationally efficient than the e-IpDFT.Taking into account the 3p variant of the e-IpDFT employed in the i-IpDFT and the same parameterization for both methods, the e-IpDFT requires 29K + 38 simple and 7K + 22 complex operations, while the simplified TD-IpDFT only requires 12K + 38 and 17 operations respectively.complex operations, a total decrease of 40% compared to the i-IpDFT.Moreover, considering the maximum overall error among the entire OOBI range, the TD-IpDFT also shows faster convergence, requiring 9 to 14 iterations less, for the same maximum error level compared to the i-IpDFT.This is shown in Fig. 15 where a performance comparison between TD-IpDFT and i-IpDFT is presented for noise levels with SNR equal to 60 and 80 dB.Both the error in estimating the correction term δ as a function of the iteration number (top) and the FE obtained for the OOBI test considering Q = 37 (bottom) are shown.It can be seen that the TD-IpDFT requires fewer iterations to achieve the same maximum error level compared to the i-IpDFT indistinctly from the noise level.Furthermore, for Q = 37, the TD-IpDFT shows a modest performance improvement throughout the OOBI range for the high-noise case.Finally, the TD-IpDFT, as opposed to the i-IpDFT, also delivers a complex signal, which allows one to determine its envelope and angle.These could be used to implement identification and correction techniques for the amplitude and phase steps, as done in [33].

VII. CONCLUSIONS
In this article, a SE technique based on the application of the IpDFT to a delayed in-quadrature complex signal has been presented.The use of the delayed in-quadrature complex signal has been found to be a simple and efficient way to mitigate the detrimental effects caused by the self-interference of the fundamental tone.This allows for a reduction in the total computational complexity of approximately 40% compared to the i-IpDFT, which relies on an iterative correction based on the e-IpDFT.This is especially significant when the OOBI routine is triggered, since a faster convergence, requiring up to 14 fewer iterations for the same maximum error level, is obtained for the higher noise case.However, an increase of around 150% in memory requirements is necessary to generate the delayed signal.The method represents an alternative to the i-IpDFT and the HT-IpDFT.Like the i-IpDFT, it satisfies all the accuracy requirements defined in the IEC/IEEE Std.60255-118, while presenting a lower computational cost and it represents an easier and faster implementation than the cascaded Hilbert filter used by the HT-IpDFT.Furthermore, the novel OOBI trigger improves on the one used by the i-IpDFT and by the HT-IpDFT, which only allowed correction of distortions equal to or greater than 10%.Future steps will include the implementation of the TD-IpDFT in an embedded device and its experimental validation.

Fig. 1 .
Fig.1.Frequency response of a 'filtered' delayed in-quadrature complex signal ȳ(n) generated considering a frequency f from a signal y(n) with normalized frequency f 0 [pu] = f 0 /f .Wide view (top) and zoomed (bottom).

2 Fig. 2 .
Fig. 2. Diagram of the TD-IpDFT Algorithm.The process is looped Q + 1 times instead of Q as in the pseudocode in Algorithm 2 since the IpDFT function on line 4 is nested within the loop.

Fig. 5 .
Fig. 5. Boxplot representation of the error in estimating the correction term δ (top) and the total number of iterations (bottom) for different values of λ Re .The maximum number of iterations Q is set to 200 (a) and 37 (b).

Fig. 15 .
Fig. 15.Performance comparison between TD-IpDFT and i-IpDFT for noise levels with SNR equal to 60 and 80 dB.Error in estimating the correction term δ as a function of the iteration number (top) and FE as a function of the interference tone frequency for Q = 37 (bottom).

TABLE III MAXIMUM
[4], FE AND RFE IN STATIC TESTS AND MAXIMUM LIMIT ALLOWED BY[4]

TABLE IV MAXIMUM
[4], FE AND RFE IN DYNAMIC TESTS AND MAXIMUM LIMIT ALLOWED BY[4] and 14, andtheir response, delay times, and maximum overshoot values are presented in Table V.All results correspond to the positive step cases (similar results are obtained in the case of a negative step).Both Figs. 13 and 14 represent TVE, FE, and RFE as a