Generalization of interpolation DFT algorithms and frequency estimators with high image component interference rejection

Abstract

This paper focuses on the problem of frequency estimation of noise-contaminated sinusoidal. A basic tool to solve this problem is the interpolated discrete Fourier transform (DFT) algorithms, in which the influences of the spectral leakage from negative frequency are often neglected, resulting in significant errors in estimation when the signals contained small cycles. In this paper, analytic expressions of the interference due to the image component are derived and its influences on the traditional two-point interpolated DFT algorithms are analyzed. Based on the achieved expressions, the interpolated DFT algorithms are generalized and a novel frequency estimator with high image component interference rejection is proposed. Simulation results show that the frequency errors returned by the new algorithm are very small even though only one or two cycles are obtained. Comparative studies indicate that the new algorithm also has a good performance in the noise condition. With the advantages of high precision and strong robustness against additive noise, the proposed algorithm is a good choice for frequency estimation when the negative frequency interference is the dominant error source.

1 Introduction

Spectral analysis based on the discrete Fourier transform (DFT) and implemented by the fast Fourier transform (FFT) has been widely used in many fields for several decades. However, there are two drawbacks in the classic version: the spectral leakage effect (SLE) caused by the lack of periodicity and the picket fence effect (PFE) due to the frequency sampling [1]. These effects will lead to significant errors in spectral analysis such as parameter estimation [2, 3]. In order to obtain accurate estimates of signal parameters, a lot of solutions were proposed [4–14]. The interpolation discrete Fourier transform (IpDFT) algorithm is one of the most popular algorithms.

Early in 1970s, interpolation algorithms based on the moduli of two FFT spectral bins were presented by Rife et al. [4]. Various improved algorithms have been put forward in the following decades, such as the weighted interpolated DFT (WIDFT) proposed by Agrez [9] and the multi-point interpolated DFT approach (WMlpDFT) by Belega and Dallet [12]. Specifically, simple analytical solutions can be obtained when the maximum sidelobe decay windows (MSDW, also known as Rife–Vincent class I windows) are adopted [11]. However, the algorithms mentioned above are all established on a very important assumption that the leakage coming from the image component plays a minor role and could be ignored [11–13]. In fact, if signals contain a small number of cycles, the negative frequency component usually will exercise great influences on the estimators [15]. Although the multi-point IpDFT methods can reduce the sensitivity to some extent, frequency estimation errors due to the spectral interference from the image component still remain significant [15], especially for the signals with a few cycles. Furthermore, it should be emphasized that the reduction of systematic errors by the weighted interpolated DFT or multi-point IpDFT methods is at the expense of worse noise properties [9, 12]. In engineering practice, noise properties usually are more important than systematic errors. The accuracy could be improved with longer records as well, but the cost is increased response time [16]. As a result, we often have to make a tradeoff between the estimation accuracy and the overall system responsiveness. Consequently, it is of great significance to propose a new and simple algorithm by which accurate parameter estimates can be achieved with short intervals, especially for those fields where real-time response is required [16–18].

Recently, Belega et al. proposed an improved three-point IpDFT which exhibits a high rejection capability with respect to the interference from negative frequency [15]. In this paper, we proposed a novel frequency estimator by which the leakage coming from image component can be further reduced compared with the algorithm proposed by Belega et al. More importantly, it also keeps good noise properties due to a two-point-based mechanism. The remaining parts are organized as follows. In Section 2, we present a concise summary of the traditional interpolation algorithms. Symbols and basic equations used throughout the paper are defined in this section as well. In Section 3, the analytical expressions of the interference from the negative frequency are derived. With the properties of the derived expressions, we generalize the interpolation DFT algorithms. In Section 4, the frequency estimators with high interference rejection of image component are proposed. In Section 5, some computer simulations are carried out and the performance of the proposed algorithms are compared with some other state-of-the-art IpDFT methods. Finally, main conclusions are drawn in Section 6.

2 Theoretical background

In order to explain the basis of the frequency-estimating procedure, first, let x _raw(n) be the samples of a discrete cosine wave in the form

$$ {x}_{\mathrm{raw}}(n)={A}_0\kern0.5em \cos \left(2\pi {f}_0n\varDelta T+{\varphi}_0\right)\kern3em n=0,1,2,\dots $$

(1)

where A ₀, f ₀, and φ ₀ are the amplitude, frequency, and phase, respectively. ΔT denotes the sampling interval and n is the index of the samples. Sampling rate f _s  = 1/ΔT is supposed to fulfill the Nyquist sampling theorem so that aliasing of spectrum does not occur. When N samples are acquired, f ₀ is normalized by the frequency resolution Δf = f _s /N and is expressed as

