Improved windowed phase difference method for frequency estimation of distorted power signals

Abstract

Frequency is one of the most significant indicators of power system, which needs accurate estimation to provide a reliable basis for monitoring, control, and protection. Owning to the efficiency, DFT is widely used in the frequency measurement of power systems. However, the inherent defects, that is, the spectral leakage and picket fence effect, decrease the accuracy of DFT. One of the variants, named the windowed phase difference method (WPD), addresses the problem using the phase variation of two windowed sequences in the frequency domain and achieves accurate frequency estimation. Suffering from the ignored negative component and contamination of noise, there are still undesired errors in WPD-based frequency estimation. To reduce the influence of the undesired errors, this paper first analyzes the systematic error of WPD introduced by spectrum leakage and proposes an improved WPD (IWPD) through systematic error compensation to achieve more accurate frequency estimations. Then, the windowing effect on IWPD-based frequency estimation of a noisy signal is studied by deducing the theoretical expression of frequency estimation variance. Finally, the proposed IWPD is validated by extensive computer simulations and practical experiments.

Introduction

The increasing proportion of renewable energy in power systems not only aggravates the imbalance between the power supply and demand but also deteriorates the power quality, because of the penetration of power electronic equipment and non-linear loads. ¹ As one of the products, frequency fluctuations bring great challenges for the safety and stable operation of power systems.^2,3 To address this issue, high-precision frequency estimation is needed to monitor the status of the power system, which can provide a reliable basis for power system control and protection.^4–6

In recent decades, several methods have been utilized to provide accurate frequency estimation, such as discrete Fourier transform (DFT),^7–9 Phase-locked loop, ¹⁰ Prony’s-method, ¹¹ wavelet transform, ¹² least-squared technique, ¹³ Kalman filter, ¹⁴ Taylor method, ¹⁵ demodulation technique, ¹⁶ and artificial neural network, ¹⁷ etc. Among them, DFT is widely used for its efficient implementation, the fast Fourier transform (FFT). With coherent sampling, DFT can estimate the frequency with extremely high accuracy. ¹⁸ However, acquiring the signal coherently always encounters intricate problems due to the frequency fluctuation in practice. Under non-coherent sampling, the inherent defects of DFT, for example, the spectral leakage and picket fence effect, will be produced, which are detrimental contributions for accurate frequency estimation. ¹⁹ Many tricks, such as the window functions and spectral correction algorithms, have been applied to DFT to reduce the errors introduced by those inherent defects, which can be called the windowed DFT methods. ²⁰ Weighting the discrete signal by window function in the time domain can effectively suppress the spectral leakage. ²¹ In addition, it is known that the accuracy of windowed DFT on frequency estimation is mainly determined by the window adopted. Thereby, various windows, Hanning window, Blackman-Harris (BH) window, Nuttall window, Maximum sidelobe Decay window (MDW), and the convolution windows, etc., are used. The tradeoff between the major lobe width and the sidelobe behaviors is unavoidable when selecting the window for frequency estimation. ²¹

A popular spectral correction method to eliminate the picket fence effect is the interpolation in the frequency domain. In which, the dual-spectrum lines and multi-spectrum lines are widely utilized thanks to their good performance on spectral correction. ²² However, getting the window’s frequency estimation interpolation polynomial involves solving high-order equations which may increase the computational burden. The Windowed Phase Difference method (WPD) utilizes the inherent relation between the estimated frequency and the phase difference of two groups of DFT. Compared to the interpolated DFT, since WPD avoids solving high-order equations, it has the merits of simple calculation and easy embedding. In Wen et al., ²⁰ a discrete phase difference correction algorithm based on a desirable sidelobe window is proposed for harmonic analysis, and the result validates that it can considerably reduce spectral leakage and harmonic interference. In Xia and Liu, ²³ the systematic error of the phase angle of unwindowed DFT coefficients in a single-phase system is analyzed. Although a “quasi-positive-sequence-DFT” is proposed in Xia and Liu ²³ to compensate for the systematic error, the frequency estimation precision is still limited without the help of windows.

As well known, in actual measurement, the inherent noise and the noise introduced by the sampling system, including the transducers, converters, conditioning circuits, etc. are unavoidable. Since the noise will decrease the accuracy of frequency estimation, evaluation of the frequency estimator’s performance under noisy conditions is mandatory. ²⁴ The windowing effect on the bias and variance of the DFT-based frequency estimation is investigated in Schuster et al. ²⁵ The influence of noise on the power system frequency estimation is studied in Wen et al.,^22,26 and Belega and Petri. ²⁷ The difference is that ²² concentrates on the triple-line interpolated DFT, while ²⁶ focuses on a double-line interpolated DFT with a triangular self-convolution window. In addition, the noise effect of the phase difference method without weighting effect (rectangle window) on frequency estimation is analyzed in Hwang and Liu. ²⁸ As mentioned before, the performance of the rectangle window (without windowing) is not so satisfying in frequency measurement under incoherent sampling. Therefore, it is necessary to analyze the systematic error of the WPD-based frequency estimator and the frequency variance under noisy conditions.

