The Periodic Characteristics of Harmonic Measurement Errors with the Initial Sampling Time

2.1 Power system harmonic algorithm based on windowed DFT

The signal model for harmonic analysis can be expressed by

where Am and φm are individually the amplitude and the initial phase of the mth component, ω0=2πf0 is the fundamental frequency. With equally sampling space Ts for nearly p periods (satisfying the Nyquist Sampling Theorem), where p is determined by the mainlobe of the adopted window, the signal is then digitized with length N

Then the windowed sequence is

in which w(n) is the sampled window function.

The DFT of the xW(n) can be given by [6]

where λ is the synchronous error, W(ω) is the frequency response of the window, Bm=Amejφm/2j, and Bm∗ stands for the complex conjugate of Bm.

From (4), it can be known that all components in the signal make their contribution to the DFT’s of the mth component because of the leakage caused by the limited window width and unsynchronized sampling. With the leakage coming from the negative part of the mth harmonic and other components neglected [6],

(5) can be used to estimate the amplitude of the mth harmonic, where G=∑n=1N−1w(n). It immediately follows that if λ=0, we have Aˉm=Am. As the synchronous error always exists practically and λ is always a nonzero value, the estimated value Aˉm≠Am. Windowed algorithm is the most popular way to improve the algorithm accuracy, the stability and consistency for each measurement with the same signal should also become a concern.

2.2 The sensibility of the algorithm precision with the initial-sampling time

For the sake of convenience, use the periodical exponential signal model

where the fundamental frequency is 50 Hz.

Assume the synchronous error is λ for one fundamental period, Hanning window is adopted, and the signal is only consisted of the fundamental component and the 3rd harmonic. Consider the amplitude of the fundamental component, which can be expressed by

where X1=A1ejφ1W(−2λ)/G,X3=A3ejφ3W(−4−6λ)/G and G=∑n=1N−1w(n). Note that Aˉ1 is the magnitude of the additive result of vector X1、X3.

If λ=0, we have Aˉ1=A1.

Then if the synchronous error λ≠0 is fixed, we define W(−2λ)=a1ejb1, W(−4−6λ)=a3ejb3, and notice that the magnitudes of X1 and X3 remain constant while the time t changes。

The initial-sampling time has strong relationship with the initial-sampling phase, at the moment t=t0=0, the initial-sampling phase of the fundamental component is φ1 and angle of vector X1 is φ1+b1, so the 3rd’s is φ3 and X3’s is φ3+b3. Supposing that the angle difference of the two vectors X1 and X3 is φt0=φ1+b1−φ3+b3≠0, we can find that there is always at a t=tm that the initial-sampling phase satisfied 2π×50tm+φ1+b1=2π×150tm+φ3+b3, because of the randomicity of the initial-sampling time. It can be easily calculated that the angle difference φtm=0 at the moment tm=φt0/200π. The result is shown in Figure 1 for example, where the dashed lines stand for X1=220×ejπ/2,X3=35×ej0, and their additive results (X₁+X₃) at time t₀, and the solid line is the vector of (X1+X3)attm=0.01. It means that the calculated amplitude of the fundamental component is different at time t₀ and t_m, and the same result can be found for the 3rd component.

Figure 1:

The amplitude of the fundamental component at t₀ and t_m.

For the Mth signal model and still only the amplitude of the fundamental component is considered, from (4) it immediately follows that

where X˜2=∑k=2MXk,Xk=AkejφkW(1−k2−2kλ)/G.

It can also be found that the magnitudes of vector X1,…,XM never change, while their angles always change with the initial-sampling time t, if the synchronous error λ≠0 is fixed. Note that if the instant is changed from t0 to t0+Δt, the angle variety of vector Xk is Δφtk=k×2π×50×Δt. Especially if Δt=l×0.02 and l is integer, the angles of vector X1t0,…,XMt0 at instant t0 are equal to the angles of vector X1tm,…,XMtm at instant tm=t0+Δt respectively. It means that X1,…,XM are period vector changed with the variable t, and the period is the same as the signal (1). When interharmonic components are concerned in the case, the angle variety of interharmonic vector, with frequency pik/qik×f0,isΔφtik=pik/qik×2π×50×Δt, and p_ik., q_ik are integers. The same result can be found, and the period is ∏qik×0.02 for the time, in which ∏q_ik means the lease common multiple of denominator q_ik. Thus it can be concluded that the estimated amplitude Aˉm of the mth component by windowed DFT are changed periodically with the initial-sampling phase or the initial-sampling time.

