The Influence of Power Network Disturbances on Short Delayed Estimation of Fundamental Frequency Based on IpDFT Method with GMSD Windows

Abstract

This paper presents an application of the IpDFT spectrum interpolation method to estimate the fundamental frequency of a power waveform. Zero-crossing method (ZC) with signal prefiltering was used as a reference method. Test models of disturbances were applied, based on real disturbances recorded in power networks, including voltage harmonics and interharmonics, transient overvoltages, frequency spikes, dips and noise. It was determined that the IpDFT method is characterized by much better dynamic parameters with better estimation precision. In an example, in the presence of interharmonics, the frequency estimation error was three times larger for the reference method than that for the IpDFT method. Furthermore, during the occurrence of fast transient overvoltages, the IpDFT method reached its original accuracy about three times faster than the ZC method. Finally, using IpDFT, it was possible to identify the type of disturbances: impulsive, step changes of frequency or voltage dips.

1. Introduction

Monitoring the frequency of voltage waveforms in power supply networks is one of the most important parameters from the point of view of, e.g., power quality measurements [1,2] and power network protection automation devices, especially those occurring in inverters used in the photovoltaic industry [3,4,5]. In the first industry, incorrect estimation of the fundamental frequency of the waveform results in a so-called spectrum leakage, and consequently, in an increased uncertainty of measurement of spectral parameters such as harmonics and interharmonics. On the other hand, incorrect estimation of the fundamental frequency in the second industry results in large power losses in the case of frequency mismatch between two different power systems and their phase synchronization [6]. This is important both in the case of cooperation among large power systems and also when connecting a large number of micro photovoltaic power plants to the power grid [4]. Inverters, which are a key component of such systems, must be equipped with disturbance-tolerant algorithms that estimate the instantaneous frequency of the electrical grid. According to [3], the requirements for power generation modules require the use of the so-called limited frequency sensitive mode at overfrequency (LFSM-O) and limited frequency sensitive mode at underfrequency (LFSM-U), which consist in a dynamic change of the generated energy as a function of changes in the grid frequency. Moreover, the aim of the protection automatics in photovoltaic systems is to protect the systems against unwanted switching between grid-connected and islanding operating modes [4,7]. For this purpose, the rate-of-change-of-frequency (ROCOF) method is used, among others. This method requires calculating the frequency over a number of cycles and comparing with a specified trigger threshold [3,8].

There are many algorithms for estimating the fundamental frequency of power waveforms. These methods can be divided into time, filter, frequency and other methods (e.g., time─frequency methods).

One of the basic temporal methods is the zero-crossing method. This method, due to its sensitivity to noise found in real power waveforms, is supplemented by prefiltering or regression least mean squares (LMS) modeling methods [9]. These models [10,11] and their extended least squares (ELS) extensions [12] are also used independently of the zero-crossing method.

Filtering methods include the use of Kalman filters [12,13,14,15]. In the ROCOF method, Kalman filters have been used, for example, to reduce the nondetection zone (NDZ) [8]. Furthermore, a combination of Kalman filtering and the least squares (LS) method is used to improve the dynamic properties of accurate frequency estimation [16].

Spectral methods provide the estimation of the fundamental frequency of the power waveform based on the spectrum processing obtained by discrete Fourier transform (DFT). In order to improve the accuracy of the spectrum estimation, its interpolation IpDFT (which will be discussed later in this paper) or other modifications, such as the so-called smart DFT (SDFT) with the additional use of the complex-valued least squares (CLS-SDFT) method [17], are used. The DFT transformation is also used repeatedly [18]. Spectrum interpolation can also be based on non-Fourier methods, e.g., Prony’s method to model the signal using exponentially damped sinusoids [19] or the Aboutanios and Mulgrew’s (HAM) method [20] and the cooperative weighted least squares (WLS) method representing HAM-WLS [20]. Spectrum estimation is also performed using time─frequency analyses e.g., discrete wavelet transform (DWT) [21] or Wigner─Ville transform [22].

Finally, the methods can be applied to the analysis of single- and three-phase waveforms [23,24,25].

IpDFT spectrum interpolation methods deserve special attention; the first one was developed by Rife and Vincent in 1970 [26], but their rapid development has occurred over the past two decades. These methods combine the advantageous properties of using nonrectangular time windows, a FFT algorithm and further processing of the resulting DFT spectrum to reduce errors due to the discrete nature of the resulting spectrum. The use of such time windows as Hann (often referred to as Hanning windows), Hamming, Kaiser, Chebyshev and many others, allows a significant reduction of the phenomenon of spectrum leakage, originating from harmonics and interharmonics in the conditions of nonsynchronous (noncoherent) sampling with the period of the measured signal. On the other hand, the use of the FFT algorithm in a modern DSP system with a signal processor that is optimized for fast execution of such an algorithm makes the computation time much smaller than the signal measurement time, while maintaining a low cost of the system. The time window method combined with the FFT algorithm belongs to classical Fourier analysis and has a reputation for being fast but not very accurate due to the so-called picket fence effect. This effect is the cyclic nature of the frequency and amplitude estimation errors when these parameters are determined from the local maximum of the raw DFT spectrum. Then, the maximum error of frequency estimation is half of the DFT computational resolution (i.e., ±0.5/NT—half the inverse of the measurement time) and does not depend on the applied data window. For typical time windows, the maximum amplitude estimation errors caused by picket fence effect range from −9% to −22% [27]. The essence of IpDFT methods is to further process the raw DFT spectrum to eliminate errors caused by the discrete nature of the DFT spectrum (picket fence effect) without significantly increasing the computation time. In this way, IpDFT methods combine the beneficial properties of classical Fourier analysis methods (high speed and low cost) and non-Fourier methods (usually much more accurate, but much slower). The latest strands of IpDFT methods are those that take into account the conjugate component in the spectrum, developed especially for applications with short measurement times—of the order of a few signal periods [28,29,30,31]. Their latest version [31] is further discussed in a more detailed way in Section 3 and its basic flowchart for a real-time grid-monitoring system is presented in Figure 1.

The method presented in [31] and in Figure 1 has the following advantages:

• It allows for the estimation in short measurement times with high accuracy because the conjugate component in the spectrum (i.e., a component with a negative frequency resulting from the mathematical properties of the Fourier transform) is taken into account during derivation of the estimating formula;
• The estimation formula used in the calculations only applies a GMSD window, the FFT algorithm and a simple interpolation formula, which involves three points of the FFT-resulting spectrum;
• It allows for a cheap implementation, because the calculation time depends mainly on the calculation time of the FFT algorithm, hence modern cheap signal processors and microcontrollers optimized for FFT are applicable;
• It can be used for signals with a large THD coefficient, because GMSD windows significantly eliminate its influence on the accuracy of the estimation.

The basic comparison of the presented method with parametric methods is presented in

Table 1. Qualitative comparison of the IpDFT method with GMSD windows (which takes into account the conjugate component during derivation of the estimation formula) vs. parametric methods depending on the normalized frequency λ₁ (the number of signal periods in the measurement window).

Measurement Window Duration	IpDFT-Based Frequency Estimation Methods		Parametric Methods
	IpDFT Method Which Takes into Account Conjugate Component	IpDFT Methods Which Do Not Take into Account Conjugate Component	Parametric Methods: Prony LS, TLS (Total Least Squares), ESPRIT
Very short window (λ1 < 1)	Applicable only for low level of noise	Not applicable due to great number of systematic errors	Applicable, especially for high resolution methods (e.g., ESPRIT)
Short window (λ1 ≈ 1…3)	Applicable Good accuracy (systematic errors below the level caused by noise) Cheap implementation	Not applicable due to great number of systematic errors	Applicable, but expensive in practice due to the great number of calculations required (even by a few orders for N > 1000)
Long window (λ1 >> 3)	Applicable Good accuracy Cheap implementation	Applicable Good accuracy Cheap implementation	Applicable Good accuracy Expensive implementation

. More quantitative details are provided, e.g., in [29,30]. It shows that the presented method has very good properties for short, delayed estimation of the fundamental frequency, as for λ₁ ≈ 1…3 it still maintains its high accuracy and at the same time it does not require many calculations.

The paper focuses on the influence of the most important disturbances on the estimation results. It is the first systematical study of this aspect of the method [31]. The models of grid disturbances are presented in Section 2 and the interpolated discrete Fourier transform is summarized in detail in Section 3. The important role of the proper choice of the method’s parameters is described in detail in Section 4, as well as a short description of the characteristics of the zero-crossing (ZC) method used in the paper as a reference method. The simulations performed in the MATLAB software environment are presented in Section 5, and they include detailed information about the features of the method when a grid signal contains various disturbances, as is usually the case. The conclusions are presented in Section 6.

2. Disturbances and Its Models in Power Networks

2.1. Disturbances in Power Networks and Power Quality

The voltage electrical waveform in power networks is by definition a strictly deterministic signal. It is modeled by a single sinusoidal waveform with nominal rms value, symmetric with respect to the 0 V level and nominal frequency. In three-phase networks, these waveforms are additionally shifted in individual phases by 120°.

In fact, this waveform additionally contains disturbances of deterministic and nonstationary character. The first group includes cyclic or permanent disturbances that are specific to the supply network, for example, the continuous presence of the 5th harmonic, long-lasting voltage underdeviation or asymmetry. The second group includes disturbances appearing every indefinite time, or whose character evolves in time in an unpredictable manner, e.g., voltage overvoltages caused by lightning, momentary changes of the fundamental frequency under the influence of rapid changes of the load or oscillatory disturbances occurring during activation of capacitor banks to compensate reactive power.

Some of the disturbances occurring in the supply voltage waveform have a direct impact on the power quality and on the operation and safety of electrical equipment. These include all kinds of voltage fluctuations, harmonics, interharmonics, frequency fluctuations and supply voltage asymmetry. The remainder at least affects the user comfort of electrical equipment. These disturbances can include, inter alia: the phenomenon of flickering caused by a series of changes in the rms value of the voltage of the light source supply, dips and interruptions in power supply. Definitions and acceptable levels of most of the mentioned disturbances are defined in relevant standards of power quality, such as EN501610 [32], and electromagnetic compatibility from the IEC 61000 series [33].

Based on the authors’ own research, various types of interference were acquired, which can be divided into three groups:

-

Temporary disturbances—of incidental nature;

-

Sustained disturbances—occurring in the power network for a longer period of time;

-

Mixed disturbances—containing the two above mentioned types of disturbances.

Figure 2 as well as Figure 3 and Figure 4 present examples of disturbances from the first group: transient overvoltage, voltage fluctuations and dips and instantaneous frequency changes, respectively.

The duration of the first one is about 3 ms and it has an oscillatory character with a frequency of about 4 kHz.

The second disturbance is an example of a voltage drop of about 2.4% and two consecutive dips of 12% and 24%.

The third disturbance represents changes in the instantaneous frequency of the voltage. These changes are often of a sustained disturbance nature, which include long-term deviations of the power network frequency within ±1% for 99.5% of the year [32]. However, dynamic fluctuations of the fundamental frequency are also encountered, caused by rapid changes in the power network load, which include the frequency fluctuation of 0.3% shown in Figure 3.

The second group of disturbances—sustained disturbances—often involves, among other things, the occurrence of harmonics and interharmonics, which are generated by nonlinear loads, in particular by voltage converters, available, for example, in any IT equipment. These disturbances are parameterized, for example, by so-called spectral power quality parameters, which include [34]:

-

Total distortion factor of harmonic subgroups—Equation (1),

$$THD{S}_U=\sqrt{{\displaystyle \sum }_{h=h_{min}}^{h_{max}}{\left(\frac{U_{sg,h}}{U_{sg,1}}\right)}^2}$$

(1)

where:

$$U_{sg,h}^2={\displaystyle \sum }_{k=-1}^1{U}_{C,\left(N\times h\right)+k}^2$$

(2)

$U_{sg,h}^2$—harmonic subgroup of the h order.

$U_{sg,1}^2$—first-harmonic subgroup.

$U_{C,l}$—rms value of the spectral voltage component of the l order.

N—number of periods of the fundamental component.

h—harmonic order.

h_min—minimum harmonic order: 2.

h_max—maximum harmonic order: 50.

-

Centered subgroups of interharmonics—Equation (3),

$$U_{isg,h}^2={\displaystyle \sum }_{k=2}^{N-2}{U}_{C,\left(N\times h\right)+k}^2,$$

(3)

where the interpretation of individual parameters is as for Equation (2).

