Harmonic Phasor Estimation Method Considering Dense Interharmonic Interference

Due to the limitation of frequency resolution and the spectrum leakage caused by signal windowing, the spectrums of harmonic and interharmonic components with close frequencies overlap each other. When the dense interharmonic (DI) components are close to the harmonic spectrum peaks, the harmonic phasor estimation accuracy is seriously reduced. To address this problem, a harmonic phasor estimation method considering DI interference is proposed in this paper. Firstly, based on the spectral characteristics of the dense frequency signal, the phase and amplitude characteristics are used to determine whether DI interference exists in the signal. Secondly, an autoregressive model is established by using the autocorrelation of the signal. Data extrapolation is performed on the basis of the sampling sequence to improve the frequency resolution and eliminate the interharmonic interference. Finally, the estimated values of harmonic phasor, frequency and rate of change of frequency are obtained. The simulation and some experimental results demonstrate that the proposed method can accurately estimate the parameters of harmonic phasors when DIs exist in the signal, and has a certain anti-noise capability and dynamic performance.


Introduction
With a large amount of renewable energy and nonlinear loads connected to the grid, harmonics and interharmonics are becoming increasingly abundant, posing a great threat to the safe and stable operation of the grid [1,2]. On the one hand, some large industrial equipment can generate dense interharmonics (DIs) with close frequencies [3]. At the same time, many harmonic components are considered as interharmonics due to the asynchronous sampling [4]. When these signals are mixed together, serious spectral interference is generated, leading to distortion of the measurement of harmonic parameters. Therefore, it is significant to study the estimation method for harmonic phasor under DI interference.
At present, the estimation methods of harmonic phasor can be divided into two categories: non-signal-transformation-based and signal-transformation-based methods. The non-signal-transformation-based methods include model-based and mode-decompositionbased methods. Among them, the model-based methods, such as the Prony algorithm [5], least squares method and rotational invariance technique (ESPRIT) [6], are employed. On the one hand, these methods have a deficiency of adaptive capability and cannot follow the dynamic signal to change the model coefficients in time. On the other hand, these methods will cause large errors when the signal to noise ratio (SNR) is low. The mode-decomposition-based methods include Hilbert transform (HT) and variational mode decomposition (VMD) [7]. However, these methods require a pre-set number of modes to be decomposed. If the number is not set properly, the decomposition results will be biased.
The signal-transformation-based methods can be further divided into static signal estimation and dynamic signal estimation methods. The static signal estimation methods include discrete Fourier transform (DFT), S-transform, wavelet transform, etc. Among them, DFT is simple to implement, has better stability and suitability, and has gained wide neural networks [28], although they have a high accuracy, all require some prior information and their applications are restricted.
To address the above problems, this paper proposes a harmonic phasor estimation method considering DI interference. Firstly, the estimation model for harmonic phasor is established. Based on the characteristics of the spectral lines with spectral interference, we use the phase and amplitude characteristics of the spectral lines to determine whether there are DIs in the signal. Secondly, for the signal containing DI interference, an autoregressive model is established using the sampled values based on the autocorrelation of harmonic data. Moreover, data extrapolation is performed to improve the frequency resolution, so that all interharmonic components can be estimated and eliminated. Finally, the estimates of harmonic phasor, frequency and ROCOF are calculated. Simulation and experimental results demonstrate that under the conditions of asynchronous sampling and DI interference, the proposed method can accurately estimate harmonic phasors. In addition, it has a certain dynamic performance and anti-noise capability, and presents superiority in comparison with other methods. The proposed method belongs to the dynamic signal estimation methods. Since the Taylor-Fourier transform is used, the proposed method can also be classified as signal-transformation-based methods.
This paper is structured as follows. Section 2 presents the estimation model for harmonic phasors and the interference caused by DIs. Section 3 describes the proposed method in detail. Section 4 shows the simulation and experimental results. Finally, Section 5 concludes the whole paper.

Harmonic Phasor Estimation Model
In general, the actual signal containing harmonics can be represented in the time domain as where a h (t) and ϕ h (t) are the amplitude and phase of the h-th harmonic, respectively. f h is the harmonic frequency. H is the maximum number of harmonics. According to the IEEE standard [29], the harmonic phasor P h (t) corresponding to f h is expressed as follows: The Taylor-Fourier (TF) model can describe the changes of amplitude and phase in the observation time. Therefore, a Taylor expansion model of order K h is used to represent P h (t). By substituting the phasor P h (t) after the Taylor expansion into the expression of s(t), we can obtain where p h (k) is the TF coefficient, which is the k-th order derivative of P h (t) at time t = 0. * is the conjugate operator. For different h-th harmonics, the order of Taylor expansion may be different. Equation (3) establishes the TFM model, which provides the basis for various analyses, including the estimation of harmonic phasors.

