Frequency estimation of the weighted real tones or resolved multiple tones by iterative interpolation DFT algorithm

Digital Signal Processing 41 (2015) 118–129

Contents lists available at ScienceDirect

Digital Signal Processing www.elsevier.com/locate/dsp

Frequency estimation of the weighted real tones or resolved multiple tones by iterative interpolation DFT algorithm Jiufei Luo ∗ , Zhijiang Xie, Ming Xie Department of Mechanical Engineering, Chongqing University, Chongqing 400044, People’s Republic of China

a r t i c l e

i n f o

Article history: Available online 10 March 2015 Keywords: Parameter estimation Discrete Fourier transform Spectrum analysis Weighting function

a b s t r a c t Accurate frequency estimation of sinusoidal signal is commonly required in a large number of engineering practice scenarios. In this paper, an iterative frequency estimation algorithm is proposed based on the interpolation of Fourier coefficients of weighted samples. This algorithm is nearly applicable for all conventional window functions. Systematic errors for various windows are presented and the performance of the proposed algorithm is investigated in the presence of white Gaussian noise. The simulation results demonstrate that errors caused by a mistaken location of the spectral line can be significantly reduced. The proposed algorithm is straightforward to implement, has high precision, good compatibility, and is robust against additive noise, all of which make it ideal for accurate frequency estimation in spectral analysis. © 2015 Elsevier Inc. All rights reserved.

1. Introduction A fundamental problem encountered in many fields of signal processing (e.g., communication systems, biomedical signal processing, vocal analysis, vibration analysis, power system analysis) is the accurate estimation of frequency, amplitude, and phase for a mono-tone or multi-tone signal corrupted by random noise. Many algorithms have been presented in the literature and can generally be divided into two categories: time-domain and frequencydomain [1]. Frequency-domain algorithms are based on the windowed discrete Fourier transform (DFT), as implemented by the fast Fourier transform (FFT), and are commonly preferred in consideration of their simple implementation and low computational burden. However, due to the FFT picket fence effect (PFE) and spectral leakage effect (SLE), significant bias would be introduced into the estimates of frequency, amplitude, and in particular, phase, if the sample frame is not an integer multiple of signal period [2,3]. Previous studies have shown that the maximum estimation error for the amplitude of a single sinusoid is up to −3.92 dB for the rectangle window and −1.42 dB for the Hanning window [4]. The absolute maximum error for a phase can be as high as π /2 [5]. In previous works [1,6–37], windowed interpolation algorithms were proposed in which interference caused by SLE was reduced by windowing, and errors introduced by PFE were eliminated by interpolation.

*

Corresponding author. E-mail address: jiufl[email protected] (J. Luo).

http://dx.doi.org/10.1016/j.dsp.2015.03.002 1051-2004/© 2015 Elsevier Inc. All rights reserved.

Several decades ago, Rife et al. [6] presented a frequency estimator based on the modulus of two FFT spectral lines. Jain et al. [7] reinvestigated the problem for the rectangle window. Later, it was extended for the Hanning window [8] and for Rife–Vincent class I windows [9,10]. In 1996, the same idea with a different explanation was developed by Xie and Ding in [11], in which both rectangle window and Hanning window were studied. Quinn proposed a frequency estimator based on interpolation using Fourier coefficients [12] in 1994 and later proposed a combined estimator to minimize variance [13]. Similar estimators were explored for Hanning-windowed data [14]. Macleod presented interpolators using three Fourier coefficients [15]. A comprehensive discussion on the use of Fourier coefficients to estimate frequency can be found in [16,17] and [18]. Aboutanios proposed two iterative interpolation algorithms, each of which converged asymptotically with a variance marginally above the Cramér–Rao bound (CRB) [19]. A more general form of the complex DFT coefficients based estimator was given by Liu et al. [20]. The estimators were further improved and extended for decaying sinusoids or exponentials [21–23]. The expressions of the estimators were similar to those in [24]. Diao and Meng considered weighted damped sinusoid signals [25], while Nandi and Kundu considered signals composed of multiple sinusoids [26]. In 2009, analytic interpolation formulas for maximum sidelobe decay windows (also known as Rife– Vincent class I windows) were uniformly derived by Belega and Dallet [27]. To improve the ability of long-range leakage rejection, weighted interpolated DFT (WIDFT) was proposed by Agrez [28], in which three spectral lines were used to improve the estimation result. Belega and Dallet [1] further developed the technique

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

119

x w (n) = x(n) w (n).

(3)