$$ {\lambda}_0=\frac{f_0}{\varDelta f}=\frac{f_0}{f_s/N}, $$

(2)

where λ ₀ is the normalized frequency expressed in bins [12]. Usually, samples are weighted by a window function w _N (n) before DFT

$$ x(n)={x}_{\mathrm{raw}}(n){w}_N(n) $$

(3)

The DFT of N weighted samples at the spectral line k is given by

$$ X(k)=\frac{A_0}{2}{e}^{j{\varphi}_0}{W}_N\left(k-{\lambda}_0\right)+\frac{A_0}{2}{e}^{-j{\varphi}_0}{W}_N\left(k+{\lambda}_0\right), $$

(4)

where W _N (·) denotes the DTFT of the window w _N (n). If the bin number with the largest magnitude is $l$, then the largest magnitude is given by

$$ \left|X(l)\right|=\left|\frac{A_0}{2}{e}^{j{\varphi}_0}{W}_N\left(l-{\lambda}_0\right)+\frac{A_0}{2}{e}^{-j{\varphi}_0}{W}_N\left(l+{\lambda}_0\right)\right| $$

(5)

The second largest magnitude is given by

$$ \left|X\left(l\pm 1\right)\right|=\left|\frac{A_0}{2}{e}^{j{\varphi}_0}{W}_N\left(l-{\lambda}_0\pm 1\right)+\frac{A_0}{2}{e}^{-j{\varphi}_0}{W}_N\left(l+{\lambda}_0\pm 1\right)\right|, $$

(6)

where it takes the negative if l − λ ₀ > 0 and the positive if l − λ ₀ < 0. The second terms on the right in (5) and (6) represent the image component in the spectrum. The interference is very small compared with the first term as long as l is far from zero frequency and Nyquist frequency. In this situation, W _N (l + λ ₀) can be ignored and (5) is reduced to

$$ \left|X(l)\right|=\frac{A_0}{2}\left|{W}_N\left(l-{\lambda}_0\right)\right| $$

(7)

Similarly, (6) can be rewritten as

$$ \left|X\left(l\pm 1\right)\right|\cong \frac{A_0}{2}\left|{W}_N\left(l-{\lambda}_0\pm 1\right)\right| $$

(8)

For the two-point (2p) interpolation algorithm [6, 8, 11], we introduce α defined as the ratio of the two largest magnitudes

$$ \alpha =\frac{\left|X\left(l\pm 1\right)\right|}{\left|X(l)\right|}\cong \frac{W_N\left(l-{\lambda}_0\pm 1\right)}{W_N\left(l-{\lambda}_0\right)} $$

(9)

It can be clearly seen from Eq. (9) that the ratio α only depends on the normalized frequency λ ₀ if the data window is already known. In particular, λ ₀ can be determined by simple and explicit forms when the maximum sidelobe decay windows are used [11],

$$ {\lambda}_0=l\pm \frac{H\alpha -H+1}{\alpha +1} $$

(10)

In (10), H denotes the number of terms in the maximum sidelobe decay window. For other windows, λ ₀ can be obtained by the polynomial approximation method [13].

Similarly, with proper combination of three or more spectral lines, a ratio α can be obtained which only depends on the selected window and λ ₀ [9, 12]. Once the data window is chosen, λ ₀ can be solely determined by α,

$$ {\lambda}_0=h\left(\alpha \right) $$

(11)

If the maximum sidelobe decay windows are selected, α can be obtained from (10). Once α is determined, the normalized frequency can be computed and the frequency can be worked out in Eq. (2).

As shown above, the second terms on the right in Eqs. (5) and (6) representing the interferences from the image component in the spectrum have been neglected. The approximation will be reasonable if λ ₀ > 5 and λ ₀ < N/2 − 5. However, if λ ₀ were out of the specified ranges, significant errors would be generated [15]. Therefore, it is of great importance and necessary to deduce a simple algorithm that is applicable even when λ ₀ is in the extreme ranges.

3 Generalization of interpolation DFT algorithms

In Section 2, it was already known that

$$ X(l)=\frac{A_0}{2}{e}^{j{\varphi}_0}{W}_N\left(-\delta \right)+\frac{A_0}{2}{e}^{-j{\varphi}_0}{W}_N\left(2l+\delta \right), $$

(12)

where W _N (·) is the DTFT of the window w _N (n) and

$$ \delta ={\lambda}_0-l $$

(13)

It should be pointed out that due to the actual processing requirements of causality, w _N (n) is a time-shifting window, where n = 0, 1 ⋅ ⋅ ⋅ N − 1. It can be obtained by

$$ {w}_N(n)=w\left(n-N/2\right) $$

