Power System Dynamic Frequency Measurement Based on Novel Interpolated STFT Algorithm

As the product of time and frequency resolution of the Short-Time Fourier Transform (STFT) constant, the dynamic frequency measurement accuracy provided by the STFT is reduced under asynchronous sampling. This paper proposes a novel Interpolated Short-Time Fourier Transform (IpSTFT) based on the Blackman window. The Blackman window is adopted to reduce the spectral leakage, and the spectral interpolation procedure is applied to eliminate the picket-fence effects. Results of simulations are provided to show that the proposed method can improve the accuracy of dynamic frequency measurement significantly with low time consumptions.


Introduction
Frequency is an important parameter of power system protection, power quality monitoring, operation, and control of devices using digital technologies.The ideal power system sinusoidal waveform is pure, continuous, and of a constant fundamental frequency.However, due to the mismatch between power generation and load demand, frequency in practical power system varies over a small range which poses a threat to efficiency and safety of entire system [1].Thus the accurate measurement and tracking of system frequency is of utmost importance [2] and [3].
Power system frequency could be estimated commonly by zero-crossings detection, phase-locked loop, Kalman filter, least-squares methods, Prony method, Matrix Pencil, et al.However, the above mentioned methods suffer from poor dynamic performance or high computational burden for dynamic frequency measurement [4].Conventionally, the Discrete Fourier Transform (DFT) and its accompanying Fast Fourier Transform (FFT) are efficient for power system frequency measurement because of the simplicity and easy implementation.DFT-based methods generally require an accurate match between the periodicity of DFT data set and the periodicity of signal waveforms, i.e. the synchronous sampling, otherwise the spectral leakage and picket-fence effects will occur [5].The windowed interpolation FFT algorithms are proposed to reduce spectrum leakage and picket-fence effects caused by asynchronous sampling, which can significantly improve the accuracy of frequency measurement with a steady-state assumption [6].Although the DFTbased methods perform well for the infinite time case of a stationary signal, they cannot be used to resolve any temporal information associated with these fluctuations [7].Since the power system frequency is subject to small random deviations, it is required to depict variation of the signal's spectrum with time [8].Thus a series of time-frequency algorithms, including STFT, Wavelet Transform (WT), S Transform (ST), and Hilbert-Huang Transform (HHT) are applied for power system dynamic frequency measurement.However, those approaches may either suffer from low solution accuracy, cross-terms or less computational efficiency.For example, the STFT has the effect of giving some time resolution to the measurement at the expense of frequency resolution; the WT requires expensive computation in the decomposition of waveforms; the ST suffers from poor energy concentration; HHT has the weakness of cross-terms.For all above algorithms, the STFT has the lowest computational burden due to its realization through FFT.That is why the improved STFT-based approaches are called for as an important research field even now.Recently, many improved algorithms based on the STFT have been put forward.A commonly used method to improve the time resolution of STFT is the so-called variable window procedure.For instance, in [7], a variable window STFT to analyse voltage and current signals for transmission line protection is proposed.In [8], an algorithm inspired by the intersection-of-the-confidence intervals rule is proposed for obtaining suboptimal window width in the STFT.In [9], the threshold for autocorrecting the noise affected frequency estimates is chosen according to the frequency resolution determined by the STFT window size.Moreover, many studies have been conducted to improve the time-frequency resolution of STFT by combing additional operations.For example, in [10], the synchronized STFT is proposed for time-frequency analysis by using a synchronized linear transform.In [11], a novel method to reduce distortions is presented by using the STFT and the wavelet packet filter banks.In [12], a time-stretched STFT is proposed and demonstrated to overcome the limitation in the time-frequency analysis of ultrafast signals.In [13], a modified STFT is presented to estimate signals from short-time magnitude spectra.However, few investigations have been done on the application of grid dynamic frequency measurement by STFT because of the low time-frequency resolution and low accuracy.
As mentioned above, the crucial drawback of the STFT is that the product of time and frequency resolution is constant.Especially, the frequency measurement accuracy provided by STFT is reduced when the sampling frequency is not synchronous with the analogue signal's frequency.To improve the accuracy of dynamic frequency measurement provided by STFT, this paper proposes the Interpolated Short-Time Fourier Transform (IpSTFT) based on the Blackman window.In the proposed method, the Blackman window is adopted to reduce the error caused by the spectral leakage, and the spectral interpolation procedure is applied to eliminate the error caused by the picket-fence effects.Results of simulations and practical applications are provided to show that the proposed method can improve the accuracy of dynamic frequency measurement significantly.
The organization of this paper is as follows.The principles of STFT are reviewed in Sec. 2.
The concept of IpSTFT and the calculation formulas are presented in Sec. 3. Simulation result is given in Sec. 4. Finally, conclusion is drawn in Sec. 5.