into a general form called the multipoint interpolated DFT approach (WMlpDFT). Recently, an average-based lpDFT was explored [29,30], in which the weighted average of two or more estimates, instead of a single estimate, was adopted to obtain the optimal variance of the frequency estimator in noise. An overview of frequency estimation methods was given in [38]. The previously described algorithms can only be applied to certain type(s) of windows, and typically only maximum sidelobe decay windows have been shown to have simple analytic solutions. There still remain many windows whose interpolation formulas are not described, such as Hamming window, Kasier–Bessel windows, Dolph–Chebyshev windows, to name a few [37]. The Hamming window is known for its good performance in rejection of adjacent interference. The Kasier–Bessel and Dolph–Chebyshev windows are known for superior performance in multi-tone detection and high noise immunity [37]. As a result, some approximation methods were presented for those windows. One method is the linear interpolation DFT (LIDFT), which was explored for a special window in [31] and metro-logically analyzed in [32]. LIDFT was then extended for classic windows in [33] achieved by the zero padding technique. Based on an idea from Offelli and Petri [10], Duda derived the polynomial approximation interpolation algorithm [37]. The two approximation methods require coefficients to be pre-calculated, which means that coefficients for different windows must be computed and stored in memory. Further, a matrix equation must be solved with LIDFT. These requirements are memory-intensive and inconvenient, especially for adjustable windows the characteristics of which can be adjusted by one or more parameters. The aim of this paper is to present an iterative DFT interpolation algorithm that can accurately estimate frequency in a discrete spectrum. The proposed algorithm is compatible with arbitrary windows, is straightforward to implement, and does not require parameter pre-calculation. An accurate frequency estimate can be obtained after two or three iterations. Moreover, incorrect polarity estimation (IPE) [14,16,35] can almost be completely avoided in the proposed algorithm. The rest of the paper is organized as follows. Section 2 gives a summary of analytic DFT interpolation formulas. In Section 3, the basic theory of the proposed algorithm is described. In Section 4, computer simulations are used to verify the effectiveness and accuracy of the algorithm. Systematic errors for different windows are shown and the performance of the proposed algorithm in the presence of white Gaussian noise is investigated. Incorrect polarity estimation is also discussed. Conclusions are given in Section 5.

Based on the DFT definition, the DFT coefficients of the weighted signal x w (n) can be computed by

X w (k) =

N −1 

x w (n)e − j

2π N

kn

,

0 ≤ k ≤ N − 1.

Substituting (3) into (4) yields

X w (k) =

A0 2

e j ϕ0 W N (k − λ0 ) +

A 0 − j ϕ0 W N (k + λ0 ). e 2

λ0 = l + δ,

(6)

where l and δ (−0.5 ≤ δ < 0.5) are the integer and the fractional parts, respectively. If signal-to-noise ratio (SNR) is well above the so-called “threshold” (approximately −18 to −16 dB) [39], l can be readily and correctly located because of its largest magnitude as given by | X w (l)|. It is worth noting that, if SNR is below the “threshold”, the largest magnitude may be badly distorted by noise, which causes spurious peaks or outliers. In this case, the returned local maximum may occur at frequencies far from l, and the coarse search routine will yield large errors. Consequently, only an SNR above the threshold is considered in this paper. If the index of the second largest bin is l + 1, then λ0 is greater than l (δ > 0). If the index of the second largest bin is l − 1, then λ0 is less than l (δ < 0). The goal of the interpolation algorithm is to obtain an ˆ (or δˆ ). accurate estimate of λ0 (or δ ), given by λ Substituting (6) into (5), the largest magnitude is given by

       X w (l) =  A 0 e j ϕ0 W N (−δ) + A 0 e− j ϕ0 W N (2l + δ).   2 2

x(n) = A 0 cos 2π

f0 fs

 n + θ0 ,

n = 0, 1 , 2 , . . . , N − 1

(1)

where A 0 is amplitude, f 0 is frequency, and θ0 is the phase angle. f s is the sampling frequency and n is the sample index. It is assumed that f s exceeds 2 f 0 in order to satisfy the Nyquist theorem. When N samples are acquired, the ratio between frequency, f 0 , and sampling frequency, f s , can be expressed as:

f0 fs

=

λ0 N

,

(7)

The second term on the right in (7) represents the contribution of leakage from the imaginary part of the spectrum. It is assumed from now on that λ0 satisfies λ0 > 5 and λ0 < N /2 − 5 to ensure that leakage coming from the negative part of the spectrum is minimal and can be ignored. With this assumption, (7) reduces to

     X w (l) = A 0  W N (−δ).

(8)

Similarly, the second largest magnitude is

In order to explain the general theory of interpolation algorithms, let us consider, for simplicity but without loss of generality, the discrete time signal in the form



(5)

In Eq. (5), W N (k) is the Discrete Time Fourier Transform (DTFT) of the window w (n). If x(n) does not contain an integer number of periods, the normalized frequency λ0 lies between two DFT bins. Therefore, λ0 can be further expressed as:

2

2. Theoretical background

(4)

n =0

(2)

where λ0 is the number of the acquired cosine-wave cycles. Note that the parameter λ0 also represents the normalized frequency expressed in bins. After multiplying the signal samples, x(n), by the data window values, w (n), we obtain the weighted samples

     X w (l ± 1) = A 0  W N (−δ ± 1). 2

(9)

In the two-point (2p) interpolation algorithm, the ratio of the two largest magnitudes is expressed as

α=

| X w (l ± 1)| W N (−δ ± 1) = . | X w (l)| W N (δ)

(10)

Eq. (10) shows that the ratio α depends only on the selected window and the frequency deviation δ . Therefore, if the data window is known, δ can be determined solely by α . For maximum sidelobe decay windows [27], δ can be written explicitly as