(14)

w(n) is assumed to be a DFT-even window [2], which is symmetric with respect to the origin. According to the time-shifting property of DFT, we have

$$ {W}_N(k)=W(k){e}^{-jk\pi }, $$

(15)

where W(k) is the DTFT of the window w(n) and the complex exponential factor corresponds to the time shift. W(k) is also symmetric and real because w(n) is symmetric and real. Substituting Eq. (15) into Eq. (12), X(l) can be rewritten as

$$ X(l)=\frac{A_0}{2}{e}^{j\left(\delta \pi +{\varphi}_0\right)}W\left(-\delta \right)+\frac{A_0}{2}{e}^{-j\left(2l\pi +\delta \pi +{\varphi}_0\right)}W\left(2l+\delta \right) $$

(16a)

Because the period of the complex sinusoidal is 2π, Eq. (16a) can be further expressed as

$$ X(l)=\frac{A_0}{2}{e}^{j\left(\delta \pi +{\varphi}_0\right)}W\left(-\delta \right)+\frac{A_0}{2}{e}^{-j\left(\delta \pi +{\varphi}_0\right)}W\left(2l+\delta \right) $$

(16b)

The real and imaginary parts are expressed as

$$ {X}_R(l)=\frac{A_0}{2} \cos \left({\varphi}_0+\delta \pi \right)\left[W\left(-\delta \right)+W\left(2l+\delta \right)\right], $$

(17a)

and

$$ {X}_I(l)=\frac{A_0}{2} \sin \left({\varphi}_0+\delta \pi \right)\left[W\left(-\delta \right)-W\left(2l+\delta \right)\right], $$

(17b)

respectively. Accordingly, we have

$$ {X}_R\left(l\pm 1\right)=-\frac{A_0}{2} \cos \left({\varphi}_0+\delta \pi \right)\left[W\left(-\delta \pm 1\right)+W\left(2l+\delta \pm 1\right)\right], $$

(18a)

and

$$ {X}_I\left(l+1\right)=-\frac{A_0}{2} \sin \left({\varphi}_0+\delta \pi \right)\left[W\left(-\delta \pm 1\right)-W\left(2l+\delta \pm 1\right)\right] $$

(18b)

Now, we introduce two variables α _R , α _I defined as

$$ {\alpha}_R=\left|\frac{X_R\left(l\pm 1\right)}{X_R(l)}\right|=\frac{W\left(-\delta \pm 1\right)+W\left(2l+\delta \pm 1\right)}{W\left(-\delta \right)+W\left(2l+\delta \right)}, $$

(19a)

and

$$ {\alpha}_I=\left|\frac{X_I\left(l\pm 1\right)}{X_I(l)}\right|=\frac{W\left(-\delta \pm 1\right)-W\left(2l+\delta \pm 1\right)}{W\left(-\delta \right)-W\left(2l+\delta \right)} $$

(19b)

It is assumed, similar to that in Section 2, that λ ₀ is far enough from the origin so that the leakage terms W(2l + δ) and W(2l + δ ± 1) are very small compared to W(−δ) and W(−δ ± 1) and could be ignored. With this assumption, we can get the following relationship

$$ {\alpha}_R\cong {\alpha}_I\cong \alpha \cong \frac{W\left(-\delta \pm 1\right)}{W\left(-\delta \right)} $$

(20)

The relationship indicates that the ratio of the real or imaginary part can also be used for estimation if the interference from the image part can be neglected. Note that a very important phenomenon is implied from Eqs. (19a) and (19b) that the ratio is phase independent, i.e., the phase does not affect the result of frequency estimation. In other words, the frequency estimate only depends on the frequency itself.

We proceed to analyze the traditional algorithm without making approximation. The moduli of X(l) and X(l ± 1) can be obtained by the square roots of J(l) and J(l + 1), respectively, according to the results in Eqs. (17a), (17b), (18a), and (18b).

$$ J(l)=\frac{A_0^2}{4}\left[{W}^2\left(-\delta \right)+{W}^2\left(2l+\delta \right)\right]+\frac{A_0^2}{2} \cos \left(2{\varphi}_0+2\delta \pi \right)W\left(-\delta \right)W\left(2l+\delta \right) $$

(21a)

$$ J\left(l\pm 1\right)=\frac{A_0^2}{4}\left[{W}^2\left(-\delta \pm 1\right)+{W}^2\left(2l+\delta \pm 1\right)\right]+\frac{A_0^2}{2} \cos \left(2{\varphi}_0+2\delta \pi \right)W\left(-\delta \pm 1\right)W\left(2l+\delta \pm 1\right) $$

(21b)

Without making approximation, Eq. (10) becomes