Interference of DIs on Harmonic Phasor Estimation
DIs can be understood as interharmonics with close frequencies. Due to the signal windowing in the time domain, DFT produces spectral leakage, leading to the appearance of the main and side lobes of each frequency component. When two adjacent spectral peaks of interharmonics are close to each other, the main and side lobes of spectral peak A are overlapped onto the main and side lobes of spectral peak B. This will cause serious spectral interference problems and thus lead to deviations in the positions of the spectral peaks of the dense components in the spectrum. Especially when the interharmonic spectral peaks are close to the harmonic spectral peaks, the estimation accuracy of the harmonic phasor will be seriously affected. Figure 1 shows the spectrum diagram of the signal x(t) = 0.5cos(2π × 95.5t) + cos(2π × 100t) + 0.5cos(2π × 102.8t). The sampling frequency is 5 kHz and the number of sampling points is 1024. From Figure 1, two DIs exist near the 2nd harmonic, causing serious spectral interferences. Only two spectral peaks can be observed, and their corresponding frequencies and amplitudes deviate significantly from the true values. In order for the calculated spectrum to reveal the adjacent spectral peaks, it is necessary that the frequency difference between the two adjacent components is greater than or equal to the frequency resolution. Although this could be achieved by increasing the number of sampling points, that is, increasing the observation interval, this would reduce the dynamic tracking capability for the signal. The proposed method can improve the frequency resolution while keeping the number of sampling points unchanged, thus allowing the estimation of harmonic phasors under DI interference.

Interference of DIs on Harmonic Phasor Estimation
DIs can be understood as interharmonics with close frequencies. Due to the signal windowing in the time domain, DFT produces spectral leakage, leading to the appearance of the main and side lobes of each frequency component. When two adjacent spectral peaks of interharmonics are close to each other, the main and side lobes of spectral peak A are overlapped onto the main and side lobes of spectral peak B. This will cause serious spectral interference problems and thus lead to deviations in the positions of the spectral peaks of the dense components in the spectrum. Especially when the interharmonic spectral peaks are close to the harmonic spectral peaks, the estimation accuracy of the harmonic phasor will be seriously affected. Figure 1 shows the spectrum diagram of the signal x(t) = 0.5cos(2π × 95.5t) + cos(2π × 100t) + 0.5cos(2π × 102.8t). The sampling frequency is 5 kHz and the number of sampling points is 1024. From Figure 1, two DIs exist near the 2nd harmonic, causing serious spectral interferences. Only two spectral peaks can be observed, and their corresponding frequencies and amplitudes deviate significantly from the true values. In order for the calculated spectrum to reveal the adjacent spectral peaks, it is necessary that the frequency difference between the two adjacent components is greater than or equal to the frequency resolution. Although this could be achieved by increasing the number of sampling points, that is, increasing the observation interval, this would reduce the dynamic tracking capability for the signal. The proposed method can improve the frequency resolution while keeping the number of sampling points unchanged, thus allowing the estimation of harmonic phasors under DI interference.

Determination of Spectral Interference
The presence of interference between spectral lines can be determined by the phase and amplitude characteristics of the spectral lines [30]. Under the condition that the spectral lines are accurate, the phase difference between adjacent spectral lines is 180°. If the phase difference is not 180°, it means that there is a spectral interference. In practical application, two spectral lines A and B near the spectral peak can be selected, and the phase difference Δθ is calculated.
where ε1 is the set error threshold. If (4) is satisfied, it is necessary to further determine whether there is a spectral interference by the amplitude characteristic. If it is not satisfied, then there is a spectral interference. For the determination of the amplitude characteristic, three adjacent spectral lines, A, B and C, near the spectral peak are selected. The corrections are done by A and B as well as B and C, respectively. If the results are the same, there is no spectral interference. Firstly, the ratio of adjacent spectral lines is calculated.

Determination of Spectral Interference
The presence of interference between spectral lines can be determined by the phase and amplitude characteristics of the spectral lines [30]. Under the condition that the spectral lines are accurate, the phase difference between adjacent spectral lines is 180 • . If the phase difference is not 180 • , it means that there is a spectral interference. In practical application, two spectral lines A and B near the spectral peak can be selected, and the phase difference ∆θ is calculated.
where ε 1 is the set error threshold. If (4) is satisfied, it is necessary to further determine whether there is a spectral interference by the amplitude characteristic. If it is not satisfied, then there is a spectral interference. For the determination of the amplitude characteristic, three adjacent spectral lines, A, B and C, near the spectral peak are selected. The corrections are done by A and B as well as B and C, respectively. If the results are the same, there is no spectral interference. Firstly, the ratio of adjacent spectral lines is calculated.
where y A , y B and y C are the corresponding DFT amplitudes of spectral lines A, B and C, respectively. In this paper, the rectangular window is used as an example. The frequency correction ∆k of the rectangular window can be expressed as ∆k = (v − 1)/(v + 1), where ∆k is the frequency correction of the spectral line and v is the ratio of the amplitudes of two adjacent spectral lines. The frequency corrections ∆f A and ∆f B of the spectral line A and B can be obtained by bringing v 1 and v 2 into the expression of ∆k, respectively. The difference of the two corrections is used as a criterion. Since the difference between the numbers of the spectral lines A and B is 1, the correction difference ∆P is calculated.
where ε 2 is the set error threshold. If ε 2 is 0.01, it means that the error of ∆P is 0.01 times of the frequency resolution. If (6) is satisfied, there is no spectral interference.