δ=±

H α − ( H − 1)

α+1

,

(11)

where H denotes for the number of terms in the maximum sidelobe decay window. For other windows, the polynomial approximation is suggested [37].

120

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

In the three-point (3p) interpolation algorithm, the ratio of the three largest magnitudes is defined as

α=

| X w (l)| + | X w (l − 1)| | W N (−δ)| + | W N (−δ − 1)| = . | X w (l)| + | X w (l + 1)| | W N (−δ)| + | W N (−δ + 1)|

(12)

Similarly, when the maximum sidelobe decay windows are used, δ can be explicitly obtained [27]

δ=H

1−α 1+α

(13)

.

The weighted multipoint interpolated DFT, which is a generalization of the three-point interpolation algorithm, is described in [27]. There also exist variants based on maximum sidelobe decay windows. For the sake of conciseness, the complex DFT coefficientbased estimators and the variants of the modulus-based estimators are not reiterated here. Once δ is determined, the normalized frequency can be computed and the frequency f 0 can be determined using (2). Amplitude and phase can be calculated using the corresponding formula derived in [37]. As discussed earlier, simple interpolation formulas are only valid for maximum sidelobe decay windows. In practice, it is often desired to use Rife–Vincent class II windows or Rife–Vincent class III windows. Non-cosine windows or adjustable windows are also commonly used, such as Kaiser–Bessel and Dolph–Chebyshev windows. Windows with a narrow mainlobe and low sidelobes are important for detecting closely spaced frequency components, particularly in cases of short data records [37]. Studies show windows with a narrow mainlobe also have better noise immunity [4]. Windows with rapidly attenuated sidelobes are often preferred to compress the leakage from components far away [4]. Therefore, it is important to design a simple algorithm that is applicable for various windows. 3. Proposed algorithm In Section 2, it was already shown that the coarse search returns l, which is the index of the bin with the largest magnitude. Further, we define the Fourier coefficient

V q (m) =

N −1 

x w (n)e − j2π

lm +q n N

(14)

,

n =0

where l0 = l. According to the definition of (14), we have | X w (l)| = | V 0 (0)|, | X w (l + 1)| = | V 1 (0)| and | X w (l − 1)| = | V −1 (0)|. Therefore, we can obtain an estimate l1 = l0 + δ , where δ is determined by (10) and (11) with H = 2. l1 is close to the true frequency when the Hanning window is selected [8]. Other windows will introduce an undesirable error. To improve accuracy, we can use formula (11) again in which the ratio defined in (10) is replaced by α = | V 0.5 (1)|/| V −0.5 (1)|(| V 0.5 (1)| ≤ | V −0.5 (1)|) or α = | V −0.5 (1)|/| V 0.5 (1)|(| V 0.5 (1)| > | V −0.5 (1)|). The two Fourier coefficients V 0.5 (1) and V −0.5 (1) are computed according to (14). A new frequency estimate is obtained by

l2 = l1 − 0.5 +

    V 0.5 (1) ≤  V −0.5 (1) ,

(15a)

    V 0.5 (1) >  V −0.5 (1) .

(15b)

2α − 1 

α+1

or

l2 = l1 + 0.5 −

2α − 1 

α+1

Substituting α into (15a), (15b) and after some algebraic manipulations, we can achieve that

l2 = l1 +

3 | V 0.5 (1)| − | V −0.5 (1)| 2 | V 0.5 (1)| + | V −0.5 (1)|

,

(16)

for both (15a) and (15b). Further iterations can be generalized as

l m +1 = l m +

3 | V 0.5 (m)| − | V −0.5 (m)| 2 | V 0.5 (m)| + | V −0.5 (m)|

,

m = 1, 2, 3, . . . ,

(17)

until |lm − lm+1 | < τ . Finally, a precise estimate can be obtained through the iterative interpolation DFT algorithm (IIpDFT). Eq. (17) seems to be similar to Algorithm 2 proposed in [19], with the difference being that Eq. (17) is based on the Hanning window. However, it should be pointed out that the algorithm derived by Aboutanios and Mulgrew in [19] did not involve the technique of windowing, and all the derivations were based on the rectangle window. Nonetheless, in spectrum analysis windowing is extensively adopted to suppress spectral leakage. In addition, the initial conditions of the two algorithms are different. When Eq. (17) was directly used for other classic windows, there were some problems. First, the iterative algorithm is only applicable for some windows because of the convergence problem. For windows with large scalloping loss (SL), the algorithm may diverge. Scalloping loss is defined as the ratio of coherent gain for a frequency component located halfway between DFT bins to the coherent gain for a frequency component located exactly at a DFT bin [27],

SL =

| W (0.5)| , W (0)

(18)

where W (·) is the DTFT of window w (n) expressed in bins. Second, convergence is slow with windows that have low scalloping loss (SL < 1.30 dB) or high scalloping loss (SL > 1.60 dB). This increases the computation time required to converge to the true signal frequency. In order to overcome these shortcomings, an improved iterative algorithm (AcIIpDFT) is proposed, in which a suitable power of spectral lines was adopted. It is already known that an accurate frequency estimate can be obtained by Eq. (17) when the Hanning window is used. If we can make the spectrum of other windows similar to that of the Hanning window, it would greatly increase the convergence using Eq. (17). We modify the ratio α in (10) and the iterative algorithm in (17) as shown below,