Examples of harmonic disturbances registered in an office building are presented in Figure 5. The values of particular harmonic subgroups are shown in

Table 2. Parameters of measured harmonics from THDS_U = 10.3%.

Parameter	i—Harmonic Number
	1	3	5	7	9	11	13	15	17	19	21	23
$f_i=i·f_1\left[\mathrm{Hz}\right]$	50	150	250	350	450	550	650	750	850	950	1050	1150
$U_{sg,i \mathrm{rms} }\left[\mathrm{V}\right]$	225	0.6	15.2	14.5	1.2	6.4	4.1	0.8	4.6	3.6	0.9	0.6
$U_{sg,i \mathrm{rms}} [\%]$	100	0.29	6.74	6.47	0.56	2.91	1.82	0.37	2.04	1.65	0.41	0.26

. Deformation of the waveform in the area of minima and maxima is characteristic for this type of disturbances. Moreover, this type of disturbance usually has a quasi-stationary character and lasts for a minimum of several dozen periods or persists permanently. As a rule, the values of even harmonics are minimal, while the dominant ones are the 5th, 7th and 11th harmonics. The total distortion factor of the harmonic subgroups should remain below 8% for 95% of the week [32]. The mentioned example represents the exceedance of this value: 10.3%. When interharmonics are present, the distortion factor of the centered interharmonic subgroups has not yet been defined in the relevant standard.

Noise is also a disturbance that affects the distortion of the waveform spectrum. Figure 6 shows an example of recorded noise-like interference with a signal-to-noise (SNR) of 27 dB. Disturbances of this type are also usually of sustained character.

An example of the third group of disturbances, concerning a combination of both transient and sustained disturbances, is presented in Figure 7. It shows a harmonic disturbance with superimposed transient overvoltage of oscillatory nature.

2.2. Disturbance Modeling

Based on the real disturbances occurring in power grids and on the test procedures included in the standards for power quality and electromagnetic compatibility, as well as in the relevant documents related to the requirements imposed on inverters in PV microgrids, the following disturbance models are proposed. These models will allow us to investigate the influence of certain types of disturbances on the accuracy of the instantaneous frequency measurement.

2.2.1. Harmonics and Interharmonics

The parameters of the signal models representing the harmonic and interharmonic disturbances are presented in

Table 3. Harmonic parameters of the signal from THDS_U = 10.3%.

Parameter	i—Harmonic Number
	1	3	5	7	9	11	13	15	17	19	21	23
$f_i=i·f_1\left[\mathrm{Hz}\right]$	50	150	250	350	450	550	650	750	850	950	1050	1150
$U_{sg,i} \left[\mathrm{V}\right]$	225	0.6	15.2	14.5	1.2	6.4	4.1	0.8	4.6	3.6	0.9	0.6
${\phi }_i \left[°\right]$	0	0	0	0	0	0	0	0	0	0	0	0

and

Table 4. Interharmonic signal parameters.

Parameter	i—Interharmonic Number
	1	3.2	3.3	3.4	3.5	3.6	3.7	3.8	3.9
$f_i=i·f_1\left[\mathrm{Hz}\right]$	50	160	165	170	175	180	185	190	195
$U_{ig,i}\left[\mathrm{V}\right]$	230	1.5	2	4.5	12	17	4	3	2
${\phi }_i \left[°\right]$	0	0	0	0	0	0	0	0	0

, respectively. Due to the specificity of the power waveforms, harmonics are practically nonexistent and they have not been taken into account. The amplitudes of the individual harmonics and interharmonics have been selected based on actual tests. The THDS_U parameter exceeds the acceptable threshold defined in the standards—8% [32].

Graphical representations of the waveforms with harmonic and interharmonic disturbance models are shown in Figure 8a and Figure 8b, respectively.

2.2.2. Transient Overvoltages

The transient overvoltages were modeled using exponentially damped oscillatory disturbances. Two models were proposed: short duration and long duration. These models were expressed as (4) and (5):

$$s_{\mathrm{I}}(t)=2000{\mathrm{e}}^{-0.0583t}\sin (2\mathsf{\pi }{10}^{6}t)$$

(4)

$$s_{\mathrm{II}}(t)=500{\mathrm{e}}^{-0.00029278t}\sin (2\mathsf{\pi }{10}^{3}t)$$

(5)

In the first case, the oscillation frequency is equal to 1 MHz, the amplitude is 2000 V, and the duration after which the waveform is damped 100 times is 15 µs.

In the second case, these parameters are, respectively: 1 kHz, 500 V and 6 ms.

In the case of a short-duration disturbance, its amplitude was limited to 400 V due to the fact that the input systems of the measuring equipment are practically unable to transfer such a steep disturbance and in fact, the final waveform to be acquired has an amplitude of 400 V.

Figure 9a shows the signal with the disturbance expressed by (4), and Figure 9b shows the superimposed disturbance on the power waveform. Figure 10 shows the above for the disturbance defined by Equation (5).

2.2.3. Changes in Fundamental Frequency

The model of the signal containing fluctuations of the basic component of the power waveform consists of five waveforms with the rms value of 230 V and the following frequencies (6):