$$ \overline{\alpha}=\frac{\left|X\left(l\pm 1\right)\right|}{\left|X(l)\right|}=\frac{\sqrt{J\left(l\pm 1\right)}}{\sqrt{J(l)}} $$

(22)

Obviously, $\overline{\alpha}$ is phase dependent. Also, for Hanning window, (namely two-term MSD window) it can be proved that (see Appendix 1)

$$ {\lambda}_R\ge \overline{\lambda}\ge {\lambda}_I, $$

(23)

where λ _R , λ _I , and $\overline{\lambda}$ are corresponding estimation values by α _R , α _I and $\overline{\alpha}$. If δ ≠ 0, the first equality sign is only for φ ₀ + δπ = 0 or π and the second equality sign is only for φ ₀ + δπ = ± π/2. Systematic error will be introduced no matter which of the three ratios is used in formula (11). The difference is that α _R is only affected by the real part of negative frequency and α _I is only affected by the imaginary part of the negative frequency component, while $\overline{\alpha}$ is affected by both the real and imaginary parts.

Inspired by modulus-based ratio, we can extend the ratio to a more general notation. With a proper combination of the real and imaginary parts of X(l) and X(l ± 1), we can get various kinds of ratios. For example, $\widehat{\alpha}$ can be defined as

$$ \widehat{\alpha}=\frac{\left|{X}_R\left(l\pm 1\right)\left|+\right|{X}_I\left(l\pm 1\right)\right|}{\left|{X}_R(l)\left|+\right|{X}_I(l)\right|} $$

(24)

Similar to $\overline{\alpha}$, we have

$$ \widehat{\alpha}\cong \frac{W\left(-\delta \pm 1\right)}{W_N\left(-\delta \right)}, $$

(25)

and for Hanning window, we can also obtain

$$ {\lambda}_R\ge \widehat{\lambda}\ge {\lambda}_I, $$

(26)

where $\widehat{\lambda}$ is estimation value by $\widehat{\alpha}$. It is indicated that a new estimator can be obtained if α in formula (11) is replaced by $\widehat{\alpha}$. The proof of Eq. (26) is similar to that of Eq. (23). Frequency estimation errors of the four estimators when 1 < λ ₀ < 2 (samples are weighted by Hanning window) are shown in Figs. 1, 2, 3, and 4. In Figs. 1 and 2, it can be seen that for a certain frequency, the frequency error is constant in spite of the variable phase. It is confirmed that the frequency estimators, based on the α _R , α _I separately, are independent of the signal phase, while the estimators based on $\overline{\alpha},\widehat{\alpha}$ are the functions of both the frequency and the phase. Comparing the estimated errors in Figs. 1, 2, 3, and 4, we can also see that the simulation results coincide well with the theoretical analysis in Eqs. (23) and (26). The difference between the two estimators in Figs. 3 and 4 is that the developing trend of errors in Fig. 3 is smoother than that in Fig. 4. We now have to emphasize again that maximum absolute errors are obtained when the phase and the offset satisfy the relationship φ ₀ + δπ = 0 or φ ₀ + δπ = π/2. For a certain offset, the maximum error is equal to the error shown in Figs. 1 and 2. In addition, we can also construct other types of ratio to create new estimators by combining the real and imaginary parts of X(l) and X(l ± 1) properly.

Fig. 1

Frequency estimation errors when Hanning window is used (α = α _R )

Fig. 2

Frequency estimation errors when Hanning window is used (α = α _I )

Fig. 3

Frequency estimation errors when Hanning window is used ($\alpha =\overline{\alpha}$)

Fig. 4

Frequency estimation errors when Hanning window is used ($\alpha =\widehat{\alpha}$)

4 Algorithms with high image component interference rejection

As we can see in Figs. 1, 2, 3, and 4, the maximum error due to negative frequency can reach as high as 0.04 for all the estimators. Simulation results show that it can be even up to nearly 0.1 for the rectangle window and up to 0.2 for the three-term maximum decay window. As a result, it is of great significance to reduce the interference resulting from the negative frequency. A new interpolated algorithm which has strong resistance against interference from the negative frequency component is proposed in this section.

4.1 Simple algorithms with high image component interference rejection

Recalling Eqs. (19a) and (19b) and observing Figs. 1 and 2, we can infer that there is nearly a complementary relationship between the two estimators. A proper combination of the two estimators can probably result in an intrinsic rejection of negative frequency leakage. Hence, we introduce α ₁ defined as the arithmetic mean value of α _R and α _I ,