Spectrum Refinement Algorithm
For the spectrum containing DI interference, this paper proposes a spectrum refinement algorithm to improve the frequency resolution while keeping the total sampling time unchanged. Therefore, the spectral peaks of DIs can be distinguished and their parameters can be estimated. The basic idea of the proposed method is to build an autoregressive (AR) model based on the existing sampling sequence. Data extrapolation is performed based on the established AR model, thus the sequence space is expanded. Afterwards, DFT analysis is performed on the expanded sequence space. Therefore, the frequency resolution is improved because the sampling frequency remains unchanged. When the frequency resolution is less than the frequency difference of adjacent interharmonics, the interharmonic parameters can be accurately estimated. The harmonic phasors can then be further estimated.
The AR model for the sequence x(n) can be expressed by the following equation [31]: where p is the order of the AR model. ε n is a random error term with a mean of 0 and a standard deviation of σ. The premise of the AR model is that the sequence must have an autocorrelation property. Considering that harmonic data in power systems have certain autocorrelations [32], the signal can be modeled using (7). First, the order p is calculated using the singular value decomposition method. The autocorrelation function r x (k) of the sequence x(n) is calculated as follows: where N is the number of the sampled sequence x(n). The autocorrelation coefficient of x(n) can be expressed as r x (1). The autocorrelation matrix R of x(n) is further calculated.
The matrix R is diagonalized using the singular value decomposition as follows: where U and V are unitary matrices. λ 1 , λ 2 , . . . , λ N are the eigenvalues of matrix R and satisfy λ 1 ≥ λ 2 ≥ . . . ≥ λ N ≥ 0. Then, the order of the AR model can be obtained by the following equation: where T h is the set threshold value. Next, the AR model coefficients a(i) and the random error σ 2 are calculated. The forward prediction error f p (n), and the backward prediction error g p (n) are defined as follows: The energy of the error should be minimized when the forward and backward prediction is optimal. In addition, it is necessary to satisfy the orthogonality of the forward and backward prediction errors with the signal. Based on this, the objective functions P 1 (k), P 2 (k) and P 3 (k) are defined as follows: where k is an iterative variable, k = 1, 2, . . . , p. The expressions for f k (n) and g k (n) are obtained by replacing p with k in (12). The objective function J(k) is constructed as follows: The estimates of a(i) and σ 2 are obtained by the minimization of J(k). According to the Levinson recursive algorithm, the reflection coefficients β 1 , β 2 , . . . , β p are introduced. The recursive formulas of a(i) and σ 2 can be expressed as follows: where a k (i) represents the value of the k-th iteration of a(i). To obtain the estimates of the reflection coefficients, the following equation can be derived using the recurrence of the forward and backward prediction errors. Substituting (16) into (14), and let the partial derivative of J(k) with respect to β k to be 0. Thus, the estimate of β k is obtained. The initial values for the iterations in (14) are In practical application, let k = 1, after which β k , a k (i), σ k 2 , f k (n) and g k (n) are computed in turn. Then, let k = k + 1, and repeat the above steps until a p (i) and σ p 2 are calculated, which are the AR model coefficients a(i) and the random error σ 2 . After obtaining the AR model for the sequence x(n), data extrapolation can be performed. The additional number of sequence elements to be predicted, ∆N, can be expressed as where f s is the sampling frequency. ∆f and ∆f denote the frequency resolution before and after refinement, respectively. The expanded sequence is DFT analyzed and the interharmonic parameters are calculated, thus eliminating the interference of DIs on the harmonic phasor estimation.

Harmonic Phasor Estimation
After eliminating the interharmonic interference, the signal s(t) is sampled at a sampling interval T s = 1/f s . Then, the discrete representation of the signal s(t) is [33] where ω h = 2πf h is the actual angular frequency in radian units. Equation (19) can be further written in a matrix form: is the H(K h + 1) × 1 order column vector of TF coefficients. The estimation of the matrix P can be obtained by the least squares method: where H is the Hermitian operator. After obtaining the estimate of P, the estimate of the harmonic phasor P h (t) within the observation window can be obtained by Taylor formula. Based on the obtained TF coefficients, the derivative estimates of P h (t) can also be obtained as follows:P h is the estimated value of the TF coefficient. T w is the length of the observation window. In this where d 0h , d 1h and d 2h are the zero-order, first-order and second-order derivative estimates of P h (t), respectively. Re and Im are the operations to extract the real and imaginary parts of the complex numbers, respectively. In summary, the proposed harmonic phasor estimation method considering the DI interference is shown in Figure 2.
Entropy 2023, 25, 236 8 of 18 where d0h, d1h and d2h are the zero-order, first-order and second-order derivative estimates of Ph(t), respectively. Re and Im are the operations to extract the real and imaginary parts of the complex numbers, respectively. In summary, the proposed harmonic phasor estimation method considering the DI interference is shown in Figure 2.

Truncate signal s(t) with a Hanning window, and perform
DFT on the intercepted signal to obtain the spectrum