To the best of the authors’ knowledge, for the WPD, the error sourced from the negative part and the windowing effect on the measurement variance have not been investigated yet. Thus, to fill this gap, the systematic error of the WPD-based frequency estimator is first analyzed in this paper, and then an improved WPD (IWPD) is proposed to compensate for the systematic error. After that, the weighting effect of the window function on the noise is investigated, and the variance expression of frequency estimation is derived theoretically. At last, the computer simulations and experimental tests are carried out to verify the proposed systematic error and variance expressions, as well as the performance of the proposed IWPD.

The organization of this paper is as follows. In Section 2, the WPD for power frequency estimation is recalled. In Section 3, the systematic error of the WPD-based frequency estimator is analyzed and the IWPD is presented. Section 4 derives the theoretical variance and mean of WPD for frequency estimations under noisy conditions. Afterward, the simulations and experiments are carried out in Sections 5 and 6, respectively, to validate the proposed expressions and the IWPD. Finally, Section 7 discusses the conclusion.

Frequency estimation by WPD

A power system signal corrupted by harmonics acquired by a sampling frequency $f_s$ , can be presented as:

$$v(n)=\sum _{h=1}^H{A}_h\sin \left(2\pi n{f}_h/f_{\mathrm{s}}+{\theta }_h\right)$$ (1)

where $A_h$ , ${\theta }_h$ , and $f_h$ are the amplitude, phase and frequency of the $h\mathrm{th}$ harmonic, respectively, H is the highest harmonic order in the signal. $f_s=D{f}_{\mathrm{n}}$ , where D is an integer, $f_{\mathrm{n}}$ represents the nominal frequency of power systems.

To reduce the spectrum leakage caused by the truncation and incoherent sampling, an N-samples window function w(n) is employed to truncate v(n) in the range [0, N-1] in the time domain as:

$$x\left(n\right)=v\left(n\right)·w\left(n\right)$$ (2)

The DFT coefficients of $x(n)$ can be calculated by equation (3).

$$X(k)=\sum _{h=1}^H\frac{A_h}{2j}[e^{j{\theta }_h}{W}_f(k-k_h)-e^{-j{\theta }_h}{W}_f(k+k_h)]$$ (3)

where $j$ is the imaginary unit, $k=0,1,2$ , . . ., $N-1$ , is the spectral line index corresponding to the frequency $k\mathrm{\Delta }f$ , $\mathrm{\Delta }f=f_s/N$ means the frequency resolution, $k_h=f_h/\mathrm{\Delta }f$ represents the position where $f_h$ should be in the spectrum; $W_f$ is the DFT of $w(n)$ .

Note that $k_h$ can be represented by two parts due to the incoherent sampling by (4).

$$k_h=\frac{f_h}{\mathrm{\Delta }f}=l_h-{\lambda }_h$$ (4)

where $l_h$ is the integer part, ${\lambda }_h(-0.5\le {\lambda }_h\le 0.5)$ is the fractional part. This is to say that, the spectral line $k_h$ exactly corresponding to the $h\mathrm{th}$ harmonic may lie between two spectral lines, that is, the $l_h\mathrm{th}$ and the $(l_h-1)\mathrm{th}$ , or the $(l_h-1)\mathrm{th}$ and the $(l_h+1)\mathrm{th}$ .

As for the power system frequency $f_1$ , we can obtain that $k_1=l_1-{\lambda }_1=f_1/\mathrm{\Delta }f$ , where $l_1$ is the position of the maximum value in the spectrum. By neglecting the effect of harmonics and negative fundamental frequency, the complex spectrum of $l_1\mathrm{th}$ spectral line can be simplified as:

$$X\left(l_1\right)\approx \frac{A_1}{2j}{e}^{j{\theta }_1}{W}_f(l_1-k_1)$$ (5)

According to the principle of DFT, the phase $\phi$ of the spectral line $k_1$ can be expressed as:

$$\phi ={\theta }_1+\pi {\lambda }_1$$ (6)

It is worth noting that $\pi {\lambda }_1$ represents the phase difference caused by incoherent sampling, which is called incoherent sampling phase difference (ISPD). Besides, it can be seen from (4) that the ISPD is only determined by $\mathrm{\Delta }f$ which depends on $f_s$ and window length N.

When the signal is truncated by the same window $w(n)$ over the range $[M,N-1+M]$ , viz., there is a time delay $M$ between the two truncated groups, the DFT of new truncated samples $x_M(n)$ is calculated by:

$$\begin{array}{c}\multicolumn{1}{c}{X_M(k)=\sum _{h=1}^H\frac{A_h}{2j}[e^{j{\theta }_{M-h}}{W}_f(k-k_h)}\\ \multicolumn{1}{c}{-e^{-j{\theta }_{M-h}}{W}_f(k+k_h)]}\end{array}$$ (7)

where ${\theta }_{M-h}$ represents the initial phase of the new truncated samples shown as:

$${\theta }_{M-h}={\theta }_h+2\pi h{f}_{1}M/f_{\mathrm{s}}$$ (8)

It is straightforward that the difference of the initial phase between the two truncated samples is $2\pi f_{1}M/f_s$ .

Similarly, the complex spectrum of the $l_1\mathrm{th}$ spectral line of $X_M$ is calculated by:

$$X_M\left(l_1\right)=\frac{A_1}{2j}{e}^{j{\theta }_{M-1}}{W}_f(l_1-k_1)$$ (9)

The phase of the $l_1\mathrm{th}$ spectral line in its spectrum is:

$$\phi \prime ={\theta }_{M-1}+\pi {\lambda }_1={\theta }_1+\pi {\lambda }_1+2\pi f_{1}M/f_{\mathrm{s}}$$ (10)

As can be seen from the comparison of (6) and (10), the phase difference between the $l_1\mathrm{th}$ spectral lines of the two truncated groups can be simplified as:

$$\mathrm{\Delta }\phi =2\pi f_{1}M/f_{\mathrm{s}}$$ (11)

Thus, the frequency of the power system is estimated by the following equation:

$$f_1=\mathrm{\Delta }\phi f_{\mathrm{s}}/2\pi M$$ (12)

where $\mathrm{\Delta }\phi$ can be easily acquired by the values of $l_1\mathrm{th}$ spectral lines of the two sequences through equation (13).

$$\mathrm{\Delta }\phi =\angle [X_M\left(l_1\right)·X^{*}\left(l_1\right)]$$ (13)

where * denotes the complex conjugate, ∠ expresses the angle of a complex number.

The process of the WPD for power system frequency measurement is shown in Figure 1.

Figure 1.

Schematic diagram of WPD for power system frequency estimation.

Improved WPD by systematic error compensation

The operation in (14), which ignores the influence of harmonics and negative fundamental frequency components, leads to an error in the frequency estimation. To lower this undesired error, the IWPD based on systematic error compensation is presented in this section.

According to the expression of DFT coefficient given in (3), the $l_1\mathrm{th}$ spectral line can be calculated as:

$$X(l_1)=\sum _{h=1}^H\frac{A_h}{2j}[e^{j{\theta }_h}{W}_f(l_1-k_h)-e^{-j{\theta }_h}{W}_f(l_1+k_h)]$$ (14)

Comparing (14) and (5), it can be seen that, except for the component $(l_1-k_1)$ , the effect of other components, that is, $(l_1-k_h,2\le h\le H)$ , $(l_1+k_h,1\le h\le H)$ , is ignored. However, due to the frequency deviation, the interferences shown in (15) cannot be attenuated to 0, which produces an error on the phase estimation of the $l_1\mathrm{th}$ spectral line.

$$\begin{array}{c}\multicolumn{1}{c}{R\left(l_1\right)=\sum _{h=2}^H\frac{A_h}{2j}{e}^{j{\theta }_h}{W}_f(l_1-k_h)}\\ \multicolumn{1}{c}{-\sum _{h=1}^H\frac{A_h}{2j}{e}^{-j{\theta }_h}{W}_f(l_1+k_h)}\end{array}$$ (15)

The phase error $\varphi$ caused by $R(l_1)$ can be calculated by:

$$\varphi =\angle \left[X\left(l_1\right)\right]-\angle \left[\frac{A_1}{2j}{e}^{j{\theta }_1}{W}_f(l_1-k_1)\right]$$ (16)

Define ${\varphi }_M$ as the phase error caused by the residual of $X_M(l_1)$ , the systematic error of frequency estimation $f_e$ produced by the residual thereby can be expressed according to (12) as:

$$f_{\mathrm{e}}=\frac{({\varphi }_M-\varphi )f_{\mathrm{s}}}{2\pi M}$$ (17)

According to Zhang et al., ²⁹ the expression of $\varphi$ concerning the time index $n$ can be presented by equation (18).

$$\begin{array}{c}\multicolumn{1}{c}{\varphi \left(n\right)=\sum _{i=1}^H{a}_i\cos (2\pi (f_1+f_i)n/f_{\mathrm{s}}+{\varsigma }_i)}\\ \multicolumn{1}{c}{+\sum _{i=2}^H{b}_i\cos (2\pi (f_1-f_i)n/f_{\mathrm{s}}+{\psi }_i)}\end{array}$$ (18)

where $a_i$ , $b_i$ , ${\zeta }_i$ , and ${\psi }_i$ are the parameters determined by the modulation effect between $i\mathrm{th}$ harmonic and the fundamental. These parameters are not discussed since they have no affection for the proposed compensation method.