Short-Time Fourier Transform
In order to solve the shortcomings of Fourier transform, Gabor proposed the STFT in 1946, which can be used for time-frequency analysis of non-stationary signals where the FT is not adaptable.The STFT which is an extension of FT provides an insight in the time-evolution of each signal components by decomposing the time-varying signal into time-frequency domain components [14] and [15].And the continuous STFT is expressed as where x(t) is the continuous-time signal, w(t) is the window employed.Its corresponding discrete form is described as where x(k) is the discrete signal, w(n) is the discrete form of window employed, N is the total number of samples.When the w(k − n) is sliding along the x(k) followed k varies, the discrete signal is truncated and localized into a plurality of locally stationary signal, and then the windowed small signal is analysed utilizing FFT [16] and [17].
The result of STFT is a two-dimension complex matrix, and the magnitude of STFT can be expressed as The column of P (m, n) corresponds to the sampling time, each column vector is the spectrum of the corresponding moment.Similarly, row corresponds to the frequency bin and each row vector represents spectrum of one fixed frequency bin in time distribution.Matrix elements are the spectral amplitude.According to magnitude P (m, n), the time-frequency characteristic can be acquired.
It is known that the use of STFT will result in a trade-off between frequency resolution and time resolution, which is determined by the adopted window w(t) from Eq. (1) [18].Obviously, with different windows as well as flexibility of different window length, STFT can be used to analyse non-stationary signals with different time-frequency resolution [19].Notice, for a given frame (or window) length, the trade-off between frequency resolution and time resolution can be optimized by applying a rectangular window to each frame, which will provide the best frequency resolution at a cost of higher side lobes compared to other windows.Considering that the desirable sidelobe performances of adopted windows, i.e., a low-peak sidelobe level and a high-sidelobe roll off rate, are indicative of a small spectral leakage, which are very important for suppressing harmonic interference in frequency measurement [20] and [21].Thus, it is of interest to improve the frequency resolution without loss the desirable sidelobe performance.