α=

| X w (l ± 1)|κ | X w (l)|κ

(19)

and

l m +1 = l m +

3 Y 0.5 (m) − Y −0.5 (m) 2 Y 0.5 (m) + Y −0.5 (m)

,

m = 1, 2, 3, . . . ,

(20)

where Y q (m) = | V q (m)|κ and κ = SLHan /SL w (n) . SLHan and SL w (n) represent the scalloping loss of the Hanning window and the corresponding window used in (3), respectively. This power can make the spectrum of any classic window and the spectrum of the Hanning window fit well with each other in the middle of their mainlobes so that convergence can be accelerated. Figs. 1 and 2 show iterations in IIpDFT as a function of δ for various windows, including Riesz window (Riesz), Boham window (Boham), Blackman window (Blackman), minimum-3 term window (Min-3-term), minimum-4 term window (Min-4-term), Gaussian windows (GS), Kaiser–Bessel windows (KB), Dolph–Chebyshev windows (DC), Riemann window (Riemann), triangle window (Triangle), Hanning window (Hanning) and Hamming window (Hamming) [4,40]. The simulation was performed in the absence of noise. The frequency was scanned from 63.5 to 64.5 Hz in order to minimize the influence of negative frequency. Phase was set to zero. The number of samples was 128 and the sampling rate was 128 Hz. Figs. 3 and 4 show iterations using the improved algorithm for comparisons. Results indicate that the number of iterations required to achieve the stopping criterion (τ = 10−4 ) in AcIIpDFT are significantly less than in IIpDFT for the same window. The number of iterations in IIpDFT for most windows ranges from 5 to 10, and

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Fig. 1. Iterations for windows with low SL (IIpDFT,

τ = 10−4 ).

Fig. 2. Iterations for windows with high SL (IIpDFT,

τ = 10−4 ).

Fig. 3. Iterations for windows with low SL (AcIIpDFT,

τ = 10−4 ).

Fig. 4. Iterations for windows with high SL (AcIIpDFT,

τ = 10−4 ).

121

122

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Fig. 5. Frequency error (M = 0).

Fig. 6. Frequency error (M = 1).

sometimes exceeds 10 to reach up to 30. In contrast, the number of iterations in AcIIpDFT for all windows is no more than 2 and satisfactory results can be obtained after only a few iterations. The improved algorithm accelerates convergence and decreases computation time. Moreover, the algorithm based on (18) and (19) can be applied to the rectangle window, which is not compatible with IIpDFT. Given that the rectangle window has the highest SL, we can reasonably infer that the improved algorithm is appropriate for all conventional windows. Although it is difficult to prove mathematically, this fact has been validated by numerical testing. Fig. 4 shows the number of iterations for the following windows: the rectangle window, Dolph–Chebyshev window (α = 1.5), and Kaiser–Bessel window (β = 0.2), which are all exceptions in IIpDFT. No more than 3 iterations are required to reach the stopping criteria, except with the rectangle window, which requires 4 to 6 iterations. The AcIIpDFT algorithm can be summarized as follows: (1) Perform a coarse search to determine l, the index of the bin with the largest magnitude. (2) Obtain an initial condition using the two-point interpolation algorithm according to (19) and (11) with H = 2, and arrive at l1 = l0 + δ , where l0 = l.

 N −1

lm +q

− j2π n N (3) Define V q (m) = , Y q (m) = | V q (m)|κ , n=0 x w (n)e κ = SLHan /SL w (n) and perform iterations:

l m +1 = l m +

3 Y 0.5 (m) − Y −0.5 (m) 2 Y 0.5 (m) + Y −0.5 (m)

,

m = 1, 2, 3, . . .

(4) Stop iterations when |l M − l M +1 | < τ , where τ is the stopping criterion determined by the required accuracy. ˆ = lM . (5) Finally, we have λ

4. Simulation results In the previous section, we derived formulas that produce an estimate of frequency. In order to verify the effectiveness and accuracy of the proposed algorithm, results of the above-mentioned procedure are applied to simulated signals. First, the dependency of precision on the number of iterations is discussed. Furthermore, we analyze systematic errors and noise sensitivity. Finally, we investigate the problem of incorrect polarity estimation occurring at a nearly synchronous sampling condition. Only the accelerated algorithm is considered in all simulations. The cosine wave amplitude, A, was set to 1, the number of analyzed samples, N, was set to 1024, and the sampling rate was set to 1024 so that frequency resolution was 1. For the sake of brevity, frequency estimates were scaled by frequency resolution and expressed in bins. 4.1. Influence of iterations As stated in Section 3, the number of iterations required in AcIIpDFT to achieve the stopping criterion (τ = 10−4 ) for all described windows is no more than three, except for the rectangle window. Note that the stopping criterion is not the same as the estimation error. We now study the estimation error of frequency for different numbers of iterations, which provides a foundation for practical implementation. δ is varied from −0.5 to 0.5 in steps of 0.05. For each δ , we generated 72 instances with different phases. This means that a phase scan was done in the interval (−π , π ] in steps of π /36. The maximum absolute value of the difference between true and estimated frequency was selected as the maximum frequency estimation error (MFEE). The estimation error was investigated for zero to three iterations and various windows as shown in Figs. 5–8.

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