According to (18), it can be observed that the variation period of $\varphi$ is approximated to $1/2f_1$ when the second harmonic interference is missing in signal $v(n)$ , otherwise, the period will roughly be $1/f_1$ . As a result, summing up the periodically changing errors in one whole period can effectively offset the effects of errors, ³⁰ and then a frequency with systematic error compensation is acquired as:

$$\begin{array}{c}\multicolumn{1}{c}{f_1=\frac{f_s}{2\pi \mathrm{MD}}\sum _{n=0}^{D-1}[{\phi }_M\left(n\right)-\phi \left(n\right)]}\\ \multicolumn{1}{c}{=\frac{f_s}{2\pi \mathrm{MD}}\sum _{n=0}^{D-1}\mathrm{\Delta }\phi \left(n\right)}\end{array}$$ (19)

To improve the calculation efficiency, based on the sampling theorem, only necessary samples in sequence $\mathrm{\Delta }\phi$ are needed for the compensation, thereby the IWPD by systematic error compensation for power frequency estimation is given by equation (20).

$$f_1=\frac{\mu f_s}{2\pi \mathrm{MD}}\sum _{n=0}^{D/\mu -1}\mathrm{\Delta }\phi \left(\mu n\right)$$ (20)

where $\mu$ is a submultiple of $D$ . Usually, the value of $\mu$ is determined by the highest order $H$ of harmonic in the estimated signal. According to the sampling theorem, $f_{n}D/\mu$ should be greater than or equal to twice the highest frequency to be compensated, that is, $f_{n}D/\mu \ge 2(H+1)f_1$ , then we get $\mu \le f_{n}D/[2(H+1)f_1]$ . For example, in the case of $f_n=50\mathrm{Hz}$ , $D=64$ , $f_s=3.2\mathrm{kHz}$ , $f_1=49.9\mathrm{Hz}$ , and $H=5$ , the highest frequency needed to be compensated is 49.9 × 6 = 299.4 Hz, thus, $\mu$ should be less than 3200/(2 × 299.4)≈5.3, which means 4 is the suitable choice of $\mu$ . However, the smaller $\mu$ is, the better the compensation effect is, and the heavier the corresponding computational burden is. Considering that the amplitudes of some harmonics are relatively small, the contribution of these components to the systematic error can be neglected. By ignoring the effect of these weak components in the compensation, a more suitable $\mu$ can be determined by re-evaluating the maximum frequency needed to be compensated.

Noise analysis of WPD on frequency estimation

The analysis of noise is mandatory due to the inevitability of noise in actual measurement. Thus, the effect of noise on frequency estimation by WPD is studied in this section.

First, some of the merits of a window function are introduced. That is, normalized peak signal gain (NPSG), normalized noise power gain (NNPG), equivalent noise bandWidth (ENBW). The definitions of the mentioned merits are given by the following equations. ³¹

$$\mathrm{NPSG}=\frac{W_f\left(0\right)}{N}$$ (21)

$$\mathrm{NNPG}=\frac{{\left[W_f\left(0\right)\right]}^2}{N}$$ (22)

$$\mathrm{ENBW}=\frac{\mathrm{NNPG}}{\mathrm{NPS}{G}^2}$$ (23)

As we know, the sampled signal in practical measurement of the power systems can be expressed as:

$$y\left(n\right)=v\left(n\right)+z\left(n\right)$$ (24)

where $v(n)$ is a power signal, and $z(n)$ is the additive white Gaussian noise (AWGN) with zero mean and variance ${\sigma }^2$ . The signal-to-noise ratio (SNR) of the noisy power signal is defined as $\mathrm{SNR}=A^2/2{\sigma }^2$ , and A can be approximated to $A_1$ usually.

After being weighted by the window function $w(n)$ in the time domain, the windowed signal $x(n)$ can be presented as:

$$x\left(n\right)=y\left(n\right)·w\left(n\right)=(v\left(n\right)+z\left(n\right))·w\left(n\right)$$ (25)

while its DFT coefficients can be calculated by equation (26).

$$X\left(k\right)=V\left(k\right)+Z\left(k\right)$$ (26)

where $V(k)$ is the DFT of the windowed $v(n)$ , $Z(k)$ is the DFT of the windowed $z(n)$ . According to the characteristic of $z(n)$ , the mean and variance of $Z(k)$ can be calculated by the following equation. ²²

$$\mathrm{E}\left(Z\left(k\right)\right)=\sum _{n=0}^{N-1}\mathrm{E}\left(z\left(n\right)\right)w\left(n\right)e^{-j\frac{2\pi \mathrm{nk}}{N}}=0$$ (27)

$$\mathrm{Var}\left(Z\left(k\right)\right)=\mathrm{E}\left({\left[Z\left(k\right)\right]}^2\right)=N{\sigma }^2·\mathrm{NNPG}$$ (28)

where E() and Var() represent the expectation and variance, respectively. Besides, more equations about the statistical properties of $Z(k)$ can be acquired as follows.

$$\mathrm{E}\left(\mathrm{Im}\left(Z\left(k\right)\right)\right)=\mathrm{E}\left(\mathrm{Re}\left(Z\left(k\right)\right)\right)=\mathrm{E}\left(Z\left(k\right)\right)=0$$ (29)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Var}\left(\mathrm{Im}\left(Z\left(k\right)\right)\right)=\mathrm{Var}\left(\mathrm{Re}\left(Z\left(k\right)\right)\right)}\\ \multicolumn{1}{c}{=N{\sigma }^2·\mathrm{NNPG}/2}\end{array}$$ (30)

where Im() and Re() express the imaginary and real parts of a complex number, respectively.

Variance of the estimated frequency

After performing DFT calculation on two windowed signals having a time delay M, the phase difference $\mathrm{\Delta }\phi$ is acquired according to (13). For convenience, the index $l_1$ is omitted in the following derivation process.

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }\phi =\angle \left[(V_M+Z_M)(V^{*}+Z^{*})\right]}\\ \multicolumn{1}{c}{=\angle \left(V_M{V}^{*}\right)+\angle (1+\frac{Z_M}{V_M})+\angle (1+\frac{Z^{*}}{V^{*}})}\end{array}$$ (31)

Notice that the first term in (31) is the desired one, and the rest two terms related to $z(n)$ are interferences for phase difference estimation. Here we mark the latter two terms in (31) as $\alpha$ and $\beta$ as follows.

$$\{\begin{array}{c}\multicolumn{1}{c}{\alpha =\angle (1+\frac{Z_M}{V_M})}\\ \multicolumn{1}{c}{\beta =\angle (1+\frac{Z^{*}}{V^{*}})}\end{array}$$ (32)

Generally, the power of noise is neglected compared to the power of the signal when $\mathrm{SNR}>>1$ . In this case, the real parts of the $\alpha$ and $\beta$ are negligible. Also, it is easy to know that sin $a\approx a$ when $a$ approaches 0. Accordingly, we can obtain that:

$$\{\begin{array}{c}\multicolumn{1}{c}{\alpha \approx \angle (1+\mathrm{Im}\left(\frac{Z_M}{V_M}\right))\approx \mathrm{Im}\left(\frac{Z_M}{V_M}\right)}\\ \multicolumn{1}{c}{\beta \approx \angle (1+\mathrm{Im}\left(\frac{Z^{*}}{V^{*}}\right))\approx \mathrm{Im}\left(\frac{Z^{*}}{V^{*}}\right)}\end{array}$$ (33)

When $M$ is fixed, the desired term in (31) will be a constant, thus, the variance of $\mathrm{\Delta }\phi$ can be expressed as:

$$\mathrm{Var}\left(\mathrm{\Delta }\phi \right)=\mathrm{E}\left({\alpha }^2\right)+\mathrm{E}\left({\beta }^2\right)+2\mathrm{E}\left(\alpha \beta \right)$$ (34)

Meanwhile, the means of ${\alpha }^2$ , ${\beta }^2$ , and $\alpha \beta$ can be estimated using (9), (22), (30), which are given as below.

$$\mathrm{E}\left({\alpha }^2\right)=\mathrm{E}\left({\beta }^2\right)\approx \frac{2{\sigma }^2·\mathrm{ENBW}}{{\mathrm{NA}}_1^2}$$ (35)

$$\mathrm{E}\left(\alpha \beta \right)=\mathrm{E}(\mathrm{Im}\left(\frac{Z_M}{V_M}\right)·\mathrm{Im}\left(\frac{Z^{*}}{V^{*}}\right))$$ (36)

Since E $\left(\alpha \beta \right)$ can be considered in two scenarios, that is, $M\ge N$ and $1<M<N$ , the focus of Var $\left(\mathrm{\Delta }\phi \right)$ is on how to obtain E $\left(\alpha \beta \right)$ .

When $M\ge N$ , there is no overlap between $x_M(n)$ and $x(n)$ shown as Figure 2. In light of the uncorrelated property of AWGN, the mean of $\alpha \beta$ is 0. Thus, the variance of $\mathrm{\Delta }\phi$ is expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Var}\left\{\mathrm{\Delta }\phi \right\}=\mathrm{E}\left({\alpha }^2\right)+\mathrm{E}\left({\beta }^2\right)}\\ \multicolumn{1}{c}{=\frac{4{\sigma }^2·\mathrm{ENBW}}{{\mathrm{NA}}_1^2}}\end{array}$$ (37)

Figure 2.

Schematic diagram of WPD when $M\ge N$ .