2.3 Simulation results

In this section, we present a simulation example to show the sensibility of the amplitude precision with the initial-sampling phase. Consider the special case of the signal model for simplification, where φm=mφ0, M=35,Am=1/m and the fundamental frequency is 50 Hz. The equally sampling space is Ts=1.01×2π/128,which means that the synchronous error for one signal period is 0.01×2π. The relative error of the mth component’s estimated amplitude is defined as

The results with Hanning window, rectangular window, and Blackman window are shown in Figures 2, 3 and 4, for φ0=0∘,45∘,90∘ respectively. From figures, we can find that the initial-sampling phase does have influence on the estimation precision. The same results can be found with other windows.

Figure 2:

The relative amplitude errors of harmonics with Hanning window.

Figure 3:

The relative amplitude errors of harmonics with rectangular window.

Figure 4:

The relative amplitude errors of harmonics with Blakman window.

Figure 5 presents the normalized magnitude frequency response of these three windows. From the figure, it can be found that rectangular window has the poorest magnitude frequency response, while Blackman window has the most satisfactory characteristic, which means that Blackman window can most effectively suppress the spectral leakage and lead to more accuracy of harmonic estimation with windowed DFT than rectangular window and Hanning window. Figure 6 shows the relative amplitude errors of harmonics with these three windows for φ0=90∘. To our complete surprise, the results of harmonic amplitude estimation with Blackman window by DFT suffer much more error than rectangular window and Hanning window, and the results with rectangular window are most accurate. That is to say, the window which has better characteristics may give relatively worse result at some phases.

Figure 5:

Normalized magnitude frequency response of rectangular window, Hanning window, and Blackman window (a) main-lobes, (b) side-lobes.

Figure 6:

The relative amplitude errors of harmonics with rectangular window, Hanning window, and Blackman window for φ₀=90°.

Figure 8(a) shows the results of the rectangular window when φ0 changed within [0~360⁰], which is the further proof of the same phenomenon.

To verify the estimated amplitude of the mth component by windowed DFT are changed periodically with initial-sampling phase or the initial-sampling time. We change the initial-sampling phase φ₀ from 0° to 720° with an increment of 1°, and only consider the relative amplitude errors of 3rd harmonic for simplification. Test results with rectangular window are shown in Figure 9(a) which changed periodically as expected.

It is very well known in power systems that the amplitude of the harmonic components get attenuated significantly with the increase in harmonic order, the effect on harmonic signals beyond the 15th harmonic order are also investigated. At this time, the synchronous error for one signal period is changed to 0.03%×2π according to IEC standard, and the rectangular window and Hanning window are considered. Figures 10(a) and 11(a) shows the results of the 35th harmonic relative amplitude with the rectangular window and Hanning window individually, when φ0 changed within [0~720⁰] with 1° step. From the figures, the periodic characteristics of harmonic measurement errors with the initial sampling time can be clearly found. As the synchronous error is decreased, the relative amplitude errors in Figure 10(a) are much smaller than Figure 8(a).

2.4 IEC test results

In part B and C, we studied the sensibility of the algorithm precision with initial-sampling phase by theoretical analysis and matlab simulation. And the aim of this part is to discuss the relationship between the IEC technique accuracy with the initial-sampling time, and the signal xt=11.5sin(3×2πf0×t)+11.5sin(3×2πf0×t) is considered in the case with f₀=50.05 Hz, which consists of the 3rd harmonic and the 5th harmonic. The test system, shown in Figure 7, is mainly composed of 4 components: (1) signal generator (Electrical Power Standards Fluke 6100A), which is qualified for IEC signal generating; (2) a carefully designed 16bit AD sampling board [23], which can be used to obtain the samples from Fluke 6100A; (3) notebook (Thinkpad X220i) is used to calculate the harmonic components from the sampling data; (4) oscilloscope(TDS3000) is for the use of signal testing. Test results for the initial-sampling time changed in [0–0.02 s] with step 0.0001 s are shown in Figure 9(a), with a fixed sampling frequency 10 kHz. It can be observed that estimation errors are changed periodically, and the maximum relative error of group3 is nearly ±0.3% in the case.