123

Fig. 7. Frequency error (M = 2).

Fig. 8. Frequency error (M = 3).

Fig. 5 shows MFEE for different windows with no iterations. For Hanning window, the remaining error is quite small (less than 10−10 ), and for the Hamming window the error is less than 10−3 . The MFEE for the rest of the windows is approximately 0.02 at δ = 0. After one iteration, MFEE was notably reduced for all windows except for Hamming window a Hanning window, in which the error remained at the same level as shown in Fig. 6. With further iterations (Figs. 7 and 8), the error for certain windows continued to decrease (e.g., Blackman window and Minimum 4-term window). The residual error at this point is primarily due to interference from negative frequency. Therefore, three iterations appear to be adequate for practical implementation. 4.2. Systematic errors Certain approximations were necessary in order to derive the frequency estimation algorithm. The proposed method neglects the contribution from the negative frequency axis (Eq. (7)). In addition, as discussed in Section 4.1, the finite iterations would result in residual error, although it was small after a few iterations. As a consequence, the approximations introduced systematic modeling errors in the estimate, even in the noiseless case. The definition of systematic errors (SME) was similar to that in [41] but with a few modifications. For a certain frequency bin k, a frequency scan in the interval [k − 0.5, k + 0.5), in steps of 0.05, and a phase scan in the interval (−π , π ] in steps of π /36, was performed individually. Therefore, for each frequency bin 1440 test signals were generated. The maximum absolute difference value was the SME. Figs. 9–11 show the systematic frequency estimation errors for different windows as a function of the frequency bin. In all cases, the number of iterations was set to three. Fig. 9 shows that errors for the Hanning window are nearly the same as in [41]. Note that errors for the Blackman window become nearly a straight

line in the interval [150, 350]. This is because in this range estimation errors are determined by finite iterations and the error from negative frequency is less significant. Systematic errors of the Kaiser–Bessel windows with different selected values of parameter β are shown in Fig. 10. They are similar to the systematic errors of the 2p interpolation algorithm described in [37]. Straighttype errors for the Kaiser–Bessel windows are also apparent when β = 4.90 and β = 5.92. For Dolph–Chebyshev windows, we chose the attenuation of the sidelobes at 50, 70, 90, 110, and 130 dB, respectively. Corresponding systematic errors are shown in Fig. 11. The errors are nearly straight lines due to the uniform structure of the sidelobes. Frequency estimation errors dramatically increased when frequency was near the zero frequency or Nyquist frequency because of serious interference from the negative component. For general cases, systematic error is very small and can be neglected. 4.3. Multiple frequency components In this subsection, we considered signals containing the resolved multiple frequency components. To maintain simplicity but without loss of generality, a double frequency component model was studied, in which a frequency component λ1 gradually approaches λ2 = 255.5 in steps of 1/8. For each frequency step, a phase scan was performed in the range (−π , π ] in steps of π /36, for both λ1 and λ2 , yielding 5184 independent trials per step. The worst-case scenario was selected as the result for each step. Figs. 12 and 13 show frequency errors of λ1 and λ2 as a function of frequency difference λ. The two components show nearly the same results when signals were weighted by the same window. Figures reveal that the algorithm has poor performance when the frequency difference is small. This is due to the serious interference from window sidelobes. Even mainlobe interference is possible if the two frequency components are close enough to

124

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Fig. 9. SME for different windows (M = 3).

Fig. 10. SME for Kaiser–Bessel windows (M = 3).

Fig. 11. SME for Dolph–Chebyshev windows (M = 3).

each other. If multiple frequency components are contained in the signals, the proposed approach can be applied independently on each detected spectral peak, as long as the frequency difference between the adjacent spectral components is larger than half of the window mainlobe width. Interference from adjacent components can be significantly reduced by selecting proper windows with low sidelobes (i.e., minimum 4-term window) or with sidelobes that have a sharp falloff rate (i.e., Hanning window). 4.4. Influence of noises In the above subsections, we considered the influence of iteration number and leakage effect. In a practical measurement scenario, theoretical signals are usually distorted by noise from various uncontrollable factors. As a consequence, estimation accuracy is greatly reduced and the frequency estimate becomes a ran-

dom variable. In order to assess the algorithm’s sensitivity to white noise, signals disturbed by additive zero-mean Gaussian noise were generated. The noise level was adjusted to acquire test signals with SNR ratios ranging from −5 to 120 dB. The SNR is defined as





SNR = 10 log10 A 2 / 2σ 2 ,

(21)

where A denotes the amplitude of the component and σ 2 is the variance of white noise. True frequency was randomly selected in the range [255.5, 256.5). For each selected frequency, a phase scan was done, and the worst-case scenario was selected as the result. In total, 10,000 instances of such test signals were generated for each value of SNR, and the standard deviation (STD) of frequency estimation errors was evaluated. Figs. 12–14 show standard deviation of frequency estimation errors for different windows as a function of SNR. The Cramér–Rao lower bound (CRLB) for fre-

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Fig. 12. Frequency error for λ1 (M = 3).

Fig. 13. Frequency error for λ2 (M = 3).

Fig. 14. Standard deviation of errors for different windows (M = 3).

Fig. 15. Standard deviation of errors for Dolph–Chebyshev windows (M = 3).

125

126

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Fig. 16. Standard deviation of errors for Kaiser–Bessel windows (M = 3).

Fig. 17. Illustration of the IPE problem for 2p-interpolation algorithm.

quency estimation is also shown for comparison. The CRLB for frequency estimation of the signal (1) is given in [37]. Figs. 14–16 show that the standard deviation of frequency estimation errors decreases as SNR increases. Results from the new algorithm are in accordance with those in [22]. As shown in the figures, when SNR is low, noise is the determining factor of uncertainty, and the influence of systematic errors can be neglected. In contrast, when SNR is high, systematic errors become more pronounced. The value of SNR inflection depends on different windows. Results also reveal that at very low SNRs, those windows with a narrow mainlobe have better performance than those with a wide mainlobe. This is consistent with the equivalent noise bandwidth (ENBW) which is an assessment of window noise properties. 4.5. Study of incorrect polarity estimation An intrinsic weakness in traditional interpolation algorithms is the problem of incorrect polarity estimation [14,16,35,42]. It often occurs at a nearly coherent sampling condition (when |δ| is near zero) for 2p-interpolation based algorithms and at the worst case of noncoherent sampling condition (when |δ| is near 0.5) for 3pinterpolation based algorithms. Under such sampling conditions, correctly locating the correct spectral line is often difficult because of wideband noise or spectral leakage. For a 2-point based algorithm, if the spectral line is selected incorrectly, the true frequency bin will not be embraced in the two largest bins. Such erroneous locations will not only affect the value of α but also result in sign reversal of δ . It will lead to a significant error in the frequency estimate if the value of δ is considerable, as illustrated in Fig. 17 for the case of a 2p-interpolation based algorithm. In the case of low SNR, the ratio of mistaken location increases significantly. The ratios of mistaken locations of the lower spectral line for various windows under different noise conditions are shown in Figs. 18

and 19. Fig. 18 shows that significant incorrect locations are apparent when SNR is −5 dB, especially when |δ| is near zero. As SNR increases to 10 dB in Fig. 19, the problem is alleviated, but a high possibility of error still remains when |δ| is close to zero. Consequently, in the noise condition, the traditional 2p-interpolation based algorithms become vulnerable when a small value of |δ| is encountered. Similar phenomena are found when |δ| is near 0.5 for 3p-interpolation based algorithms. However, incorrect polarity estimations can be avoided in the proposed algorithm because iterations can automatically adjust the estimation polarity. In order to assess the algorithms against the IPE problem, theoretical signals with additive zero-mean Gaussian noise were generated. SNR was set to −5 dB so that a high possibility of incorrect polarity estimation would occur. Simulation frequency was scanned from 255.5 to 256.5, in steps of 0.025. Phase was randomly selected. For each frequency, 10,000 independent instances were produced. Samples were weighted by the Hanning window before DFT. Fig. 20 shows the standard deviation of estimates (STD) as a function of the frequency deviation δ . In addition, another five algorithms, including the estimator proposed by Grandke in 1983 (Grandke), the 2-point based estimator (Quinn1) and the 3-point based optimal estimator (Quinn2) proposed by Quinn in 2006, the estimator proposed by Agrez in 2002 (Agrez) and the estimator proposed by Macleod in 1998 (Macleod) are also displayed for comparison. The figures are very informative regarding the impact on frequency estimation due to the IPE problem. We can clearly see that the Grandke and Quinn1 have almost identical performance. As expected, STD increases with decreasing |δ|. This trend matches that of the ratio of mistaken locations. This confirms the fact that the IPE has a significant influence on the accuracy and uncertainty of frequency estimation. The 3p-interpolation based algorithms (Quinn2, Agrez and Macleod) also suffer from the defects. The worst case occurs when |δ| ≈ 0.5, while the performance is bet-

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Fig. 18. Ratio of mistaken locations in 10,000 independent trials (SNR = −5).

Fig. 19. Ratio of mistaken locations in 10,000 independent trials (SNR = 10).

Fig. 20. Standard deviation of estimates with Hanning window (SNR = −5).

Fig. 21. Comparison between different algorithms with rectangle window.

127

128

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

ter when |δ| ≈ 0. In contrast, the STD of the proposed algorithm is nearly flat in the entire range of δ . The behavior is almost independent of δ . It is implied that the proposed algorithm appears insensitive to IPE problem. Moreover, the proposed algorithm shows the smallest STD compared with all other algorithms. It has an overwhelming advantage over the traditional algorithms under noise condition. From Figs. 18 and 19, it is convinced that the conclusion also holds for other windows. To confirm this fact more clearly, we compared the proposed algorithm (AcIIpDFT(Rect)) with Quinn’s estimators including the 2-point based estimator (Quinn1(Rect)) [12] and the 3-point based optimal estimator (Quinn2(Rect)) [13] when rectangle window is used. The STDs of different estimators scaled by CRLB were shown as a function of SNR in Fig. 21. The simulation parameters are in full accord with those in Section 4.4. We can see that, at low SNRs, the proposed algorithm provides better results than Quinn’s algorithms and at high SNRs, the optimal estimator has better performance. As a result, the proposed estimator is more robust than traditional estimators when noise is the determining factor of uncertainty because of its strong resistance to IPE. 5. Conclusion This paper proposes an interpolation algorithm that provides accurate frequency estimation. The algorithm is based on iterative interpolation of DFT coefficients and can estimate frequency using only a few iterations. The proposed iterative algorithm is straightforward to understand, has the advantage of practical simplicity, keeps the cost of computation time at a tolerable level, and is compatible with almost all conventional windows. Furthermore, no precalculation is required and there is no need to know the spectrum of the data window. Efficacy and accuracy of the proposed algorithm were evaluated for different windows by simulation. The influence of negative frequency and white Gaussian noise was studied, particularly for Dolph–Chebyshev and Kaiser–Bessel windows. The problem of incorrect polarity estimation for traditional interpolation algorithms was discussed. Simulation results show that the proposed algorithm is robust against mistaken location of the spectral line, and consequently has stronger resistance to additive noise compared with current interpolation algorithms. References [1] D. Belega, D. Dallet, D. Petri, Accuracy of sine wave frequency estimation by multipoint interpolated DFT approach, IEEE Trans. Instrum. Meas. 59 (2010) 2808–2815. [2] F. Zhang, Z. Geng, W. Yuan, The algorithm of interpolating windowed FFT for harmonic analysis of electric power system, IEEE Trans. Power Deliv. 16 (2001) 160–164. [3] I. Santamaria-Caballero, C.J. Pantaleon-Prieto, J. Ibanez-Diaz, E. Gómez-Cosío, Improved procedures for estimating amplitudes and phases of harmonics with application to vibration analysis, IEEE Trans. Instrum. Meas. 47 (1998) 209–214. [4] F.J. Harris, On the use of windows for harmonic analysis with the discrete Fourier transform, Proc. IEEE 66 (1978) 51–83. [5] D.S. Huang, Phase error in fast Fourier transform analysis, Mech. Syst. Signal Process. 9 (1995) 113–118. [6] D.C. Rife, G. Vincent, Use of the discrete Fourier transform in the measurement of frequencies and levels of tones, Bell Syst. Tech. J. 49 (1970) 197–228. [7] V.K. Jain, W.L. Collins, D.C. Davis, High-accuracy analog measurements via interpolated FFT, IEEE Trans. Instrum. Meas. 28 (1979) 113–122. [8] T. Grandke, Interpolation algorithms for discrete Fourier transforms of weighted signals, IEEE Trans. Instrum. Meas. 32 (1983) 350–355. [9] G. Andria, M. Savino, A. Trotta, Windows and interpolation algorithms to improve electrical measurement accuracy, IEEE Trans. Instrum. Meas. 38 (1989) 856–863. [10] C. Offelli, D. Petri, Interpolation techniques for real-time multifrequency waveform analysis, IEEE Trans. Instrum. Meas. 39 (1990) 106–111. [11] M. Xie, K. Ding, Corrections for frequency, amplitude and phase in a fast Fourier transform of a harmonic signal, Mech. Syst. Signal Process. 10 (1996) 211–221.