$$ {\alpha}_1=\left({\alpha}_R+{\alpha}_I\right)/2=\frac{W\left(-\delta \pm 1\right)-W\left(2l+\delta \pm 1\right)W\left(2l+\delta \right)/W\left(-\delta \right)}{W\left(-\delta \right)-{W}^2\left(2l+\delta \right)/W\left(-\delta \right)} $$

(27a)

Now, the interference terms are W(2l + δ ± 1)W(2l + δ)/W(−δ) and W ²(2l + δ)/W(−δ), much smaller than the original terms W(2l + δ ± 1) and W(2l + δ). Most errors are eliminated by a simple arithmetic average operation. Similarly, we can also introduce α ₂, α ₃, and α ₄ defined as the geometric mean value, the harmonic mean value, and the quadratic mean of α _R and α _I , respectively, as indicated in

$$ {\alpha}_2=\sqrt{\alpha_R{\alpha}_I}=\sqrt{\frac{W^2\left(-\delta \pm 1\right)-{W}^2\left(2l+\delta \pm 1\right)}{W^2\left(-\delta \right)-{W}^2\left(2l+\delta \right)}}, $$

(27b)

$$ {\alpha}_3=\frac{2}{1/{\alpha}_R+1/{\alpha}_I}=\frac{\alpha_2^2}{\alpha_1}, $$

(27c)

and

$$ {\alpha}_4=\sqrt{\left({\alpha}_R^2+{\alpha}_I^2\right)/2}=\sqrt{2{\alpha}_1^2-{\alpha}_2^2} $$

(27d)

Similar to α ₁, the interference terms of α ₂ become much smaller as well. α ₃ and α ₄ can be written as simple functions of α ₁ and α ₂, respectively. The values of the four means are approximately equal. To demonstrate the ability of negative frequency leakage rejection, the frequency errors of the four estimators are displayed when 1 < λ ₀ < 2. Note that the four estimators are also phase independent because no phase information is involved in the ratios shown in Eqs. (27a)–(27d). Consequently, the phase was just set to zero. Figures 5 and 6 show the estimation results as a function of the frequency. It can be seen that the frequency errors sharply decrease compared with the modulus-based algorithms. In particular, the remaining error for the one based on the harmonic mean value is less than 10⁻³, which is small enough for the engineering practice.

Fig. 5

Frequency estimation errors based on different means of α _R and α _I when Hanning window is used

Fig. 6

Frequency estimation errors based on different means of α _R and α _I when three-term MSD window is used

4.2 Further improved interpolation algorithms with slide DFT

However, there is a serious defect in the above algorithms in which the weighed ratio is used. The algorithms may become quite vulnerable if cos(φ ₀ + δπ) ≈ 0 or sin(φ ₀ + δπ) ≈ 0. Under such two circumstances, the imaginary parts or the real parts would be so small that even a small disturbance would lead to a dramatic change in α _R or α _I , resulting of significant errors in the final frequency estimates. Theoretical cosine waves corrupted by low-level random noise were generated to confirm the defect, and the vulnerability of two frequency estimators based on α _R and α _I is shown in Figs. 7 and 8, respectively. It is clearly shown that radical changes appear when cos(φ ₀ + δπ) ≈ 0 for α _R and sin(φ ₀ + δπ) ≈ 0 for α _I . Sharp peaks along these lines are observed in the surface of the frequency errors. They look like knife-edge arêtes so we call them “Luo–arêtes.” In contrast, for the frequencies that are not in the vicinity of these lines, errors are very small and remain stable.

Fig. 7

The phenomenon of “Luo–arêtes” in the α _R -based estimator (Hanning window)

Fig. 8

The phenomenon of “Luo–arêtes” in the α _I -based estimator (Hanning window)

To avoid the “Luo–arêtes,” a further improved algorithm has been proposed. The fact that α _R and α _I (including various kinds of mean values of the two) are phase independent and the change of phase has no influences on frequency estimates help us get a robust ratio by time-shifting technique. Now, consider a discrete sequence with N samples x(0), x(1) ⋅ ⋅ ⋅ x(N − 1) and assume that its observed phase is φ ₁. After delaying for L samples, we get a time-shifted sequence x(L), x(L + 1) ⋅ ⋅ ⋅ x(N + L − 1) and its observed phase φ ₂ can be expressed as

$$ {\varphi}_2={\varphi}_1+2\pi {f}_0L\varDelta T={\varphi}_1+2\pi {f}_0L/{f}_s={\varphi}_1+2\pi {\lambda}_0L/N $$

(28)

In the above equation, ΔT is the sampling interval, f ₀ is the theoretical frequency and λ ₀ is the normalized frequency scaled by frequency resolution. As λ ₀ is unknown, we can use its largest bin number l, instead.