$$y(t)=\{\begin{array}{cc}42.50 \mathrm{Hz}& \mathrm{for} t<0.353 \left[\mathrm{s}\right]\\ 46.25 \mathrm{Hz}& \mathrm{for} 0.353<t<0.677 \left[\mathrm{s}\right]\\ 50.00 \mathrm{Hz}& \mathrm{for} 0.677<t<0.977 \left[\mathrm{s}\right]\\ 53.375 \mathrm{Hz}& \mathrm{for} 0.977<t<1.258 \left[\mathrm{s}\right]\\ 57.50 \mathrm{Hz}& \mathrm{for} 1.258<t<1.600 \left[\mathrm{s}\right]\end{array}$$

(6)

Each waveform contains 16 complete periods at the frequencies specified above.

2.2.4. Multicomponent Disturbances: Harmonics, Dips, Transient Overvoltages

Finally, a model consisting of harmonics, two transient overvoltages and three voltage dips was used to illustrate the response of the tested frequency measurement method in the presence of many different types of disturbances. Detailed parameters are shown in

Table 5. Multicomponent disturbance parameters.

Type of Disturbance	Parameter	Occurrence Time (s)
Harmonics	In accordance with Table 2	0–1.6
Transient overvoltages	In accordance with Formula (5)	0.6–1.2
Voltage dips	Urms = 230 V	0–0.7
	Urms = 215 V	0.7–0.9
	Urms = 181 V	0.9–1.1
	Urms = 162 V	1.1–1.3
	Urms = 230 V	1.3–1.6

. Throughout the test run, the fundamental frequency was constant at 50 Hz.

Figure 11 presents the signal with disturbances, as described in Table 5.

3. Interpolated Discrete Fourier Transform with the Use of GMSD Windows

The characteristics of the power network signal presented in Section 2 allow us to define the stationary part of this signal as:

$$y(t)={\displaystyle \sum _{m=1}^K{A}_m\sin (2\mathsf{\pi }{f}_{m}t+{\phi }_m)}$$

(7)

where A₁, f₁, and φ₁ are the parameters (amplitude, frequency and phase, respectively) of the fundamental component, and A_i, f_i, and φ_i are the harmonic and interharmonic parameters. Such a signal is found in many fields of science and technology and is often referred to as a multifrequency signal because it consists of multiple sinusoidal components. For obvious reasons, such a signal is further disturbed by the noise present in each measurement and a number of additional distortions, which for the case of the mains signal are described in Section 2. Rms values U_{i rms} of the sinusoidal components are obtained from the relation $U_{i \mathrm{rms}}=A_i/\sqrt{2}$. For a complete description of the power network signal, estimation of all A_i, f_i and φ_i parameters is required for all i, but the most important of these is the value of f₁, i.e., the frequency of the power network signal. It is nominally 50 Hz or 60 Hz, and in practice it is a value varying within the range allowed by legal regulations or exceeding this range in cases of control system failures in the process of electricity generation. Precise determination of the parameter f₁ determines the remaining parameters of the power network signal.

Many methods have been developed, parametric and nonparametric, allowing estimation of the parameters of the multifrequency signal modeled by Equation (7). Particularly noteworthy are spectrum interpolation methods (IpDFT), whose development and range of applications is constantly increasing. The application of IpDFT methods for the determination of f₁ is also studied in the context of power network signal parameter estimation for the following reasons:

-

In IpDFT methods, the computation time is much smaller than the signal measurement time, which significantly reduces the total estimation time of f₁ (the total estimation time is equal to the sum of the signal measurement time and the duration of calculations necessary to be performed after the measurement is completed);

-

IpDFT methods achieve high estimation accuracy of f₁, even with short signal measurement times, i.e., between one and three signal periods; this allows for a quick response in systems controlling power generation or disconnecting devices from the supply network;

-

In IpDFT methods, the computational complexity of the estimation algorithm is much lower than in parametric methods, which lowers the cost of the DSP system for determining the parameters of the power network signal; this is particularly important for small power generation systems (e.g., small photovoltaic installations) or for power quality monitoring by consumers.

The latest IpDFT methods are based on the complex values of the DFT spectrum obtained by the FFT algorithm and include in their solution a conjugate component, i.e., a component with a negative frequency, resulting from the mathematical properties of the Fourier transform. The use of the complex values of the DFT spectrum eliminates the drawbacks of earlier IpDFT methods, which used the DFT spectrum modules. On the other hand, the inclusion of the conjugate component allows for a significant reduction in signal measurement time, relative to earlier IpDFT methods, without increasing the systematic error values. The paper [31] presents the latest version of the IpDFT method taking into account the conjugate component and using the complex values of the DFT spectrum. This method uses windows of sin^m type, which are, as shown in [31], generalized MSD windows already used in earlier IpDFT methods. The IpDFT method presented in [31] uses GMSD windows, which are described by two, mutually equivalent, relations:

$$w_n={\displaystyle \sum _{h=0}^{H-1}{(-1)}^h{a}_h\cos \left[\frac{(2h+z)\mathsf{\pi }n}{N}+\frac{z\mathsf{\pi }}{2}\right]}, n=0,\dots ,N-1$$

(8)

$$w_n={\sin }^m\left(\frac{\mathsf{\pi }n}{N}\right)={\sin }^{2r+z}\left(\frac{\mathsf{\pi }n}{N}\right), n=0,\dots ,N-1$$

(9)

where w_n are the values of the GMSD window function for indices n = 0, …, N − 1, N is the number of samples of the GMSD window, m = 2r + z is a natural number, H = r + 1, z = 0 for even m and z = 1 for odd m, and a_h are appropriately chosen coefficients to ensure that the GMSD window has the largest sidelobes damping for a given H. Relationships (8) and (9) are equivalent, as shown in [31]. To derive the solution of the method in [31], (8) is used, and for practical use, the more convenient form (9).

The signal y(t) from (7) is processed in an A/D converter operating at a sampling rate f_s (T = 1/f_s is the sampling period) and N signal samples y_n = y(nT) are obtained at the A/D output. After introducing normalized frequencies λ_m = f_mNT in (bins) or (cycles in range (CiR)) (and in general, λ = fNT) the estimation of f₁ is equivalent to the estimation of λ₁. The use of λ in place of f and λ_m in place of f_m simplifies the mathematical model and allows a direct frequency reference of f₁ from the power network signal to the duration of the signal measurement (time window duration). For example, the value of λ₁ = 2.5 bin means that a time window of NT (N samples obtained every sampling period T) contains 2.5 periods of the power network signal.

Applying the time window (8) and (9) to the signal samples (7) reduces to multiplying y_n by w_n for all n and then the discrete-time Fourier transform (DtFT) is calculated, i.e., X_λ defined by:

$$X_{\lambda }=X(\lambda )={\displaystyle \sum _{n=0}^{N-1}{y}_n{w}_n{\mathrm{e}}^{-\mathrm{j}2\mathsf{\pi }n\lambda /N}} \mathrm{for} \mathrm{any} \lambda$$

(10)

The spectrum X_λ for natural values of λ is a commonly known discrete Fourier transform (DFT), which is usually calculated with the FFT algorithm. However, the IpDFT method from [31] requires knowledge of the spectrum X(λ) for λ = 0, 0.5, 1, 1.5, etc. (i.e., with a step of 0.5 bin), which is obtained by the FFT algorithm in so-called double zero padding technique [27].

The most important result in [31] is the dependencies on ${\hat{\lambda }}_1$ (and consequently on ${\hat{f}}_1={\hat{\lambda }}_1/NT$):

$${\hat{\lambda }}_1^2=l_1^2+\mathrm{Re}\left\{{\epsilon }_1\right\}$$

(11)

$${\epsilon }_1=(m+2)\frac{(m+2+4l_1)X_{l_1+1}+2m{X}_{l_1}+(m+2-4l_1)X_{l_1-1}}{4(X_{l_1+1}-2X_{l_1}+X_{l_1-1})}$$

(12)

where:

$$l_1=1, 1.5, 2, 2.5, \dots$$

(13)

This solution means that the spectrum interpolation procedure uses three DtFT spectrum points obtained with a step of 1 bin. The value of l₁ should be close to the value of λ₁, as this condition minimizes the variance of the estimator and hence the error due to noise in the signal.

The value ${\hat{\lambda }}_1$ differs from λ₁ due to the presence of systematic errors (which even occur for pure sine wave), the influence of harmonics and interharmonics, the presence of noise in a real measured signal and other signal disturbances. The determination of the impact of typical disturbances occurring in the power network on the obtained frequency estimation accuracy f₁ is the goal of the research presented in Section 5 in accordance with the methodology presented in Section 4.

4. Research Methodology

4.1. Assumptions for IpDFT

4.1.1. Measurement Time (Time Window Duration) NT

Based on the properties of the IpDFT method presented in [31], it was assumed in this study that the NT time of the power network signal measurement with frequency f₁ = 50 Hz is NT = 50 ms. This means that the normalized frequency λ₁ = 2.5 bin, i.e., that the measurement window (time window) covers 2.5 periods of the power network signal. The value of the measurement time NT = 50 ms was chosen in the study as a compromise value between aiming for the shortest possible measurement time and aiming for the highest possible accuracy. The IpDFT method of [31] allows accurate signal measurement for several times smaller measurement times, but then the effect of signal harmonics on the estimation accuracy increases, causing a significant reduction in accuracy. For such smaller measurement times (of the order of 1–2 signal periods) it is necessary, to maintain high accuracy, to use filtering that reduces the influence of harmonics, and this would degrade the dynamic properties of the method. As demonstrated by the studies presented in Section 5 adoption of NT = 50 ms (λ₁ = 2.5 bin) maintains high estimation accuracy without the need for additional filtering at the input of the method. In practice the time window will not cover exactly 2.5 signal periods due to the power network frequency variation f₁ and because of the characteristics of the A/D converters used, for which a change in sampling frequency f_s = 1/T is not possible continuously, but only among a limited set of values. However, the most important advantage of IpDFT methods is that there is no need to precisely synchronize the duration of the measurement window with the current frequency f₁ of the supply network. This is an essential feature of IpDFT methods that distinguishes them from methods with so-called coherent sampling. The adoption of NT = 50 ms is also the condition furthest from coherent sampling, since in coherent sampling, the duration of the measurement window must be equal to an integer multiple of the signal period. The lack of significant sensitivity of the IpDFT method tested in Section 5 to changing the duration of the time window is shown in Section 5.1, where the accuracy of the method was tested for NT = 50 ms and for two values differing by ca. ± 3%: NT = 48.44 ms and NT = 51.65 ms. In the other studies (Section 5.2, Section 5.3, Section 5.4, Section 5.5, Section 5.6 and Section 5.7), it was also found that a similar change in NT measurement time does not significantly affect the estimation results, and therefore only the case of NT = 50 ms was included in the presented graphs to maintain the readability of the presented figures.

4.1.2. Number of Signal Samples N and Sampling rate f_s = 1/T

For most studies in Section 5, the number of samples was assumed to be equal to N = 1024 and N = 2048. This means, with measurement time NT = 50 ms, that the sampling rate is, respectively f_s = 20,480 Hz and f_s = 40,960 Hz. The accuracy of the method for wide-spectrum signal disturbances (Section 5.3) is also shown for f_s = 2,621,440 Hz ≈ 2.62 MHz (for N = 2¹⁷), to show that the error for smaller values of f_s is caused by the fact that the conditions of Shannon’s sampling theorem regarding the minimum value of the sampling frequency are not satisfied. The increase in the error due to the failure to satisfy the conditions of Shannon’s sampling theorem is eliminated in practice by the use of an antialiasing filter at the input of the A/D converter and partly by the limited spectrum of the coupling systems between the network and the measuring system. With the limited frequency spectrum of the power network signal, it is possible to use smaller values of N and f_s. For the studies presented in Section 5.1, Section 5.2 and Section 5.5, the following are presented in addition to the two fundamental cases N = 1024 and 2048 (f_s = 20,480 Hz and 40,960 Hz providing NT = 50 ms) and the results for N = 128, 256 and 512 (f_s = 2560 Hz, 5120 Hz and 10,240 Hz providing NT = 50 ms). The presentation of results for five different values of N from 128 to 2048 provided a better representation of how the parameters N and f_s at the given NT affect the accuracy of the estimation. In many cases, low values of N and f_s allow to obtain satisfactory estimation accuracy under the condition of limited power network signal bandwidth (e.g., due to the use of antialiasing low-pass filter) and fulfillment of the condition of Shannon’s theorem for the minimum sampling frequency. It should also be remembered that increasing the number of samples of the signal N (i.e., increasing the sampling frequency of f_s at a constant NT measurement time value) reduces the effect of noise on the estimation accuracy, as shown in Section 5.7.

4.1.3. The m Parameter of GMSD Window

The basic parameter of the IpDFT method from [31] is, in addition to N and f_s, the m parameter of the GMSD window defined by (8) and (9). In the studies presented in Section 5 a fixed value of m = 3 was adopted as a compromise between the desire to suppress the sidelobes as much as possible (long-range spectral leakage) and to widen the main lobe width as little as possible. An increase in the value of m also results in an increase in the variance of the estimator, i.e., the error due to noise in the signal. The authors’ studies have shown that the value of m = 3 is the optimum case for many applications of power network parameter estimation and can be treated as a universal condition when there is no more specific other rationale for adopting a different value of m. A study of the effect of other values of m in the presence of the disturbances considered in Section 5, in the opinion of the authors, should be studied for individual target applications, which is beyond the scope of this paper, which aims to present the basic properties of the method for a wide class of applications.

4.1.4. The Effect of Signal Phase and Sliding Window

The phase φ₁ of the power network signal described by (7) has important implications on the accuracy of the estimation. In the case of continuous observation of the signal (Figure 12), as is the case in network parameter monitoring systems, the phase of the signal changes over time as the measurement window of duration NT includes the last N samples of the signal relative to the current instant of time (Figure 12b). Each such set of consecutive N samples can be referred to as a data frame (Figure 12c). As shown in [31], the computation time required in a typical DSP system to determine the frequency f₁ based on a single data frame is much smaller than the NT time of measuring that frame, which is a well-known advantage of IpDFT methods. Therefore, in this study, it is assumed that the computation time is negligibly small and that the current estimate f₁ appears in the output of the method with the last sample of a given frame, i.e., with the end of the measurement time of that frame. The next estimate f₁ is calculated after a sampling period T and it takes into account the last N samples of the signal. This updating of the estimate f₁ with each successive signal sample was adopted in the study in Section 5, in order to accurately demonstrate the effect of the phase φ₁ value, which changes with a small step for successive data frames. In practice, the calculation of the frequency f₁ can be performed less frequently because successive data frames differ by only one sample, i.e., they overlap to a large extent. Frequency calculation of f₁ with every kth sample reduces the computational power requirements of the DSP system by k times. For example, in the DSP system used in [31], the frequency f₁ is updated every second sample of the signal, i.e., every 2T ≈ 90.7 μs, where T = 1/f_s ≈ 45.4 μs (f_s = 22.05 kHz).

4.2. Assumptions for Zero Crossing Method

The zero-crossing (ZC) method was used as a reference method for the investigated IpDFT method from [31], since it is one of the most popular methods for estimating the frequency of a power network signal. This frequency is calculated in the ZC method based on the last two transitions of the signal through zero, thus by determining the duration of one half-period of the signal (Figure 13). It is possible to reduce the estimation error by taking into account many half-periods in one measurement instead of one, but the essential properties of the ZC method do not change significantly, and additionally the dynamic properties of the method deteriorate. Therefore, this paper is limited to the determination of a single signal half-period. In order to precisely determine the moment when the signal passes through zero, linear interpolation was used (Figure 13b), as is usually the case in the ZC method.

In the ZC method, low-pass filtering of the input signal is necessary, because without such filtering any high-frequency (impulse) disturbance introduces too large errors into the frequency estimation result. A digital filter of order 4–6 is most commonly used; however, results presented in Section 5.2 showed that higher-order filters are necessary to significantly reduce the errors of the zero-crossing method for disturbances due to interharmonics. In the studies presented in Section 5, three versions of a digital Butterwoth low-pass filter were adopted, each consisting of 1, 2, or 3 sections, respectively, of a 6th-order Butterworth filter with a cutoff frequency of f_c = 60 Hz (Figure 13a). The signal obtained as a result of such filtering is labeled in Section 5 with the labels ZC₁, ZC₂ and ZC₃ (Figure 14b and 15b), and this filtering introduces a signal delay that depends on the filter order (Figure 14b). The results of the frequency estimation of f₁ obtained using these three types of filtering are also labeled in Section 5 with the labels ZC₁, ZC₂ and ZC₃ (Figures 14f, 15d–f, 16e, 17e, 18d, 19e,g,i and 20b,c).

The ZC method is sensitive to the phase shift φ₀ of the first sample considered in the calculations with respect to the point where the signal passes through zero (Figure 13b). Therefore, in the studies in Section 5 for the zero-crossing method, the maximum value of the estimation error was determined from the set of values calculated for multiple values of φ₀.

5. Results

5.1. The Influence of Harmonics

The study of harmonics influence was carried out for a signal containing 12 consecutive odd harmonics with parameters defined in Table 2 and with THDS_U = 10.3% (Figure 9a). For the duration of the measurement window of 50 ms, the maximum estimation error depends on the phase φ₁ of the power network signal relative to the beginning of the time window, which results in a variation of the estimated value of f₁ on the time axis. Changing the measurement time by a few (or even by several percent) does not significantly affect the estimation accuracy (Figure 14c–e illustrates this effect for NT = 48.44 ms, 50 ms and 51.65 ms). In addition, the change in f₁ (as is the case in practice) at a constant NT value will not result in a significant change in accuracy, since the decisive factor in both of the above cases (i.e., changing NT at a constant f₁ and the change of f₁ at constant NT) is the value λ₁ = f₁NT = NT/T₁, and thus the ratio of the NT window duration to the signal period of T₁. The number of samples N does not affect the error values due to harmonics in the signal (Figure 14c–e). This property occurs when the condition of the sampling theorem is satisfied and when the dominant component of the estimation is the error due to harmonics in the signal. For the ZC method, a different situation occurs because the estimation error depends on the number of samples N. For N = 128. the estimation error is larger than for the IpDFT method, for N = 256 it is, depending on the phase φ₁, smaller or larger, and for N ≥ 512, the estimation error for the ZC method is almost always smaller than for the IpDFT method (Figure 14f). This ability to decrease the estimation error with increasing N for the ZC method is a beneficial property of the method, but it comes at the expense of the additional delay (Figure 14b) introduced by the low-pass filtering shown in Figure 13a. In the case of errors due to harmonics, the order of the filter used (6 for ZC₁, 12 for ZC₂ and 18 for ZC₃) does not affect the results of the ZC estimation (Figure 14f). For the IpDFT method and for the adopted parameters of this method, the estimation error does not exceed, depending on the quotient NT/T₁ the value, ca. 0.0003 Hz (Figure 14c), ca. 0.0002 Hz (Figure 14d) and ca. 0.0001 Hz (Figure 14e).

5.2. The Influence of Interharmonics

The study of the effect of harmonics was performed for a signal containing a group of interharmonics with parameters defined in Table 4 (Figure 15a). The influence of the interharmonic group on estimation accuracy of f₁ is larger than for the harmonics group. For the IpDFT method, the estimation error is less than ca. 0.0015 Hz and it does not depend, as in Section 5.1, on the number of samples N (Figure 15c). For the ZC method, in order to eliminate the influence of interharmonics, a 12th-order low-pass filter is necessary (ZC₂, Figure 15e) for N < 512, because the filter of order 6 causes a significant increase in estimation error (ZC₁, Figure 15d). For N ≥ 512, a full minimization of errors caused by interharmonics in the ZC method requires a filter order of 18 (ZC₃, Figure 15f), which is done at the expense of a higher signal delay (the signal with interharmonics after filtering is shown in Figure 15b, and the values of delays contributed by the filters ZC₁, ZC₂ and ZC₃ are provided in Figure 14b).

5.3. The Influence of Transient Overvoltages—High-Frequency Pulse

To study the effect of transient overvoltages of short duration, i.e., wide frequency band, the model (4) with oscillations of 1 MHz was used. Such a signal (Figure 16a) was added to the power network signal (Figure 16b) and the signal amplitude was limited to ± 400 V, due to the saturation effect of the input circuits and the A/D converter itself that occurs in practice. For N = 1024 (f_s = 20.48 kHz for NT = 50 ms), the maximum estimation error of f₁ is 0.1 Hz (Figure 16c). An increase of N by a factor of two to 2048 (i.e., increase by a factor of two of f_s) causes the reduction of the maximum error by ca. two times, i.e., to a value of ca. 0.05 Hz (Figure 16c). However, only adopting N = 2¹⁷, f_s > 2 × 1 MHz eliminates errors caused by failure to satisfy Shannon’s sampling theorem (Figure 16d). The comparison of the accuracy of the IpDFT method with the ZC method shows (Figure 16e) that the maximum errors are at a similar level, but the return to very low values of the estimation errors after the occurrence of the high-frequency pulse is much faster for the IpDFT method than for the ZC method. Taking into account the principle of the ZC method, it is obvious that it is not possible to use this method for this kind of disturbances without low-pass filtering at the input of the method. For the IpDFT method, however, the presented results indicate that it is advisable to apply an antialiasing filter at the input of the A/D converter in order to reduce the impact of such interference on the result of the estimation, but this is not a necessary condition—when it is not met, the maximum error is ca. 0.1 Hz (Figure 16c).

5.4. The Influence of Transient Overvoltages—Low-Frequency Pulse

To study the effect of transient overvoltages with long duration, i.e., narrow frequency band, model (5) with 1 kHz oscillations was used. Such a signal (Figure 17a) was added to the power network signal (Figure 17b). For N = 1024 and 2048, the conditions of Shannon’s theorem for sampling frequency are satisfied, and the maximum error of estimation of f₁ in the IpDFT method does not exceed 0.5 Hz and practically does not depend on N (Figure 17c). The return to very low estimation errors after a low-frequency pulse is much faster for the IpDFT method than for the ZC method (Figure 17d,e). The pulse modeled by (5) results in larger estimation errors than that modeled by (4), even though the condition of Shannon’s theorem on sampling frequency is satisfied, because the energy contributed by pulse (5) is much larger than that by pulse (4). Similar to the study in Section 5.3, it is clear from the performance of the ZC method that it is not possible to apply this method to low-frequency pulse disturbances without low-pass filtering at the input of the method. On the other hand, for the IpDFT method, the results presented show that the error rate does not exceed ca. 0.5 Hz with a rapid return to low values after the cessation of the disturbance.

5.5. The Influence of Frequency Changes

The impact of a step change in f₁ on the estimation results was examined for a signal for which the frequency changed in steps from 42.5 Hz every ca. 0.25–0.35 s (Figure 18a). Since the IpDFT method determines f₁ based on the last N samples of the signal, then after time NT the estimation error caused by the frequency step decreases to zero (Figure 18b,c). In this case, the IpDFT method also recovers much faster after a frequency step to very low values than the ZC method (Figure 18d).

5.6. The Composition of Harmonics, Pulse and Voltage Fall and Rise

To compare the influence of different types of disturbances, the estimation accuracy was determined for the signal in which harmonics, transient overvoltages and voltage dips occur, with the parameters defined in Table 5 (Figure 19a). In the IpDFT method, there are large changes in the estimation result (Figure 19d) due to the low-frequency pulse (at time instants 0.6 s and 1.2 s) and smaller changes due to the amplitude step change (at time instants 0.7, 0.9, 1.1 and 1.3 s). The determination of the type of disturbance only based on changes in the estimation result of f₁ does not seem simple, but in combination with the estimation of the signal amplitude (Figure 19f) there appears, as an effect of the application of the IpDFT method, an additional possibility of real-time detection and classification of the types of disturbances of the power network signal. Even the simplest estimation of the power network signal amplitude seems to be sufficient to perform such detection and classification. For example, Figure 19f shows the result obtained after determining the value of the complex spectrum for the determined f₁ (Figure 19d). For this purpose, one of the methods of determining one point of the spectrum for a given frequency, e.g., Goertzel’s algorithm, can be used. The plot of the estimation error on the logarithmic scale (Figure 19h) shows good dynamic properties of the IpDFT method, i.e., a fast return to very small estimation errors after the disturbance disappears.

The ZC method gives more than ten times larger estimation errors of f₁ than the IpDFT method for dip-type disturbances (Figure 19e). The ZC method also has significantly worse dynamic properties in this case (Figure 19g,i).

5.7. The Influence of Noise

The effect of Gaussian white noise on frequency estimation errors in the IpDFT method of [31] is described in more detail in [31]. For comparison with the ZC method, sample noise was added to a signal composed of several disturbances and shown in Section 5.6 and Figure 19a. The estimation results of such a signal with example noise are shown in Figure 20. According to the properties of the method described in [31] for a pure sinusoid with noise, the IpDFT method is unbiased, i.e., the error component due to noise is dominant. Moreover, the properties of the method in [31] show that the variance of the frequency estimator decreases with increasing number of N samples. The results in Figure 20 confirm these properties. At the moments of interference mentioned in Section 5.6, the estimation errors caused by these disturbances are larger than those caused by noise in the signal (Figure 20a). The properties of the IpDFT method for the noisy signal are also clearly more favorable than in the ZC method (Figure 20b,c) for both values of N (1024 and 2048) and regardless of the filter order in the ZC method.

6. Conclusions

The paper presents the use of the IpDFT spectrum interpolation method for estimating the fundamental frequency of power waveforms. The zero-crossing method with pre-filtering by three types of low-pass filters was used as a reference method. The waveforms created on the basis of recorded disturbances occurring in low-voltage power networks were used as test models of signals. The models included: voltage harmonics and interharmonics, oscillatory transient overvoltages exponentially damped with high and low frequency pulses, frequency spikes and disturbances consisting of harmonics, voltage dips, transient overvoltages and noise. The results of frequency estimation in the presence of modeled disturbances were analyzed as a function of the number of samples—sampling frequency, measurement window length and in the case of the ZC method for three types of filters.

Conducted research indicates that there is no significant influence of nonsynchronous sampling in the case of the IpDFT method, as it is in case of the ZC method. For a larger number of samples (above 512), the frequency estimation error for the ZC method is smaller than for IpDFT. A simplified summary of the obtained results are demonstrated in

Table 6. Maximum frequency estimation level errors for the grid disturbances described in Section 5. For two cases of N (1024 and 2048), NT = 50ms and for the presented IpDFT method with the use of GMSD window with parameter m = 3 and the reference zero-crossing method (ZC1 for one 6th-order filter section, ZC2 and ZC3—for two and three, respectively, filter sections).

Type of Disturbance	Type of Method	For Number of Samples N		Data Source
		1024	2048
The influence of harmonics (THDSu = 10.3%)	IpDFT	3.4 × 10−4 Hz		Figure 14f
	ZC1…3	2.0 × 10−6 Hz	1.8 × 10−7 Hz	Figure 14f
The influence of interharmonics	IpDFT	0.0016 Hz		Figure 15d–f
	ZC1	0.0052 Hz		Figure 15d
	ZC2	9.2 × 10−6 Hz	8.1 × 10−6 Hz	Figure 15e
	ZC3	1.4 × 10−6 Hz	2.4 × 10−7 Hz	Figure 15f
The influence of transient overvoltages—high-frequency pulse	IpDFT	0.10 Hz in pulse (above 0.001 Hz for 43 ms 1)	0.051 Hz in pulse (above 0.001 Hz for 41 ms 1)	Figure 16e
	ZC1	0.14 Hz in pulse (above 0.001 Hz for 60 ms 1)	0.068 Hz in pulse (above 0.001 Hz for 50 ms 1)	Figure 16e
	ZC2	0.13 Hz in pulse (above 0.001 Hz for 80 ms 1)	0.067 Hz in pulse (above 0.001 Hz for 70 ms 1)	Figure 16e
	ZC3	0.14 Hz in pulse (above 0.001 Hz for 90 ms 1)	0.071 Hz in pulse (above 0.001 Hz for 80 ms 1)	Figure 16e
The influence of transients overvoltages—low-frequency pulse	IpDFT	0.41 Hz in pulse (above 0.001 Hz for 45 ms 1)		Figure 17d,e
	ZC1	0.55 Hz in pulse (above 0.001 Hz for 60 ms 1)		Figure 17e
	ZC2	0.54 Hz in pulse (above 0.001 Hz for 90 ms 1)		Figure 17e
	ZC3	0.58 Hz in pulse (above 0.001 Hz for 100 ms 1)		Figure 17e
The influence of frequency changes	IpDFT	3.8 Hz in pulse (above 0.001 Hz for 45 ms 1)		Figure 18c,d
	ZC1	4.1 Hz in pulse (above 0.001 Hz for 100 ms 1)		Figure 18d
	ZC3	5.0 Hz in pulse (above 0.001 Hz for 150 ms 1)		Figure 18d
The influence of voltage fall and rise	IpDFT	0.13 Hz in pulse (above 0.001 Hz for 41 ms 1)		Figure 19g–i
	ZC1	2.3 Hz in pulse (above 0.001 Hz for 81 ms 1)		Figure 19g
	ZC3	3.8 Hz in pulse (above 0.001 Hz for 120 ms 1)		Figure 19i
The influence of noise	IpDFT	0.014 Hz	0.008 Hz	Figure 20a
	ZC1…3	0.020 Hz	0.012 Hz	Figure 20b,c

¹ The time in milliseconds of the temporary increase of the error above the level 0.001 Hz in the transient state caused by the disturbance. This time allows us to present the dynamic property of the compared methods.

.

The effect of harmonics (Section 5.1) in the ZC method is minimal, due to the use of input prefilters, while for the IpDFT method it does not exceed 0.34 mHz for the harmonic disturbance model used. The precision of fundamental frequency estimation is similar for both methods.

In the presence of interharmonic interference in the signal model (Section 5.2), the frequency estimation error increases from 0.34 mHz to a maximum of 1.6 mHz in the case of IpDFT. For the ZC method, better accuracies are obtained with higher filter orders (minimum 12th-order) and larger sample numbers (minimum 256). For example, for the filter order 6 (ZC1), the estimation precision is about 5.2 mHz, representing over three times higher frequency estimation error than for IpDFT.

In the presence of high-frequency damped oscillatory pulses (Section 5.3), the frequency estimation errors are at a comparable level when the disturbance occurs. In contrast, once the disturbance disappears, the IpDFT method reaches its initial accuracy much faster. The IpDFT method reaches its initial accuracy after only 43 ms (N = 1024), while the ZC method with the lowest-order filter reaches it after about 60 ms (N = 1024 ZC1) or 90 ms (N = 1024 ZC3). It is also noticeable that there is a characteristic short oscillation in the frequency estimated waveform for IpDFT immediately after the disturbance, which can be used to identify the disturbance.

Low-frequency pulses (Section 5.4) have a greater impact on the fundamental frequency estimation than do high-frequency pulses. In the former case, the estimation error is 0.1 Hz (N = 1024), compared to 0.41 Hz (N = 1024) in the latter case. The value of the error depends primarily on the energy of the disturbance. Comparing with the ZC method, one observes, similarly to the high-frequency pulse, a significantly increased recovery time to the original accuracy after the occurrence of the pulse. The situation is very similar when stepping the frequency (Section 5.5). The time to return to the initial accuracy is at least three times faster for the IpDFT method than for the ZC method.

Finally, the IpDFT method has better dynamic performance under all types of disturbances (Section 5.6) with good estimation of the fundamental frequency. Moreover, based on the analysis of frequency estimation error variations and instantaneous amplitude variations of the fundamental frequency component, it is possible to determine the type of disturbances: pulse type, rms variations and frequency variations. Tests in the presence of noise (Section 5.7) also prove better dynamic properties of the IpDFT method in comparison to ZC, and the frequency estimation accuracy is higher for the first method.

In conclusion, the IpDFT method may find application in protection systems and measuring devices monitoring rapid changes of the fundamental frequency, even in the presence of significant disturbances of the power waveform voltage.