[12] B.G. Quinn, Estimating frequency by interpolation using Fourier coefficients, IEEE Trans. Signal Process. 42 (1994) 1264–1268. [13] B.G. Quinn, Estimation of frequency, amplitude, and phase from the DFT of a time series, IEEE Trans. Signal Process. 45 (1997) 814–817. [14] B. Quinn, Frequency estimation using tapered data, in: Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP 2006, IEEE, 2006, pp. 73–76. [15] M.D. Macleod, Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones, IEEE Trans. Signal Process. 46 (1998) 141–148. [16] B.G. Quinn, E.J. Hannan, The Estimation and Tracking of Frequency, Cambridge University Press, 2001, pp. 180–206. [17] B.G. Quinn, Recent advances in rapid frequency estimation, Digit. Signal Process. 19 (2009) 942–948. [18] B.G. Quinn, The estimation of frequency, in: T.S. Rao, S.S. Rao, C.R. Rao (Eds.), Handbook of Statistics: Time Series Analysis: Methods and Applications, Elsevier Press, 2012, pp. 585–621. [19] E. Aboutanios, B. Mulgrew, Iterative frequency estimation by interpolation on Fourier coefficients, IEEE Trans. Signal Process. 53 (2005) 1237–1242. [20] Y. Liu, Z. Nie, Z. Zhao, Q.H. Liu, Generalization of iterative Fourier interpolation algorithms for single frequency estimation, Digit. Signal Process. 21 (2011) 141–149. [21] E. Aboutanios, Generalised DFT-based estimators of the frequency of a complex exponential in noise, in: 2010 3rd International Congress on Image and Signal Processing, CISP, IEEE, 2010, pp. 2998–3002. [22] E. Aboutanios, Estimation of the frequency and decay factor of a decaying exponential in noise, IEEE Trans. Signal Process. 58 (2010) 501–509. [23] E. Aboutanios, Estimating the parameters of sinusoids and decaying sinusoids in noise, IEEE Instrum. Meas. Mag. 14 (2011) 8–14. [24] M. Bertocco, C. Offelli, D. Petri, Analysis of damped sinusoidal signals via a frequency-domain interpolation algorithm, IEEE Trans. Instrum. Meas. 43 (1994) 245–250. [25] R. Diao, Q. Meng, An interpolation algorithm for discrete Fourier transforms of weighted damped sinusoidal signals, IEEE Trans. Instrum. Meas. 63 (6) (2014) 1503–1515. [26] S. Nandi, D. Kundu, An efficient and fast algorithm for estimating the parameters of sinusoidal signals, Sankhya 68 (2006) 283–306. [27] D. Belega, D. Dallet, Multifrequency signal analysis by interpolated DFT method with maximum sidelobe decay windows, Measurement 42 (2009) 420–426. [28] D. Agrez, Weighted multipoint interpolated DFT to improve amplitude estimation of multifrequency signal, IEEE Trans. Instrum. Meas. 51 (2002) 287–292. [29] D. Belega, Accuracy analysis of the normalized frequency estimation of a discrete-time sine-wave by the average-based interpolated DFT method, Measurement 46 (2013) 593–602. [30] K.F. Chen, Y.F. Li, Combining the Hanning windowed interpolated FFT in both directions, Comput. Phys. Commun. 178 (2008) 924–928. [31] J. Borkowski, LIDFT—the DFT linear interpolation method, IEEE Trans. Instrum. Meas. 49 (2000) 741–745. [32] J. Borkowski, J. Mroczka, Metrological analysis of the LIDFT method, IEEE Trans. Instrum. Meas. 51 (2002) 67–71. [33] J. Borkowski, J. Mroczka, LIDFT method with classic data windows and zero padding in multifrequency signal analysis, Measurement 43 (2010) 1595–1602. [34] E. Jacobsen, P. Kootsookos, Fast, accurate frequency estimators [DSP tips & tricks], IEEE Signal Process. Mag. 24 (2007) 123–125. [35] Y.F. Li, K.F. Chen, Eliminating the picket fence effect of the fast Fourier transform, Comput. Phys. Commun. 178 (2008) 486–491. [36] K.F. Chen, S.L. Mei, Composite interpolated fast Fourier transform with the Hanning window, IEEE Trans. Instrum. Meas. 59 (2010) 1571–1579. [37] K. Duda, DFT interpolation algorithm for Kaiser–Bessel and Dolph–Chebyshev windows, IEEE Trans. Instrum. Meas. 60 (2011) 784–790. ´ K. Duda, Frequency and damping estimation methods—an overview, [38] T. Zielinski, Metrol. Meas. Syst. 18 (2011) 505–528. [39] D. Rife, R.R. Boorstyn, Single tone parameter estimation from discrete-time observations, IEEE Trans. Inf. Theory 20 (1974) 591–598. [40] A. Nuttall, Some windows with very good sidelobe behavior, IEEE Trans. Acoust. Speech Signal Process. 29 (1981) 84–91. [41] J. Schoukens, R. Pintelon, H. Van Hamme, The interpolated fast Fourier transform: a comparative study, IEEE Trans. Instrum. Meas. 41 (1992) 226–232. [42] E. Aboutanios, A modified dichotomous search frequency estimator, IEEE Signal Process. Lett. 11 (2004) 186–188.

Jiufei Luo was born in Sichuan Province of China in 1987. He received the M.S. degree in biomedical engineering from Chongqing University of Technology, Chongqing city of China, in 2012. He is currently a Doctor in the Laboratory for Digital Measurement Systems at the Faculty of Mechanical Engineering of Chongqing University. His fields of interest include applications of digital signal processing and analysis, equipment condition monitoring and fault diagnosis.

J. Luo et al. / Digital Signal Processing 41 (2015) 118–129

Zhijiang Xie was born in Hunan Province of China in 1963. After undergraduate course graduation in 1983, he received M.S. degree and the Ph.D. degree in mechanical engineering from Chongqing University, where he is currently a Professor with the Mechanical Engineering School of Chongqing University. His research activities, focus on mechanotronics, mechanical–electrical integration, digital and analog signal processing, mechanical and electrical integration equipment condition monitoring and

129

fault diagnosis. His interests also include machine vision, as well as computer aided test theory and technology. Ming Xie received the M.S. degree in mechanical engineering from Chongqing University in 1988. He is currently an associate professor in Chongqing University of Technology. His fields of interest include applications of digital signal processing and analysis.