$$ L\approx N\frac{\varphi_2-{\varphi}_1}{2\pi l} $$

(29)

The actual observed phase of the time delay sequence is φ ' ₂ = φ ₁ + 2πlL/N and the phase error is

$$ {\varphi}_2\hbox{'}-{\varphi}_2=\varDelta {\varphi}_2=2\pi \delta L/N $$

(30)

Generally, we have l > 1 and φ ₂ − φ ₁ < π/2 so that L/N < 1/4. Considering δ ∈ [−0.5, 0.5), the absolute phase error is less than π/4. For α _R , if φ ₂ was set to 0 or π, φ ' ₂ was limited in the range of (−π/4, π/4) or (3π/4, 5π/4). For α _I , if φ ₂ was set to π/2 or − π/2, φ ' ₂ was limited in the range of (π/4, 3π/4) or (−3π/4, − π/4). It is indicated that we can adjust the observed phase by the time-shifting technique so that we can get a robust ratio to avoid the “Luo–arêtes.” In addition, we can obtain the spectral lines l and l ± 1 of the time-delayed sequence by means of the sliding discrete Fourier transform (SDFT) [19, 20]. It will be more efficient compared with another separate FFT or DFT.

5 Comparison with other state-of-the-art methods

In this section, some computer simulations were conducted to verify the effectiveness and the accuracy of the proposed algorithms. For conciseness, only the results of the algorithm based on the harmonic mean of α _R and α _I were shown. All simulation results are returned by the algorithm proposed in Section 4.2. In addition, the results of the traditional IpDFT algorithm (2pIpDFT) [5, 8, 11], the classical three-point IpDFT algorithm (3pIpDFT) [9, 12], Quinn’s two-point-based estimator (Quinn1, only applicable for Hanning window), and Quinn’s three-point-based optimal estimator (Quinn2, only applicable for Hanning window) [21–23], the estimator proposed by Macleod in 1998 (Macleod, only applicable for Hanning window) [24], the three-point complex spectrum-based estimator (Jacobsen, only applicable for Hanning window) [25] referred by Jacobsen and Kootsookos in 2007, and the improved three-point IpDFT algorithm with high image frequency interference rejection capability recently proposed by Belega (Belega14) [15], were also displayed for comparison. In all the considered IpDFT-based algorithms, both the two-term and three-term MSD windows were adopted. For simplicity but without loss of generality, parameters used in all the simulation experiments were as follows, the amplitude of the cosine wave A = 1, the number of samples N = 512, and the sampling rate f _s  = 512. The results shown in this section were scaled by frequency resolution and expressed in bins.

5.1 Theoretical cosine wave without noise

To demonstrate the excellent rejection capability against the interference from negative frequency, theoretical signals contain a small number of cycles. The normalized frequency λ ₀ is varied in the range (1, 11) with a step of 1/8. For each frequency, the phase θ is varied in the range [−π, π) with a step of π/72. The maximum absolute frequency errors |δ|_max are shown in Figs. 9 and 10 as a function of λ ₀ for the two-term and three-term MSD windows, respectively. It is clearly revealed that the errors due to negative frequency interference are remarkably reduced. In general, both the proposed method and the improved three-point IpDFT method outperform other estimators throughout the entire range of considered λ ₀. When the Hanning window is adopted, Jacobsen’s estimator has the worst performance. The traditional IpDFT algorithm, Quinn’s two-point-based estimator, and Quinn’s three-point-based optimal estimator have similar trend. The classical three-point IpDFT algorithm and Macleod’s estimator provide better results than the above three estimators. If λ ₀ < 4.5 (especially λ ₀ < 1.5) and δ > 0, the new method provides better performance than the improved three-point IpDFT and the opposite holds if δ < 0. When λ ₀ > 4.5, the improved three-point IpDFT shows a small advantage than the proposed algorithm. When the three-term MSD window is adopted, the new estimator has overwhelming advantage over the rest. It should be pointed out that regardless of the adopted window type and the value of λ ₀, the maximum frequency error of the new algorithm never goes above 10⁻³, even though only one or two cycles are obtained.

Fig. 9

Maximum estimation error as a function of λ ₀ for different estimators without noise (Hanning window)

Fig. 10

Maximum estimation error as a function of λ ₀ for different estimators without noise (three-term MSD window)

5.2 Theoretical cosine wave corrupted by additive noise

In this subsection, we considered the ideal cosine wave contaminated with the additive Gaussian noise. Similar to [15], we investigated the RMSE of estimates returned by the considered estimators as a function of λ ₀ for certain SNRs. For each frequency, 50,000 instances were generated with a random phase. For conciseness, four algorithms are considered, including the traditional IpDFT algorithm (2pIpDFT), the classical three-point IpDFT algorithm (3pIpDFT), the improved three-point IpDFT algorithm with high image frequency interference rejection capability (Belega14), and the proposed algorithm. Results with SNR = 30 dB and 50 dB for Hanning window and three-term MSD window were shown in Figs. 11, 12, 13, and 14, respectively. The results for Hanning window agree well with those in [15]. As shown in Figs. 11 and 12, the estimated RMSE of the proposed method is always in a low level with a small fluctuation and essentially the novel method has no competitor for λ ₀ < 3. When λ ₀ becomes larger, the estimated RMSE of four estimators tend to be in a similar level and it is interesting to find that the two-point-based algorithms show a better performance at the worst incoherent sampling condition (δ ≈ 0.5), while the three-point-based algorithms provide better results when λ ₀ is close to integer values (δ ≈ 0).

Fig. 11

RMSE for different algorithms (SNR = 30 dB, Hanning window)

Fig. 12

RMSE for different algorithms (SNR = 50 dB, Hanning window)

Fig. 13

RMSE for different algorithms (SNR = 30 dB, three-term MSD window)

Fig. 14

RMSE for different algorithms (SNR = 50 dB, three-term MSD window)

When the three-term MSD window is adopted, the overall trend is similar to that of Hanning window. For λ ₀ < 2, the performance of the traditional 2pIpDFT and 3pIpDFT is worse than Hanning window due to the wider mainlobe of the three-term MSD window. For λ ₀ > 2, they have better performances because of the faster sidelobe decay rate. For the same SNR, the RMSE values of the two algorithms employing three-term MSD window decrease faster than that employing Hanning window before reaching the ultimate stable level. Meanwhile, the ultimate stable level is higher than that employing Hanning window because of the worse noise properties of the three-term MSD window. The performance of the proposed algorithm maintains its superiority to the other three algorithms for λ ₀ < 3. To sum up, the traditional two-point- or three-point-based algorithms are very good choices because of their simplicity when random noise has a significant influence on the uncertainty of estimates. The proposed algorithm is strongly recommended when the image frequency interference is the main error source especially when a small number of cycles are contained in samples.

6 Conclusions

Frequency estimation by the IpDFT method is studied in this paper. We have quantitatively analyzed the influences of the interference from the image component and generalized the interpolated DFT algorithms. Based on the analysis, novel frequency estimators are proposed, which have strong rejection against the high image component interference. Accuracy of the novel algorithm has been confirmed by simulations. Comparative studies reveal that the proposed algorithm has a better performance than the traditional algorithms when the spectral interference from negative frequency component is significant, especially for the cycles less than one and a half. This proposed algorithm is simple to understand, easy to implement, and very suitable for real-time analysis.

Appendices

Appendices

1.1 Appendix 1

According to the conclusion in Eq. (10), we can get the frequency estimation by

$$ \lambda =l\pm \frac{2\alpha -1}{\alpha +1}, $$

(A1)

when Hanning window is used.

Accordingly, we can obtain frequency estimations

$$ {\lambda}_R=l\pm \frac{2{\alpha}_R-1}{\alpha_R+1}=l\pm \left(2-\frac{3}{\alpha_R+1}\right), $$

(A2a)

$$ {\lambda}_I=l\pm \frac{2{\alpha}_I-1}{\alpha_I+1}=l\pm \left(2-\frac{3}{\alpha_I+1}\right), $$

(A2b)

and

$$ \overline{\lambda}=l\pm \frac{2\overline{\alpha}-1}{\overline{\alpha}+1}=l\pm \left(2-\frac{3}{\overline{\alpha}+1}\right), $$

(A2c)

by three kinds of ratio α _R , α _I , and $\overline{\alpha}$, respectively.

Further, combining Eq. (A2a) with Eq. (A2c) and Eq. (A2b) with Eq. (A2c), we get the following two equations

$$ {\widehat{\lambda}}_R-\overline{\lambda}=\pm \left(\frac{3}{\overline{\alpha}+1}-\frac{3}{\alpha_R+1}\right)=\pm 3\frac{\alpha_R-\overline{\alpha}}{\left(\overline{\alpha}+1\right)\left({\alpha}_R+1\right)} $$

(A3a)

and

$$ {\widehat{\lambda}}_I-\overline{\lambda}=\pm \left(\frac{3}{\overline{\alpha}+1}-\frac{3}{\alpha_I+1}\right)=\pm 3\frac{\alpha_I-\overline{\alpha}}{\left(\overline{\alpha}+1\right)\left({\alpha}_I+1\right)} $$

(A3b)

The above equations indicate the following conclusions.

• (1) For δ > 0, since ${\alpha}_R>\overline{\alpha}$, ${\alpha}_I<\overline{\alpha}$ (Appendix 2) and it takes the positive in (A3a) and (A3b), we get ${\widehat{\lambda}}_R>\overline{\lambda}$ and ${\widehat{\lambda}}_I<\overline{\lambda}$.

• (2) For δ < 0, since ${\alpha}_R<\overline{\alpha}$, ${\alpha}_I>\overline{\alpha}$ (Appendix 2) and it takes the negative in (A3a) and (A3b),we get ${\widehat{\lambda}}_R>\overline{\lambda}$ and ${\widehat{\lambda}}_I<\overline{\lambda}$.

Finally, for both δ > 0 and δ < 0, we have ${\widehat{\lambda}}_R>\overline{\lambda}$ and ${\widehat{\lambda}}_I<\overline{\lambda}$. Furthermore, if δ ≠ 0, when φ ₀ + δπ = 0 or π, we have ${\widehat{\lambda}}_R=\overline{\lambda}$; and when φ ₀ + δπ = ± π/2, we have ${\widehat{\lambda}}_I=\overline{\lambda}$.

1.2 Appendix 2

It is already known that ${\alpha}_R=\frac{W_1+{W}_{2l+1}}{W_0+{W}_{2l}}$, ${\alpha}_I=\frac{W_1-{W}_{2l+1}}{W_0-{W}_{2l}}$, $\overline{\alpha}=\sqrt{\frac{W_1^2+{W}_{2l+1}^2+2\mu {W}_1{W}_{2l+1}}{W_0^2+{W}_{2l}^2+2\mu {W}_0{W}_{2l}}}$ from Section 3, where W ₀ = W(−δ), W ₁ = W(−δ ± 1), W _2l  = W(2l + δ), W _2l + 1 = W(2l + δ ± 1) and u = cos(2φ ₀ + 2δπ) according to Eqs. (19a), (19b), and (22). We define

$$ E={\alpha}_R^2-{\overline{\alpha}}^2=\frac{W_1^2+{W}_{2l+1}^2+2{W}_1{W}_{2l+1}}{W_0^2+{W}_{2l}^2+2{W}_0{W}_{2l}}-\frac{W_1^2+{W}_{2l+1}^2+2u{W}_1{W}_{2l+1}}{W_0^2+{W}_{2l}^2+2u{W}_0{W}_{2l}}, $$

(B1a)

and

$$ F={\alpha}_I^2-{\overline{\alpha}}^2=\frac{W_1^2+{W}_{2l+1}^2-2{W}_1{W}_{2l+1}}{W_0^2+{W}_{2l}^2-2{W}_0{W}_{2l}}-\frac{W_1^2+{W}_{2l+1}^2+2u{W}_1{W}_{2l+1}}{W_0^2+{W}_{2l}^2+2u{W}_0{W}_{2l}} $$

(B1b)

After some algebraic manipulation, we obtain

$$ E=2\left(1-u\right)\frac{G}{{\left({W}_0+{W}_{2l}\right)}^2\left({W}_0^2+{W}_{2l}^2+2u{W}_0\kern0.5em {W}_{2l}\right)}, $$

(B2a)

and

$$ F=2\left(1+u\right)\frac{-G}{{\left({W}_0-{W}_{2l}\right)}^2\left({W}_0^2+{W}_{2l}^2+2u{W}_0\kern0.5em {W}_{2l}\right)}, $$

(B2b)

where G = (W ₀  W ₁ − W _2l W _2l + 1)(W ₀  W _2l + 1 − W ₁  W _2l ). When Hanning window is used, we know that W ₀ > W ₁ > 0, W ₁ > > |W _2l |, and W ₁ > > |W _2l + 1|.

• (1) For δ > 0, we have W _2l  < 0 and W _2l + 1 > 0 for Hanning window, so G > 0. Then, we have ${\alpha}_R>\overline{\alpha}$ and ${\alpha}_I<\overline{\alpha}$.

• (2) For δ < 0, we have W _2l  > 0 and W _2l + 1 < 0 for Hanning window, so G < 0. Then, we have ${\alpha}_R<\overline{\alpha}$ and ${\alpha}_I>\overline{\alpha}$.

It should be pointed that for δ ≠ 0, the requirement for ${\alpha}_R=\overline{\alpha}$ is that μ = 1 and the requirement for ${\alpha}_I=\overline{\alpha}$ is that μ = − 1. That means that φ ₀ + δπ = 0 or π for ${\alpha}_R=\overline{\alpha}$ and φ ₀ + δπ = ± π/2 for ${\alpha}_I=\overline{\alpha}$.