While $1<M<N$ , the overlap between $x_M(n)$ and $x(n)$ are $N-M$ samples shown as Figure 3. However, it is difficult to evaluate the effect of the overlap due to the different weights of the window functions. To solve this problem, an index $\eta$ named overlap correlation (OC) is introduced, ³² which is defined as:

$${\eta }_{\rho }=\frac{1}{N·\mathrm{NNPG}}\sum _{n=0}^{\rho N-1}\left[w\left(n\right)w(n+(1-\rho )N)\right]$$ (38)

Figure 3.

Schematic diagram of WPD when $1<M<N$ .

where $\rho =(N-M)/N$ is the overlap rate of the two windows.

By applying OC, we can easily obtain the E $\left(\alpha \beta \right)$ as:

$$\mathrm{E}\left(\alpha \beta \right)\approx -\frac{2{\sigma }^2·\mathrm{ENBW}·{\eta }_{\rho }}{{\mathrm{NA}}_1^2}$$ (39)

Then, the variance of phase difference $\mathrm{\Delta }\phi$ is calculated as:

$$\mathrm{Var}\left\{\mathrm{\Delta }\phi \right\}=\frac{4{\sigma }^2·\mathrm{ENBW}·(1-{\eta }_{\rho })}{{\mathrm{NA}}_1^2}$$ (40)

From the comparison of (40) and (37), it can be found that (37) is the result of (40) when ${\eta }_{\rho }=0$ $(M\ge N)$ . Thus, Var $\left(\mathrm{\Delta }\phi \right)$ can be expressed regardless of the overlap as (40).

According to the relationship between frequency and phase difference, the variance of frequency estimation is obtained as:

$$\mathrm{Var}\left\{f_1\right\}=\frac{{f_s}^2·\mathrm{ENBW}}{2{\pi }^2{M}^{2}N·\mathrm{SNR}}·(1-{\eta }_{\rho })$$ (41)

Based on (41), it can be seen that the frequency variance of WPD is in direct proportion to ENBW. Meanwhile, it is inversely proportional to SNR and N. What is more, the variance will be significantly reduced while M increases.

Mean of the estimated frequency

To analyze the windowing effect on the mean of frequency estimations, (31) is rewritten as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }\phi =\angle \left[(V_M+Z_M)(V^{*}+Z^{*})\right]}\\ \multicolumn{1}{c}{=[\angle \left(V_M{V}^{*}\right)+\angle (1+Z_M/V_M)+\angle (1+Z^{*}/V^{*})]}\end{array}$$ (42)

The mean of phase difference $\mathrm{\Delta }\phi$ is estimated as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{E}(\mathrm{\Delta }\phi \prime )=\mathrm{\Delta }\phi +\angle (1+\frac{\mathrm{E}\left(Z_M\right)}{V_M}+\frac{\mathrm{E}\left(Z^{*}\right)}{V^{*}}+\frac{\mathrm{E}\left(Z_M{Z}^{*}\right)}{V_M{V}^{*}})}\\ \multicolumn{1}{c}{=\mathrm{\Delta }\phi -\frac{2{\sigma }^2·\mathrm{ENBW}}{{\mathrm{NA}}_1^2}·{\mu }_{\rho }·\sin \left(\mathrm{\Delta }\phi \right)}\end{array}$$ (43)

Then the mean of the estimated frequency is obtained as:

$$\mathrm{E}\left({f\prime }_1\right)=f_1-\frac{f_s·\mathrm{ENBW}}{2\pi \mathrm{NM}·\mathrm{SNR}}·{\mu }_{\rho }·\sin \left(\mathrm{\Delta }\phi \right)$$ (44)

When $M<N$ , there will be a bias for frequency estimation caused by the second term of (44) which is related to the overlap. That means for an unbiased measurement, the overlap should be 0 $(M\ge N)$ . Generally, compared to the systematic error shown in (17), the bias caused by the overlap is negligible when SNR is large enough, thus the mean of frequency estimations under noisy conditions is approximated as:

$$\mathrm{E}\left({f\prime }_1\right)=f_1$$ (45)

Simulations tests and results

This section contains the simulation analysis to validate the proposed systematic error, variance, and mean expressions of WPD for frequency estimation, as well as the performance of the proposed IWPD.

Simulation of systematic error of frequency estimation

To verify the correctness of the proposed systematic error expression (17), a pure voltage signal distorted by harmonics is employed. The sampling frequency $f_s$ is 3.2 kHz. The parameters of the harmonic content are shown in Table 1. The influence of the delay $M$ and frequency deviation on the systematic error of frequency estimation are considered.

Table 1.

Parameters of harmonics.

$h$	1	2	3	4	5
$A_h(V)$	220	2	15	3	7
${\theta }_h(\mathrm{rad})$	1.15	0.7	0.1	2.75	1.45

Effect of time delay on systematic error

Figure 4 depicts the variations of frequency versus different delay $M$ . The simulated conditions are set as the fundamental frequency $f_1=49.8\mathrm{Hz}$ , the time delay $M$ varies from 1 to 256, with steps of 3. Length $N$ of the adopted Blackman-Harris window (noted as BH) and 4-term 3-order Nuttall window (noted as Nuttall) is 128. Moreover, the result of the pure fundamental wave is also provided to illustrate the effect of harmonic on the systematic error of frequency estimation.

Figure 4.

The effect of delay on systematic error of frequency estimation by WPD.

The agreement between theoretical and simulated systematic error shown in Figure 4 confirms the correctness of the proposed expression of systematic error of frequency estimation. Besides, it can be seen that as the time delay $M$ increases, the overall trend of systematic error becomes increasingly smaller. The influence of harmonics is significant, both for the system error change period and the system error range. The intervals between two adjacent minima, that is, $P$ and $P_h$ , are approximated to 64 and 32 which correspond to the cases with and without harmonics, respectively. 64 samples per period means that the variation frequency is 50 Hz while 32 corresponds to 100Hz, and this validates the statement: variation period of systematic error is approximated to $1/2f_1$ when there is no second harmonic interference, otherwise, the period will roughly be $1/f_1$ .

Effect of frequency deviation on systematic error

The influence of frequency deviation on the systematic error of the WPD-based frequency estimator is evaluated by changing the fundamental frequency $f_1$ of the signal shown in Table 1 from 49.5 to 50.5 Hz by step. Time delay $M$ is fixed to 10, and window length $N$ is set to 128 and 256, respectively. Besides, the Blackman-Harris window, 4-term 3-order Nuttall window, Hanning window, and 4-term MDW are also employed. The theoretical and simulation systematic errors of the fundamental frequency are compared in Figure 5.

Figure 5.

The comparison of theoretical and simulation systematic errors of fundamental frequency estimation with different frequency deviations.

From Figure 5, the consistency between theoretical and simulated systematic errors validates the effectiveness of the proposed expression. ¹⁷ Other than this, the following conclusions can be summarized as (1) the systematic error will be reduced as the frequency deviation decreases. (2) An increment of window length can decrease the systematic error.

Performance of the IWPD

To evaluate the performance of the IWPD, the signal shown in Table 1 with $f_1=49.7\mathrm{Hz}$ is adopted. The simulation settings are $f_s=1.6\mathrm{kHz}$ , $M=1$ , $N=64$ , Blackman-Harris window is adopted. The results of the estimated $f_1$ by the IWPD based on different $\mu$ are drawn in Figure 6. In addition, the I2pDFT ²¹ and I3pDFT ³³ are employed to make a comparison.

Figure 6.

Performance of systematic error compensation for frequency estimation.

In Figure 6, the frequency estimations, sliding along $n$ , obtained by the IWPD based on $\mu$ =1, 2, 4, 8, and 16 are compared to the WPD. As can be seen, the estimations of IWPD are more accurate than those of WPD. Moreover, a greater $\mu$ generally results in a larger error, but it implies a lower computational burden. Thus, there is a tradeoff between accuracy and efficiency when using the IWPD. As for the comparison between IWPD, I2pDFT, and I3pDFT, it is noteworthy that, although the IWPD provides less accurate results than I2pDFT and I3pDFT with $N=128$ , it only requires 81 $(\mu =16)$ , 89 $(\mu =8)$ , 93 $(\mu =4)$ , 95 $(\mu =2)$ , and 96 $(\mu =1)$ samples which are less than 128 needed by I2pDFT and I3pDFT.

Validation of variance and mean of frequency estimation

To assess the correctness of the proposed variance and mean expressions shown in (41) and (45), simulations applied on a single-tone frequency voltage signal in the presence of the zero-mean additive white Gaussian noise are carried out. The parameters of the single-tone frequency signal are $A_1=220\mathrm{V}$ , phase ${\theta }_1=0.1\mathrm{rad}$ . Without loss of generality, each result is based on 5000 independent runs.

Effect of SNR

First, simulations are performed with noise at different SNRs ranging from 10 to 100dB, while the frequency of sine-wave is set as $f_1=49.9\mathrm{Hz}$ and 50.3 Hz. The length $N$ of the adopted Hanning window, Blackman-Harris window, 4-term 3-order Nuttall window, and 4-term MDW is 128. The ENBW of the adopted windows are 1.500, 2.004, 2.125, and 2.310, respectively. The time delay $M$ is set to 64 and 150, which correspond to the cases of $M<N$ and $M\ge N$ , respectively. The variance and mean results are shown in Figures 7 and 8. Meanwhile, the variances of case $M=150$ are compared to its Cramer-Rao Lower Bound (CRLB) of an unbiased estimator. ³⁴

Figure 7.

Variances and means of frequency estimations versus different SNRs: (a) variances of frequency, $M<N$ and (b) means of frequency, $M<N$ .

Figure 8.

Variances and means of frequency estimations versus different SNRs: (a) variances of frequency, $M>N$ and (b) means of frequency, $M>N$ .

Figure 7 shows that the variances and means of $f_1$ in the simulations match their theoretical values well. It can be found in Figures 7(a) and 8(a) that the variance is completely inversely proportional to the SNR and is independent of the frequency deviation. The difference between the variances based on different windows under the same simulation conditions depends mainly on the ENBW of the adopted window, as can be seen in (41).

As for the mean of the frequency estimations shown in Figures 7(b) and 8(b), the bias reaches 0.003 Hz only when the SNR is lower than 20 dB. Moreover, we can see that the increase in SNR leads to a reduction of the bias. In detail, the bias can be lower than 0.0004 Hz when SNR>40, which can be seen in the pictures.

Effect of the delay M

As a key parameter of WPD, it is important to evaluate the effect of time delay $M$ on the frequency variance. To this end, simulations are performed with window length $N=256$ , while $M$ increases from 1 to 512, with steps of 7. The settings of simulations include two cases, that is, $M<N$ and $M\ge N$ . The employed single-tone frequency signal is generated with $f_1=49.9\mathrm{Hz}$ , $A_1=220\mathrm{V}$ , and phase ${\theta }_1=0.1$ rad, while it is sampled under $f_s=3200\mathrm{Hz}$ . SNRs of the test signals are set to 53 and 67 dB, respectively. The frequency variances based on the Hanning window, Blackman-Harris window, 4-term 3-order Nuttall window, and 4-term MDW are demonstrated in Figure 9.

Figure 9.

Variances of frequency estimates versus different delay $M$ and SNRs.

As shown in Figure 9, higher SNR produces lower frequency variance, and it shows good agreement between the theoretical and simulation variance results versus different $M$ . It is worth noting that when $M<N$ , the frequency variance decreases faster with the increase of $M$ than that of $M\ge N$ . This trend is mainly determined by the variation of OC. In addition, the differences in the results based on different windows are caused by ENBW.

Experimental tests and results

In this section, practical tests are carried out through a PC-based platform shown in Figure 10 to evaluate the performance of the IWPD. The pure sinusoidal signals are generated by calibrator FLUKE 6105A, conditioned by the preprocessing circuit, sampled by the PCI-6221 data acquisition board, and then stored and processed in the PC. The parameters of the signal are, $A_1=220\mathrm{V}$ , and ${\theta }_1=0$ rad. The sampling frequency is set to 6.4 kHz, $N=256$ , $D=128$ , $n_1=1$ , and $M=1$ . The results of the practical experiments based on the Blackman-Harris window and 4-term 3-order Nuttall window are described in Figure 11. Also, the I2pDFT and I2pDFT with length $N$ 512 are adopted to make a comparison. The frequency estimations are based on 2000 independent measurements.

Figure 10.

Hardware platform of the experiments.

Figure 11.

Performance of the IWPD: (a) Blackman-Harris window and (b) 4-term 3-order Nuttall window.

As shown in Figure 11, $\mu =1$ , 4, 16, 32, and 64 are adopted. For the pure signal, the source of systematic error is the negative fundamental which leads the period of the phase error roughly to be $1/2f_1$ with respect to $n$ . In this case, $\mu =64$ cannot fulfill the requirement to compensate for the systematic error. As for the other values of $\mu$ , it can be seen clearly that $\mu =32$ is enough to complete the compensation task, and the increment of $\mu$ does not raise the performance of the compensation. The same situation can be found both in Figure 11(a) and (b). Compared to I2pDFT and I3pDFT, IWPD offers the same accuracy level of frequency estimations, but IWPD requires fewer samples, for example, I2pDFT and I3pDFT need 512 samples while IWPD uses 321, 353, 369, 381, and 384 samples, corresponding to $\mu =64$ , 32, 16, 4 and 1, respectively.

Conclusion

The performance of WPD on the power frequency estimation is affected by the systematic error and the noise inevitable in practical measurement. This paper firstly analyzed the systematic error of frequency estimation based on WPD, and then the IWPD through systematic error compensation is proposed. After that, the variance of frequency estimation under noisy conditions is studied by introducing the ENBW and OC of the window function. The IWPD and the proposed expressions are verified by the simulations and practical measurements. Some significant findings in this paper can be drawn as follows:

• (1) The proposed IWPD overpasses the WPD with the same adopted window, and its accuracy level is related to the value of $\mu$ , that is, the greater the $\mu$ , the lower the accuracy level.

• (2) According to the derived variance expression, it can be known that the frequency variance is inversely proportional to SNR, whereas it is proportional to ENBW of the adopted window.

Although the error compensation improved the accuracy of the WPD-based frequency estimations, the samples needed will be increased simultaneously. Therefore, error compensation with fewer samples can be the focus of future work.