The Proposed Approach
Assume that the power system frequency is dynamic and it is presented as where A 0 , f , ϕ 0 are the amplitude, frequency and initial phase of the signal, respectively, while the f is a time variable.
The discrete form obtained from x(t) with the sampling frequency f s is described as where k is the discrete time index (0 ≤ k < N ).
In this paper, Blackman window which is a three order raise cosine window is adopted in the STFT, and it can be expressed as where L is the window length and i = 0, 1, . . ., L − 1.
The DFT of the Blackman window can be represented as where W R (λ) is the DFT of the rectangle window and λ = 0, 1, . . ., L − 1.
The Blackman-window-based STFT can thus be obtained by substituting Eq. ( 6) into Eq.( 2), which is (8) Equation ( 8) can be rewritten in the form of matrix, as follows: The α m which named frequency vector is described as and the β n which is called window vector can be given by where m = 1, 2, . . ., (L/2) + 1, n = 1, 2, . . ., fix(N − D)/(L − D), in which fix() means the rounding function and D is the duration between two times of FFT when calculating the STFT.H k m corresponds to the transformation kernel that is defined as By ignoring the effect of negative frequencies, the value of α m β n can be written as where W E represents the spectral of window employed, k 0 (n) is the exact spectral line corresponds to signal frequency f .
In the case of synchronous sampling, the effective position of signal frequency f is sitting right at the frequency bin of column n, i.e., k 0 (n) = f L/f s is an integer, where the frequency can be calculated by using the corresponding frequency bins.However, it is almost impossible to achieve synchronous sampling.With non-synchronous sampling, k 0 (n) can be written in two parts, i.e., k 0 (n) = θ(n) + ξ(n), where θ(n) and ξ(n) are the integer and fractional parts, respectively.The integer part θ(n) can be determined by applying a simple peak search procedure to the column n of P B−ST F T (m, n) which is described as follows Assume that the | α m1 β n | represents the local largest magnitude in P B−ST F T (m, n) of column n, therefore, we can obtain that θ(n) = m 1 (n).Notice that the direct utilization of the integer part θ(n) for frequency measurement with non-synchronous sampling would produce undesired errors.
In this paper, the spectral interpolation procedure is applied on the matrix P B−ST F T (m, n) which could improve the accuracy of dynamic frequency measurement.In the proposed interpolation procedure, the fraction part ξ(n) is calculated via using spectral lines of the local largest and second largest magnitude in of column n, respectively.Thus the interpolation coefficient ψ(n) can be demonstrated as Substituting α m1 , α m2 and β n into Eq.( 15), the coefficient ψ(n) can be rewritten in Eq. (16).
where the symbol '±' must be the same at one time.
So the result of coefficient ξ(n) can be calculated by the inverse of Eq. ( 16) which could be represented as where the ψ(n) is independent variable and the ξ(n) is the dependent one.
In order to reduce the error, we apply the polynomial approximation method to calculate the ξ(n).The polynomial approximation formula of the coefficient ξ(n) can be described as: The accuracy is determined by the order i of the polynomial, and for usual, when the order i is 7, the polynomial is effective and sufficient for most implementation, which can be written compactly in matrix notation as where C is the coefficient vector of window employed that can be described as and the vector Q is demonstrated as The result of C is estimated by least-squaresmethod-based equation which can be represented as where e is the error between the test data and results calculated by Eq. ( 17).We can obtain the vector C of different windows by solving Eq. ( 22).
The coefficient vector C B of Blackman window is [1.96043163, 0.15277325, 0.0742583, 0.04998548].So that we can obtain the polynomial of coefficient ξ(n) based on the Blackman window as follows: Finally, the Blackman-window-based dynamic frequency f d (n) of the signal is obtained as From Eq. ( 14), it is known that if the frequency fluctuates in a small range around its normal value which is 50 Hz in China, as well as the window length L is determined, we can lower the time consumption significantly by calculating the values of 2 or 3 rows only (usually m 1 , m 2 , and the third largest if necessary), i.e., the amount of calculation will reduce to 4/L or 6/L compared with the original.The reason is that the spectral line, which corresponds to the band where frequency varies, may just lie between two of the largest three spectral lines.
The workflow of the proposed IpSTFT is shown in Fig. 1.

Simulations
In order to evaluate the time-frequency performance of the proposed method, two kinds of signals are implemented in this section.One is frequency-hopping signal while the other is chirp.To make comparisons, the Blackman-window-based FFT, Morlet-based WT, Blackman-window-based STFT and ST are adopted.The parameter D is set to 24 while the length L of Blackman window is 256.

Frequency Measurement of Frequency-Hopping Signal
The simulated frequency-hopping signal can be written as where the range of fundamental frequency f is from 49.5 Hz to 50.5 Hz with varying in steps of 0.1 Hz.The variations of f and the hopping time range are tabulated in Tab. 1.
Tab. 1: Frequency variation of simulated signal.The sampling frequency f s of the simulated signal is 2000 Hz and the total length is 8 • 10 3 .The time domain waveform of the simulated signal is shown in Fig. 2.

1)
Frequency Measurement Based on FFT The result of FFT with Blackman window is a single value 50.2139Hz which is meaningless and cannot represent the variation of the simulated signal.That means the FFT is not suitable for dynamic frequency measurement.

2) Frequency Measurement Based on WT
Wavelet transform in this section uses the Morlet Continuous Wavelet Transform whose scale are 256 and 512, respectively.The peak frequency (frequency bin which the maximum value in each column corresponding to from the matrix of time-frequency analysis) characteristics are plotted in Fig. 3.In Fig. 3, the blue and green lines represent the results of scale 512 and 256 WT respectively.It can be seen from the graph apparently that the 256-scale WT is constrained whose result is a constant over the entire time range.Similarly, the 512-scale WT cannot give a specification about variations even though the result is better.From Fig. 3, we know that the WT hardly meets the requirement of dynamic frequency measurement on account of the low time-frequency resolution.