Start
Calculate the order of AR model according to (8)− (11) Initialization: calculate the initial value of the iteration according to (17), let k = 1 Calculate a k (i) and σ k 2 according to (15) k = p ?
Calculate f k (n) and g k (n) according to (16) Obtain the AR model coefficients a(i) and the random error σ 2 Perform data extrapolation according to (18) to estimate the interharmonic parameters Construct matrices S, Φ, P according to (20) Obtain an estimate of matrix P by least squares method Calculate the estimates of harmonic phasor P h (t), harmonic frequency f h and ROCOF h according to (22)

Simulation Analysis
In this section, the performance of harmonic phasor estimation of the proposed method and other methods are compared and analyzed under different simulation conditions. The comparison methods are as follows: Method 1 is IpDFT, Method 2 is the Prony algorithm, Method 3 is CS-TFM, and Method 4 is SIFE. The proposed method is considered as Method 5. According to the evaluation criteria for PMU in IEEE Std C37.118.1-2011, this paper uses the total vector error (TVE), frequency error (FE), and ROCOF error (RFE) to measure the effectiveness of harmonic phasor estimation. To simulate the worst scenarios, this paper refers to the measurement standards of both P and M class PMUs. The simulation scenarios include: steady-state test, dense interharmonics test, wideband noise test, frequency deviation test, dynamic-state test and experimental signal test.

Steady-State Test
The test signal can be expressed by the following equation: where f0 is set to 50 Hz, φ1 and φh are set to random numbers uniformly distributed between (0, 2π). The amplitude of each harmonic is set to 10% of the fundamental

Simulation Analysis
In this section, the performance of harmonic phasor estimation of the proposed method and other methods are compared and analyzed under different simulation conditions. The comparison methods are as follows: Method 1 is IpDFT, Method 2 is the Prony algorithm, Method 3 is CS-TFM, and Method 4 is SIFE. The proposed method is considered as Method 5. According to the evaluation criteria for PMU in IEEE Std C37.118.1-2011, this paper uses the total vector error (TVE), frequency error (FE), and ROCOF error (RFE) to measure the effectiveness of harmonic phasor estimation. To simulate the worst scenarios, this paper refers to the measurement standards of both P and M class PMUs. The simulation scenarios include: steady-state test, dense interharmonics test, wideband noise test, frequency deviation test, dynamic-state test and experimental signal test.

Steady-State Test
The test signal can be expressed by the following equation: where f 0 is set to 50 Hz, ϕ 1 and ϕ h are set to random numbers uniformly distributed between (0, 2π). The amplitude of each harmonic is set to 10% of the fundamental amplitude. According to [23,34], harmonics after the 13th order in the utility grid have a low proportion of the fundamental. Therefore, the maximum number of harmonics in the simulation is 13. The sampling frequency is 5 kHz and the length of the window function is 6 fundamental periods. In the subsequent simulation analysis, the length of the window function will be kept constant at 0.12 s. Because it is suitable for the dynamic characteristics of varying network signals. The signals shown in (25) are detected using Methods 1 to 5, respectively. Figure 3 shows the maximum TVEs, FEs, and RFEs for each method, where the harmonic order of 1 indicates the fundamental wave. amplitude. According to [34,23], harmonics after the 13th order in the utility grid have a low proportion of the fundamental. Therefore, the maximum number of harmonics in the simulation is 13. The sampling frequency is 5 kHz and the length of the window function is 6 fundamental periods. In the subsequent simulation analysis, the length of the window function will be kept constant at 0.12 s. Because it is suitable for the dynamic characteristics of varying network signals. The signals shown in (25) are detected using Methods 1 to 5, respectively. Figure 3 shows the maximum TVEs, FEs, and RFEs for each method, where the harmonic order of 1 indicates the fundamental wave.  As shown in Figure 3, all the methods can provide estimates with errors close to zero at nominal frequencies. The maximum TVE and FE limits in [34] for the synchrophasor measurement are 1% and 0.005 Hz, respectively, when harmonic interference is included. The errors of the proposed method satisfy the requirements completely. Since the signal model contains multiple harmonics, Method 1 has problems of the spectrum leakage, fence effect and mutual harmonic interference. Method 2 requires a large number of exponential functions to fit the signal, and it is difficult to design the corresponding filter. Method 3 uses the second-order Taylor formula to model the phasors, so there is a certain amount of model error. Increasing the order of the Taylor model will improve the accuracy, but the higher order increases the computational burden. Method 4 is difficult to determine the coefficients of the sinc interpolation functions at a high harmonic order. Only Method 5 gives the most accurate estimation results.

Dense Interharmonics Test
The test signal can be expressed by the following equation: The sampling frequency fs is 5 kHz. The number of sampling points is 600. The specific signal parameters are shown in Table 1. From Table 1, the signal contains fundamental, second and third harmonics. Interharmonic interference exists near the fundamental and second harmonic. Figure 4 is the spectrum of s(t) obtained after the DFT.  As shown in Figure 3, all the methods can provide estimates with errors close to zero at nominal frequencies. The maximum TVE and FE limits in [34] for the synchrophasor measurement are 1% and 0.005 Hz, respectively, when harmonic interference is included. The errors of the proposed method satisfy the requirements completely. Since the signal model contains multiple harmonics, Method 1 has problems of the spectrum leakage, fence effect and mutual harmonic interference. Method 2 requires a large number of exponential functions to fit the signal, and it is difficult to design the corresponding filter. Method 3 uses the second-order Taylor formula to model the phasors, so there is a certain amount of model error. Increasing the order of the Taylor model will improve the accuracy, but the higher order increases the computational burden. Method 4 is difficult to determine the coefficients of the sinc interpolation functions at a high harmonic order. Only Method 5 gives the most accurate estimation results.