Figure 7:

Test system.

3.1 Moving-average windowed DFT algorithm

It is shown in Section 2 that, when the window and the synchronous error are fixed, the amplitude precision of each frequency component is changed periodically with the initial-sampling time t. Since no prior information about the time t is available, the uniform distribution defined within an interval [−T/2∼T/2]T→∞ can be adopted. When the probability density of φmti, the subscript mti denoting the initial phase of the m th component at a time ti, is concerned, it can also be considered as a random variable with a uniform distribution changed within the interval [0∼2π] because of its strong relationship with the moment ti. Then we can define Aˉm(φmti) and εm(φmti) which stands for a random variable changed with φmti. To eliminate the influence on harmonic precision caused by the random variable φmti, the expected value E[Aˉm(φmti)] can be considered as the true amplitude estimation of the mth component. Note that this algorithm is similar as the well known “shifting window average DFT”, but here we adopt this traditional method to mitigate the random influence and to achieve more stable measurement precision. Also this method can only work with the signal remaining constant for nearly p+1 periods, where p is determined by the mainlobe of the adopted window.

The steps of the algorithm are as follows:

1. Obtain x(n) with length N for nearly p periods, and each period has Ns samples, then extend the sampling interval to p+1 period, thus the sampling sequence is x(n),n=0,1,…,N,…,N+Ns−1deg

2. Compute the amplitude of the mth component as Aˉm(φmt0),…,Aˉm(φmtNs−1) by windowed DFT with length N at instant n=0,…,Ns−1 individually

   (10)A¯m(φmt0)=2|∑n=0N−1x(n)w(n)exp(−j2πNmn)| /G A¯m(φmt1)=2|∑n=0N−1x(n+1)w(n)exp(−j2πNmn)| /G                             ⋮ A¯m(φmtNs−1)=2|∑n=0N−1x(n+Ns−1)w(n)exp(−j2πNmn)|/G. 

3. The expectation is

   (11)E[Aˉm(φmt)]=∫02πAˉm(φmt)×12πdφmt≈∑ti=0Ns−1Aˉm(φmti)×12π×2πNs=1Ns∑ti=0Ns−1Aˉm(φmti)

4. The calculation of the relative amplitude error of the m th component is

   (12)εWm(π,φ)=E[Aˉm(φmt)]−1/m1/m

3.2 Recursive algorithm

Direct calculation of (10) strongly increases computation burden, even the FFT can be used. Further reduction is necessary and also possible, if rectangular window or some simple cosine windows are considered. The DFT of x(n) with rectangular window can be expressed as

If X_t0(m) is obtained with FFT, X_t1(m)... can be calculated by the following simple recursive algorithm

and then Ā_m (φ_mt0), … , Ā_m (φ_mtNs-1) can be obtained individually. The computation requirements are considerably lower using (14) rather than (10). If cosine windows are adapted, for example Hanning window, X_t0(m), X_t1(m)... can be first calculated with the steps mentioned before, then the digital spectrum with Hanning window can be obtained by using (15), and finally the expectation of the amplitude can be obtained [3].

3.3 Test results

Under the same sampling condition, test results are shown in Figures 8(b), 9(b), 10(b), 11(b), 12(b) respectively. From figures, we can observed that the proposed algorithm alleviate the initial-sampling instant influence on the harmonic precision, and calculated results are more stable and consistent as expected, when the window and the synchronous error are fixed.

Figure 8:

The relative amplitude errors of harmonics with rectangular window when φ_° is changed within [0~360_°]. (a) traditional method, (b) moving-average windowed DFT method.

Figure 9:

The relative amplitude errors of the 3rd harmonic with rectangular window for different phases. (a) traditional method, (b) moving-average windowed DFT method.

Figure 10:

The relative amplitude errors of the 35th harmonic with rectangular window for different phases. (a) traditional method, (b) moving-average windowed DFT method.

Figure 11:

The relative amplitude errors of the 35th harmonic with Hanning window for different phases. (a) traditional method, (b) moving-average windowed DFT method.

Figure 12:

The relative amplitude errors of the 3rd harmonic and the 5th harmonic. (a) traditional method, (b) moving-average windowed DFT method.