3) Frequency Measurement Based on ST
The sampling frequency interval of ST applied in this section is 1 and 2. Figure 4 demonstrates the peak frequency calculated.
From the result of ST in Fig. 4, it is known that the variations are represented obviously as well as the changing moment.Meanwhile, we see that the result is better when the sampling interval is 1.However, at 50.1 Hz, 50.3 Hz, 49.9 Hz, 50.4 Hz and 49.6 Hz, there are errors between the estimations and sets.This is to say the frequency characteristic is still described incorrectly by the curves due to the inadequate timefrequency resolution.
Accordingly, we can say that the performance of ST is better than WT.Unfortunately, the variations still cannot be measured clearly and correctly to fulfil the requirement of power dynamic frequency measurement.

4) Frequency Measurement Based on STFT
The STFT used here is based on Blackman window whose window length is 1024 and shift length is 24.The simulation result is shown in Fig. 5.
In Fig. 5, the blue line on the background of calculated values being displayed by kinds of colours.In the figure, the lighter colour corresponds to a larger value of the time-frequency matrix while the darker colour corresponds to a smaller.From Fig. 5, it can be concluded that the by using STFT, it is difficult to obtain the specification of dynamic frequency.Furthermore, there is no data in the range 0∼0.26 s and 3.74∼4 s due to the window adopted as well as edge effects.Analysis above shows that due to the poor resolution, STFT cannot meet the requirement either.

5) Frequency Measurement Based on IpSTFT
The window length L and shift length D of Blackmanwindow-based IpSTFT are 256 and 24, respectively.The STFT (window length 256) here is performed as a comparison with the result illustrated in Fig. 6.The response of the proposed algorithm is plotted in Fig. 7. Figure 6 shows that the STFT cannot give a satisfied result while the window length is 256 only.However, when it comes to IpSTFT, the result changes a lot.The time-frequency characteristic given by the curve in Fig. 7 shows that over the range from 49.5 Hz to 50.5 Hz, the specification of variations can be extracted accurately.
The relative errors of frequency measured by different algorithms are provided in Tab. 2. The results with shadows indicate the minimal relative errors among the adopted algorithms.The number represents the scale adopted in WT(256), WT(512).As to ST(2) and ST(1), the number means the sampling frequency in-  From Tab. 2, we know that the relative error of IpSTFT is smaller than that of others.In band of 49.5∼50.5Hz, the relative error of IpSTFT is 2∼3 order of magnitude lower than that of WT and STFT.Compared with ST, the error is 1∼2 order of magnitude lower.Consequently, we could say that the relative error of IpSTFT is the least among the algorithms applied.From the analysis, it is known that the result of frequency measurement is accurately determined when using the proposed algorithm, as well as the stability  and reliability both meet the requirement of dynamic frequency measurement.

Frequency Measurement of Chirp Signal
The      In Fig. 8, the STFT based on Blackman with window length 256 and 1024 are depicted.It can be seen that the concentration of STFT(1024) is better than that of STFT(256), which means the time-frequency performance of STFT(1024) is better.However, both STFT(1024) and STFT(256) cannot demonstrate the accurate frequency on the time axis.Next, Fig. 9 gives the outcomes of wavelet with different scale values, but the frequency variations results are still blurs.Figure 10 presents the results of S transform which gives a perfect response to the simulated signal.However, the end effects and the worse concentration at the higher frequency, especially the huge computational burden make it difficult to implement in embedded power system dynamic frequency measurement.
As seen from Fig. 8, Fig. 9 and Fig. 10, for all the techniques that are considered, the calculated values of frequency cannot reach a precise agreement with the presented scenario.However, utilization of the proposed algorithm can solve the issue perfectly, and the simulation result is given in Fig. 11.In Fig. 11, the frequency of the given signal is calculated by IpSTFT, and the outcome is tremendously keeping pace with the variations of simulated signal which shows the accuracy of the proposed algorithm.