Dense Interharmonics Test
The test signal can be expressed by the following equation: The sampling frequency f s is 5 kHz. The number of sampling points is 600. The specific signal parameters are shown in Table 1. From Table 1, the signal contains fundamental, second and third harmonics. Interharmonic interference exists near the fundamental and second harmonic. Figure 4 is the spectrum of s(t) obtained after the DFT.  In Figure 4, only three spectral peaks can be observed, corresponding to th cies of 51.2 Hz, 102.4 Hz, and 153.6 Hz. Not only the spectral peaks of the inter cannot be detected, but also the frequency estimations of the harmonics are b autocorrelation coefficient of the above signal was calculated to be 0.997. The AR model can be used for sequence extrapolation. In Figure 4, the frequency Δf = fs/N ≈ 8.33 Hz. In this paper, the refined frequency resolution Δf′ was set The proposed spectrum refinement algorithm is applied to s(t), and the refined is obtained as shown in Figure 5. From Figure 5, the refined spectrum can completely show the spectral pe sponding to each frequency component. In addition, the amplitudes of the spe are consistent with the real values, which proves the effectiveness of the method. The final results of the estimations of signal parameters for each m shown in Table 2.  In Figure 4, only three spectral peaks can be observed, corresponding to the frequencies of 51.2 Hz, 102.4 Hz, and 153.6 Hz. Not only the spectral peaks of the interharmonics cannot be detected, but also the frequency estimations of the harmonics are biased. The autocorrelation coefficient of the above signal was calculated to be 0.997. Therefore, the AR model can be used for sequence extrapolation. In Figure 4, the frequency resolution ∆f = f s /N ≈ 8.33 Hz. In this paper, the refined frequency resolution ∆f was set to 0.1 Hz. The proposed spectrum refinement algorithm is applied to s(t), and the refined spectrum is obtained as shown in Figure 5.  In Figure 4, only three spectral peaks can be observed, corresponding to th cies of 51.2 Hz, 102.4 Hz, and 153.6 Hz. Not only the spectral peaks of the inter cannot be detected, but also the frequency estimations of the harmonics are b autocorrelation coefficient of the above signal was calculated to be 0.997. The AR model can be used for sequence extrapolation. In Figure 4, the frequency Δf = fs/N ≈ 8.33 Hz. In this paper, the refined frequency resolution Δf′ was set The proposed spectrum refinement algorithm is applied to s(t), and the refined is obtained as shown in Figure 5. From Figure 5, the refined spectrum can completely show the spectral p sponding to each frequency component. In addition, the amplitudes of the spe are consistent with the real values, which proves the effectiveness of the method. The final results of the estimations of signal parameters for each m shown in Table 2.  From Figure 5, the refined spectrum can completely show the spectral peaks corresponding to each frequency component. In addition, the amplitudes of the spectral peaks are consistent with the real values, which proves the effectiveness of the proposed method. The final results of the estimations of signal parameters for each method are shown in Table 2.
From Table 2, the interharmonic parameters cannot be obtained by Method 1, and there is a large error in the estimation of the harmonic parameters. This is due to the small amplitude of the interharmonic components, as their spectral peaks are covered by the main and side lobes of harmonics. They also interfere with the spectral peaks of the harmonic phasors. The estimations of harmonics are relatively accurate for Methods 2 to 4, but the errors of interharmonics are larger. For Method 2, the presence of DIs increases the complexity of the model. The signal requires more exponential functions for fitting, and the error increases accordingly. This also has a negative impact on Method 3, which increases the computational complexity. The estimation results of Method 3 and Method 4 are closer to the true values, but there are still some deviations. The TF coefficients, filter coefficients and length of the observation window of Method 3 are difficult to determine. Method 4 only considers the first-order and zero-order derivatives of harmonic phasors, which also has certain model errors. Method 5 uses the autocorrelation of the signal for data extrapolation, which improves the frequency resolution and eliminates the spectral interference. Therefore, only Method 5 gives the most accurate estimation results.

Wideband Noise Test
In practice, the sampled signal often contains a certain amount of wideband noise. The simulation signal used in this part is expressed as follows: 0.01 cos(2π f i t + ϕ i ) + noise (27) where noise is Gaussian white noise with a mean value of 0. f i is the interharmonic frequency and the values are 48.9 Hz, 51.1 Hz, 52.15 Hz, 65.5 Hz, 103 Hz, 213.4 Hz. ϕ 1 , ϕ h and ϕ i are set as random numbers uniformly distributed between (0, 2π). Other parameters are consistent with Section 4.2. Two noise conditions are considered: 40 dB SNR and 60 dB SNR, respectively. Synchrophasor and harmonic phasor detection are performed for signal (27). The maximum TVEs, FEs and RFEs for each method are shown in Figure 6.  Table 2, the interharmonic parameters cannot be obtained by Method 1, and there is a large error in the estimation of the harmonic parameters. This is due to the small amplitude of the interharmonic components, as their spectral peaks are covered by the main and side lobes of harmonics. They also interfere with the spectral peaks of the harmonic phasors. The estimations of harmonics are relatively accurate for Methods 2 to 4, but the errors of interharmonics are larger. For Method 2, the presence of DIs increases the complexity of the model. The signal requires more exponential functions for fitting, and the error increases accordingly. This also has a negative impact on Method 3, which increases the computational complexity. The estimation results of Method 3 and Method 4 are closer to the true values, but there are still some deviations. The TF coefficients, filter coefficients and length of the observation window of Method 3 are difficult to determine. Method 4 only considers the first-order and zero-order derivatives of harmonic phasors, which also has certain model errors. Method 5 uses the autocorrelation of the signal for data extrapolation, which improves the frequency resolution and eliminates the spectral interference. Therefore, only Method 5 gives the most accurate estimation results.

Wideband Noise Test
In practice, the sampled signal often contains a certain amount of wideband noise. The simulation signal used in this part is expressed as follows:  From Figure 6, the errors of Method 2 increase rapidly, which is because Method 2 is extremely sensitive to noise. For wideband noise, it is impossible to use a suitable From Figure 6, the errors of Method 2 increase rapidly, which is because Method 2 is extremely sensitive to noise. For wideband noise, it is impossible to use a suitable exponential function model to fit the signal due to its irregularity. Therefore, accurate estimates of the synchrophasor and harmonic phasors cannot be obtained. The errors of Methods 3 and 4 also increase compared to Sections 4.1 and 4.2. This is because Methods 3 and 4 tend to use an observation window of shorter length in order to obtain a shorter algorithm response time. Therefore, the presence of noise can seriously interfere with the accuracy of the fitting of the Taylor formula and the sinc interpolation formula. Although Method 1 uses windowing and interpolation, there is no procedures to deal with the noise, and thus the error is increased. In contrast, Method 5 is able to accurately identify harmonic phasors under the condition of wideband noise and has an advantage of antinoise capability.

Frequency Deviation Test
To evaluate the performance of the proposed method under the condition of frequency deviation, the following signal was used for testing: According to the latest IEEE standard [34], f was set to vary from 45 Hz to 55 Hz in a step of 1 Hz. Other parameters were consistent with Section 4.3. For Method 1, the spectrum is scattered throughout the frequency domain for signals with varying frequencies, the spectral peaks cannot be accurately identified. Since Method 1 has high errors in the previous tests, it is not discussed in the subsequent analysis. Figure 7 shows the errors in the estimation of harmonic phasors by the proposed method at different fundamental frequencies. Figure 8 shows the estimation errors of the 1st to 7th harmonic phasors for Methods 2 to 5. The maximum TVEs, FEs and RFEs in Figure 8 are the average values of the errors obtained at different fundamental frequencies.
Entropy 2023, 25,236 12 of 18 exponential function model to fit the signal due to its irregularity. Therefore, accurate estimates of the synchrophasor and harmonic phasors cannot be obtained. The errors of Methods 3 and 4 also increase compared to Sections 4.1 and 4.2. This is because Methods 3 and 4 tend to use an observation window of shorter length in order to obtain a shorter algorithm response time. Therefore, the presence of noise can seriously interfere with the accuracy of the fitting of the Taylor formula and the sinc interpolation formula. Although Method 1 uses windowing and interpolation, there is no procedures to deal with the noise, and thus the error is increased. In contrast, Method 5 is able to accurately identify harmonic phasors under the condition of wideband noise and has an advantage of anti-noise capability.

Frequency Deviation Test
To evaluate the performance of the proposed method under the condition of frequency deviation, the following signal was used for testing: According to the latest IEEE standard [34], f was set to vary from 45 Hz to 55 Hz in a step of 1 Hz. Other parameters were consistent with Section 4.3. For Method 1, the spectrum is scattered throughout the frequency domain for signals with varying frequencies, the spectral peaks cannot be accurately identified. Since Method 1 has high errors in the previous tests, it is not discussed in the subsequent analysis. Figure 7 shows the errors in the estimation of harmonic phasors by the proposed method at different fundamental frequencies. Figure 8 shows the estimation errors of the 1st to 7th harmonic phasors for Methods 2 to 5. The maximum TVEs, FEs and RFEs in Figure 8 are the average values of the errors obtained at different fundamental frequencies.  In [34], the TVE and FE limits of M-class PMU for synchrophasor measurement under the condition of frequency deviation are 1% and 0.025 Hz, respectively. From Figure 7, even if this standard is applied to harmonic phasors, Method 5 can still satisfy this standard in the case of small frequency deviations. This is because Method 5 uses the first and second order derivatives to approximate the dynamic phasors. It has a good dynamic performance and still has a high tracking ability under the frequency deviation. Method 2 still does not overcome the drawbacks described in Sections 4.1-4.3, and the exponential function model cannot be adjusted to follow the dynamic signal. The observation window used in Method 3 may not contain sufficient harmonic information. Moreover, the varying harmonic and interharmonic components can enhance spectral interference, which In [34], the TVE and FE limits of M-class PMU for synchrophasor measurement under the condition of frequency deviation are 1% and 0.025 Hz, respectively. From Figure 7, even if this standard is applied to harmonic phasors, Method 5 can still satisfy this standard in the case of small frequency deviations. This is because Method 5 uses the first and second order derivatives to approximate the dynamic phasors. It has a good dynamic performance and still has a high tracking ability under the frequency deviation. Method 2 still does not overcome the drawbacks described in Sections 4.1-4.3, and the exponential function model cannot be adjusted to follow the dynamic signal. The observation window used in Method 3 may not contain sufficient harmonic information. Moreover, the varying harmonic and interharmonic components can enhance spectral interference, which generates large estimation errors. According to Figure 8, the estimation errors of Method 4 and Method 5 are relatively close. Further simulation conditions need to be established.

Dynamic-State Test
In this section, frequency ramp and harmonic oscillation tests will be simulated. First, the frequency ramp test is performed. Due to disturbances in the power system, the frequencies of the voltage or current signals may vary linearly. To further compare the SIFE (Method 4) and the proposed method (Method 5), the following signal was used for testing: where R1 is the slope of the fundamental frequency and is set to 1 Hz/s. f varies linearly from 45 Hz to 55 Hz. The other parameters are identical to those in Section 4.3. The test results are shown in Table 3, where h = 1 represents the fundamental wave.

Dynamic-State Test
In this section, frequency ramp and harmonic oscillation tests will be simulated. First, the frequency ramp test is performed. Due to disturbances in the power system, the frequencies of the voltage or current signals may vary linearly. To further compare the SIFE (Method 4) and the proposed method (Method 5), the following signal was used for testing: where R 1 is the slope of the fundamental frequency and is set to 1 Hz/s. f varies linearly from 45 Hz to 55 Hz. The other parameters are identical to those in Section 4.3. The test results are shown in Table 3, where h = 1 represents the fundamental wave.
Where the maximum TVEs, FEs and RFEs are the average values of the errors. In [34], the TVE, FE and RFE limits for M-class PMU under the condition of frequency ramp are 1%, 0.01 Hz and 0.2 Hz/s, respectively. From Table 3, the max TVEs of Method 4 and Method 5 are relatively close, but the max FEs and RFEs of Method 4 are much greater than Method 5. This is because Method 4 uses different center frequencies to design the measurement filters, thus the corresponding filters are selected according to the actual frequencies of the signals. The value of the center frequency is related to the fundamental frequency band and is often taken as the midpoint of the frequency band. However, this relationship is valid only when the fundamental frequency varies linearly. Based on (29), s(t) contains the squared term of time. The fundamental and harmonic frequencies change nonlinearly with time, thus the center frequency is no longer located at the midpoint of the fundamental frequency band. Therefore, the estimation errors of frequency and ROCOF are increased. In contrast, Method 5 does not complete the measurement through filters, but introduces the AR model, which improves the frequency resolution and resistance to DI interference. It achieves the accurate estimation of harmonic phasors in multifrequency signals. Next, the harmonic oscillation test is performed. The amplitude and phase of harmonics may oscillate due to load variations or subsynchronous oscillations. To further verify the effectiveness of the proposed method, the following signal was used for testing: where s m (t) = 0.1cos(2πf m t) and f m is the modulation frequency, which is set to 2 Hz. Equation (30) is equivalent to the mathematical expression of the modulated signal. Other parameters are same as in Section 4.3. The test results are shown in Table 4. As shown in Table 4, the errors of both Methods 4 and 5 increase under the condition of harmonic oscillation, with Method 4 showing a higher error increase. In [34], the TVE, FE and RFE limits are 3%, 0.06 Hz and 2.3 Hz/s, respectively, for P-class PMUs with a reporting rate of 50 Hz. Both Methods 4 and 5 are able to satisfy the TVE standard. For the estimation of synchrophasor and low harmonic phasors, Method 5 is able to satisfy the FE and RFE standard, but Method 4 cannot. In addition to the reasons mentioned above, under the harmonic oscillation condition, Method 4 requires extensive simulations to determine the parameters of each harmonic phasor. Moreover, when the fundamental frequency is 50 Hz, Method 4 cannot obtain results identical to the real value. In summary, Method 5 has a better accuracy and interference immunity compared with other methods.

Experimental Signal Test
To evaluate the practical value of the proposed method, we built a test platform in laboratory conditions as shown in Figure 9. The experimental signal is first generated using a real-time digital simulation system and further processed by a power amplifier and a signal acquisition device. The sampled data are input to the PC to obtain estimates of harmonic phasors, frequencies and ROCOFs.
Method 5 has a better accuracy and interference immunity compared with other methods.

Experimental Signal Test
To evaluate the practical value of the proposed method, we built a test platform in laboratory conditions as shown in Figure 9. The experimental signal is first generated using a real-time digital simulation system and further processed by a power amplifier and a signal acquisition device. The sampled data are input to the PC to obtain estimates of harmonic phasors, frequencies and ROCOFs.
The steady-state test was first performed. The signal parameters were set as shown in Table 5, where the phases of the frequency components were set as random numbers uniformly distributed between (0, 2π). The accuracy of the estimation results is measured using the following indexes: where x and x are the true and estimated values of the parameters, respectively. To simulate the environment noise, 40 dB SNR and 60 dB SNR of noise are added to the experimental signal. The signal is estimated using Methods 1 to 5. Figure 10 shows the estimation errors for different methods.   The steady-state test was first performed. The signal parameters were set as shown in Table 5, where the phases of the frequency components were set as random numbers uniformly distributed between (0, 2π). The accuracy of the estimation results is measured using the following indexes: where x andx are the true and estimated values of the parameters, respectively. To simulate the environment noise, 40 dB SNR and 60 dB SNR of noise are added to the experimental signal. The signal is estimated using Methods 1 to 5. Figure 10 shows the estimation errors for different methods.  From Figure 10, it is clear that Method 1 cannot obtain estimates of the interharmonics and thus cannot be plotted in Figure 10. Methods 2 to 4 can detect all frequency components. However, since Method 2 uses a series of exponential functions to model the signal, the presence of noise makes the modeling difficult. Since the interference from adjacent frequency components, spectrum leakage and fence effect, the spectral peaks of interharmonics cause severe interference to the spectral peaks of harmonics. Although Methods 3 and 4 use a limited number of data points to obtain derivative estimates of the interharmonic phasors, they cannot restore distorted harmonic phasors, so there are still large errors. In contrast, Method 5 eliminates the interharmonic interference by building an AR model. The parameter estimates of Method 5 are almost the same as the real values even under large environment noise.
To analyze the situation under frequency deviation, the fundamental frequency is set to 52.5 Hz and 55 Hz, respectively. The signal difference (SD) between the real signal and the reconstructed signal is used as the performance index. SD = |s(t) − sr(t)|, where s(t) is the experimental signal and sr(t) is the reconstructed signal. Figure 11 shows the waveforms of the experimental signals at different fundamental frequencies and the SDs obtained by the proposed method. As shown in Figure 11, with the increase of the frequency deviation, the estimation errors of the proposed method are within an acceptable range. Therefore, for unknown signals, the proposed method can be applied to estimate harmonic phasor parameters. From Figure 10, it is clear that Method 1 cannot obtain estimates of the interharmonics and thus cannot be plotted in Figure 10. Methods 2 to 4 can detect all frequency components. However, since Method 2 uses a series of exponential functions to model the signal, the presence of noise makes the modeling difficult. Since the interference from adjacent frequency components, spectrum leakage and fence effect, the spectral peaks of interharmonics cause severe interference to the spectral peaks of harmonics. Although Methods 3 and 4 use a limited number of data points to obtain derivative estimates of the interharmonic phasors, they cannot restore distorted harmonic phasors, so there are still large errors. In contrast, Method 5 eliminates the interharmonic interference by building an AR model. The parameter estimates of Method 5 are almost the same as the real values even under large environment noise.
To analyze the situation under frequency deviation, the fundamental frequency is set to 52.5 Hz and 55 Hz, respectively. The signal difference (SD) between the real signal and the reconstructed signal is used as the performance index. SD = |s(t) − s r (t)|, where s(t) is the experimental signal and s r (t) is the reconstructed signal. Figure 11 shows the waveforms of the experimental signals at different fundamental frequencies and the SDs obtained by the proposed method. As shown in Figure 11, with the increase of the frequency deviation, the estimation errors of the proposed method are within an acceptable range. Therefore, for unknown signals, the proposed method can be applied to estimate harmonic phasor parameters.

Conclusions
As a result of the time-domain windowing of the signal, the spectrums of signal components with close frequencies overlap each other, causing severe spectral interference and reducing the estimation accuracy of harmonic phasors. Although the frequency resolution can be improved by using a longer observation window, this reduces the ability to track the signal dynamically. To address this problem, this paper proposes a harmonic phasor estimation method that considers DI interference. Firstly, the characteristics of the interfered spectral lines are used to determine whether spectral interference occurs in the signal. After that, data extrapolation is performed based on the signal sampling sequence to refine the spectrum and eliminate the interharmonic interference. Finally, the estimates of harmonic phasors, frequencies and ROCOFs are obtained. The proposed method can remove the interference of DIs while keeping the observation time unchanged. Simulation analysis indicates that the proposed method has a high accuracy under wideband noise, frequency deviation and dynamic conditions (frequency ramp and harmonic oscillation). A potential limitation is that the proposed algorithm does not consider the decaying DC component in the signal. Its presence may affect the estimation accuracy of harmonic phasors. Another limitation is whether the proposed spectrum refinement algorithm has a high computational burden. Therefore, further directions are the estimation of harmonic phasors containing decaying DC components and the computational burden.