Analysis with White Noise
Measurement of dynamic frequency signal corrupted by white noise is another simulation carried out to evaluate the proposed algorithm.The signal applied here is demonstrated in Eq. ( 24  IpSTFT based on different windows, the simulation of sinusoid signal whose frequency f = 50.1 Hz without any disturbance is adopted and the result of relative error else with the absolute error are tabulated in Tab. 3. Figure 12 shows that the variances are inversely proportional to the SNR tolerably, which means that more accurate results can be achieved when there is less noise in the signal channel.That means that the accuracies of frequency measurement decrease, if the signal is corrupted by white noise.In addition, we can see that the Blackman window provides the minimum variance among the classical windows adopted, as well as the considerable robustness.And, when the improved windows are employed, i.e., 4-order TSCW (Triangle Self-Convolution Window) and 4-term MDW (Maximum sidelobe Decay rate Window), the performance of Ip-STFT would turn better on accuracy significantly.The variances of 4-term MDW are inversely proportional to the SNR totally, which represents high precision while signal suffers from weak white noise (SNR>60).
Table 3 shows the accuracy comparison of pure sinewave signal based on IpSTFT with different windows including BW(Blackman Window), HNW(Hanning Window), HMW(Hamming Window), TW(Triangle Window), 4-order TSCW, and 4-term MDW.From Tab. 3, it is obvious that the Blackman window offers the best performance among the classical windows.And if the higher accuracy is required, IpSTFT could also be satisfied by employing new windows, i.e., 4order TSCW and 4-term MDW.The simulations in this part show that IpSTFT can provide both suitable measurement bias and variances, which means a better performance in application of power system dynamic frequency analysis.

Simulation of Time Consumption
In this paper, the simulation platform consists of Matlab R2012b and Lenovo Y470 with the processor core I5-2450M@2.5 GHz.The average elapsed time of 500 times with different algorithms is shown in Tab. 4.
As can be seen from Tab. 4, the time consumption of IpSTFT else with STFT is much less than that of WT and ST.However, the resolution of STFT is fixed and the accuracy of frequency measurement provided is inadequate.Notice that the average running time of IpSTFT(256) is almost tenfold lower than that of WT and ST.The large calculation amount of WT and ST prevents them being adopted in embedded systems, because it cannot provide the basic resource what the algorithms need.Through the simulation, it is known that IpSTFT is up to power system dynamic frequency measurement due to the accuracy of frequency measurement and lower time consumption.Besides, the advantage of lower time consumption also makes it possible of being used in embedded systems for power system applications.

Conclusion
In this paper, a novel interpolated STFT algorithm for power system dynamic frequency measurement is proposed.The limitations of the traditional timefrequency analysis, which inherently require huge computational burden and the time-frequency resolution deficiency, have been resolved by the proposed algorithm.Utilizing the frequency domain interpolation procedure offers a convenient way to achieve high accuracy of dynamic frequency measurement.Simulations based on hopping signal and chirp signal are performed to indicate the effectiveness and accuracy of the proposed method in comparison with STFT, ST, and WT.In addition, the simulation under white noise is carried out to evaluate the robustness of the proposed algorithm.As a result, it can be seen that the proposed IpSTFT is suitable for measuring localized transient behaviour of power system dynamic frequency in embedded systems.Furthermore, the proposed algorithm in this paper can be developed not only to disturbance evolution assessment, but also for the classification and identification of the power quality disturbance.
Time frequency analysis of chirp signal based on STFT with constant length L = 256.Time frequency analysis of chirp signal based on STFT with constant length L = 1024.
chirp signal of sampling frequency 2000 Hz in the simulated time 0∼8 s its component increasing from 20 Hz to 90 Hz.The total number of samples is 1.6 • 10 4 , and the results of the techniques mentioned above are delineated in Fig. 8, Fig. 9, Fig. 10 and Fig. 11 .
Time frequency analysis of chirp signal based on WT with scale L = 256.
Time frequency analysis of chirp signal based on WT with scale L = 512.
Fig. 12: Variance of frequency measurement obtained by simulation with different SNR value when using different windows.
Frequency variation of simulated signal.
12: Variance of frequency measurement obtained by simulation with different SNR value when using different windows.c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING Tab.3: