Frequency estimators with high resistance to interference from image frequency components

Abstract

This paper presents an improved frequency estimation algorithm based on the interpolated discrete Fourier transform. High-accurate frequency estimation can be achieved by taking the geometric mean of two independent estimates, which are derived from the real parts of the two largest spectral bins and the imaginary parts, respectively. In situations where only a small number of sine wave cycles are observed, the ability of the algorithm to cancel interference from image frequency components results in improvements in accuracy. The residual errors of the proposed algorithm have been theoretically analyzed with maximum side-lobe decaying windows, since the windows have simple and uniform analytical expression of interpolation algorithm. The performance of the proposed algorithm was investigated using both Hanning and three-term maximum side-lobe decaying windows. A comparative analysis of systematic errors and noise sensitivity was performed between the new algorithm and traditional algorithms. Both the root mean squared error and the probability density of the errors were investigated under noisy conditions. Simulation results demonstrated that the new algorithm is not only highly resistant to interference from image components but is also resistant to the effects of random noise. The results presented in the paper are useful for identifying the best choice of algorithm in practical engineering applications.

Introduction

Spectral analysis in signal processing has been widely used in various engineering applications such as bio-signal processing, vibration signal analysis, power harmonic analysis, and radar ranging.^1–4 Spectrum domain-based parameter estimation methods have received great interest in recent decades because they are simple to implement and highly efficient.^1–4 However, there are two significant drawbacks to discrete Fourier transform (DFT) approaches that can lead to significant estimation errors: the spectral leakage effect, arising from the finite extent of the signal, and the picket fence effect, caused by discretization of the frequency domain.^5,6 The problems can be handled using the windowing technique and frequency interpolation, respectively.^3,4,7,8 Windowed interpolation discrete Fourier transform (IpDFT) algorithms have been shown to be effective in compensating for errors,^9–11 and these algorithms have been extended from two DFT spectral bins to three or more,^12,13 and to use complex DFT samples.^14–19 In 1983, the Hanning window is used by Grandke, ¹⁰ obtained the more accurate results for interpolated frequencies. In 2009, the IpDFT method with maximum side-lobe decaying (MSD) windows is studied and acquired the analytical formulas for parameters estimation of multi-frequency signals by Belega and Dallet. ¹¹

Nevertheless, attention still needs to be given to interference from image components, which has been mostly ignored by existing IpDFT algorithms.^2,3,20,21 This is a particular problem in many practical engineering applications where both speed and accuracy are critical. ²² For example, in a photovoltaic system, accurate frequency detection needs to be achieved in a short measurement time, on the order of one or two periods. However, if a sampled sine wave contains only a few cycles, the image frequency component can cause considerable interference, resulting in significant errors in DFT-based methods. ²³ For this reason, several methods have been proposed in recent years that exhibit excellent resistance to image component interference.^23–26 One of the effective algorithms was proposed by Chen et al.,²⁴ and Borkowski et al. ²⁵ further developed this approach, extending the range of application to maximum side-lobe decaying windows. In 2014, Belega et al. ²³ proposed an improved three-point interpolation algorithm that performs well in resisting image frequency component disturbances. Recently, a novel algorithm was developed by Luo et al. ²⁶ based on two-point interpolated DFT. Even though Luo’s algorithm has obvious advantages against random noise, its performance in noiseless conditions can be further improved.

This paper aims at studying how to further reduce the image component interference based on Luo’s idea in an attempt to achieve more accurate results. The rest of the paper is organized as follows. In “Theoretical foundation” section, a brief introduction of the classic interpolation DFT methods based on MSD window was presented. In “Proposed interpolation estimator” section, Luo’s algorithm was restated, and the proposed algorithm was derived. In “Performance of the proposed algorithm” section, some computer simulations were conducted to compare the performances of different algorithms, and the probable reasons of the phenomenon in simulations were also discussed. Finally, some concluding remarks were made.

Theoretical foundation

An analysis of the general theory of windowed interpolation algorithms begins by considering a discrete signal of the form

$$x(n)=A_0\cos (2\pi \frac{{\lambda }_0}{N}n+{\varphi }_0),n=0,1,2,\cdots ,N-1$$

(1)where $N$ is the time samples and $A_0$ , ${\lambda }_0$ , and ${\varphi }_0$ are the amplitude, the normalized frequency, and phase, respectively. When the signal is asynchronously sampled, ${\lambda }_0$ lies between two DFT bins. Therefore, ${\lambda }_0$ can be expressed as

$${\lambda }_0=l_0+\delta$$

(2)where $l_0$ and $\delta$ represent the integer and the fractional parts,¹³ respectively. Supposing that the largest magnitude bin in the spectrum is located at $l$ , the largest magnitude is given by ²⁶ :

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

(3)and the second largest magnitude is

$$\begin{array}{c}\left|X(l\pm 1)\right|=|\frac{A_0}{2}{e}^{j(\delta \pi +{\varphi }_0)}W(-\delta \pm 1)\\ +\frac{A_0}{2}{e}^{-j(\delta \pi +{\varphi }_0)}W(2l+\delta \pm 1)|\end{array}$$

(4)where the negative and positive signs correspond to $l>{\lambda }_0$ and $l<{\lambda }_0$ , respectively. In equations (3) and (4), $W(\cdot )$ is the discrete-time Fourier transform of the window. Usually, the integer part $l_0$ of the normalized frequency can be determined using a coarse search for the largest magnitude in the discrete spectrum, while the fractional part $\delta$ can be estimated using two or more adjacent DFT samples. ¹² If $l$ is far enough from 0 and $N/2$ , the second terms on the right in equations (3) and (4) are both very small compared with the first term, and consequently they can be ignored. In this situation, equation (3) becomes

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

(5)

Similarly, equation (4) can be rewritten as

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

(6)

In two-point interpolation algorithms, the ratio $\alpha$ is defined by

$$\alpha =\frac{|X(l\pm 1)|}{|X(l)|}\cong {\alpha }_0=\frac{W(-\delta \pm 1)}{W(-\delta )}$$

(7)

For the sake of simplicity, we focus on algorithms based on MSD windows because there is a simple and explicit relationship between $\delta$ and $\alpha$ given by

$$\delta =f(\alpha ,H)=\pm \frac{H\alpha -(H-1)}{\alpha +1}$$

(8)where the positive sign and negative sign correspond to, respectively, $\left|X(l+1)\right|>\left|X(l-1)\right|$ and $\left|X(l+1)\right|<\left|X(l-1)\right|$ _. ¹³ The parameter H is the term number of the MSD window.

Proposed interpolation estimator

It is worth noting that the frequency deviation estimate $\delta$ in equation (8) is obtained without accounting for interference from the image frequency component. As the normalized frequency ${\lambda }_0$ decreases, estimation errors will increase due to the terms $\left|W(2l+\delta )\right|$ and $\left|W(2l+\delta \pm 1)\right|$ . To reduce this detrimental effect, Luo et al. suggested a modification to the ratio $\alpha$ . ²⁶

Luo et al. suggested that the real and imaginary parts of the largest spectral samples can be expressed as ²⁶

$$X_R(l)=\frac{A_0}{2}\cos ({\phi }_0+\delta \pi )[W(-\delta )+W(2l+\delta )]$$

(9a)and

$$X_I(l)=\frac{A_0}{2}\sin ({\phi }_0+\delta \pi )[W(-\delta )-W(2l+\delta )]$$

(9b)

Similarly, the real and imaginary parts of the second largest spectral line are given by ²⁶

$$X_R(l\pm 1)=\frac{A_0}{2}\cos ({\phi }_0+\delta \pi )[W(-\delta \pm 1)+W(2l+\delta \pm 1)]$$

(10a)and

$$\begin{array}{c}{X}_I(l\pm 1)=\frac{A_0}{2}\sin ({\phi }_0+\delta \pi )[W(-\delta \pm 1)-W(2l+\delta \pm 1)]\end{array}$$

(10b)

If $\cos ({\phi }_0+\delta \pi )\ne 0$ and $\sin ({\phi }_0+\delta \pi )\ne 0$ , two variables ${\alpha }_R$ and ${\alpha }_I$ are introduced

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

(11a)and

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

(11b)

In Luo’s algorithm, the modified ratio $\widehat{\alpha }$ is the harmonic mean of ${\alpha }_R$ and ${\alpha }_I$

$$\widehat{\alpha }=\frac{2}{1/{\alpha }_R+1/{\alpha }_I}$$

(12)

Furthermore, equations (11a), (11b), and (12) can be combined to give

$$\widehat{\alpha }=\frac{W^2(-\delta \pm 1)-W^2(2l+\delta \pm 1)}{W(-\delta \pm 1)W(-\delta )-W(2l+\delta \pm 1)W(2l+\delta )}$$

(13)

We can see that the interference terms are $W^2(2l+\delta \pm 1)$ and $W(2l+\delta \pm 1)W(2l+\delta )$ , and that they are much smaller than the terms $W^2(-\delta \pm 1)$ and $W(-\delta \pm 1)W(-\delta )$ . This shows that the errors in equations (7) are largely eliminated by the harmonic mean operation. Finally, the frequency estimator is given by

$$\widehat{\lambda }=l+\widehat{\delta }=l+f(H,\widehat{\alpha })$$

(14)

However, in the simulation results of Luo et al., ²⁶ estimation errors have still not been completely removed, especially in case where the number of cycles is small. Inspired by Luo’s algorithm, we propose an improved algorithm that further reduces the remaining estimation errors. In contrast to Luo’s algorithm, in which the final frequency estimate was achieved by replacing the $\alpha$ in equation (8) with the $\widehat{\alpha }$ of equation (14), the frequency estimate in the new algorithm is computed by the geometric mean of two frequency estimates, respectively, produced from ${\alpha }_R$ and ${\alpha }_I$ . The two estimates are denoted by ${\lambda }_R$ and ${\lambda }_I$

$${\lambda }_R=l+{\delta }_R=l+f(H,{\alpha }_R)$$

(15a)and

$${\lambda }_I=l+{\delta }_I=l+f(H,{\alpha }_I)$$

(15b)

Supposing $\left|X(l+1)\right|>\left|X(l-1)\right|$ , and substituting equation (8) into equations (15a) and (15b), gives

$${\lambda }_R=l+\frac{H{\alpha }_R-(H-1)}{{\alpha }_R+1}=l+H-\frac{2H-1}{{\alpha }_R+1}$$

(16a)

$${\lambda }_I=l+\frac{H{\alpha }_I-(H-1)}{{\alpha }_I+1}=l+H-\frac{2H-1}{{\alpha }_I+1}$$

(16b)

Similarly, when $\left|X(l+1)\right|<\left|X(l-1)\right|$ , we have

$${\lambda }_R=l-\frac{H{\alpha }_R-(H-1)}{{\alpha }_R+1}=l-H+\frac{2H-1}{{\alpha }_R+1}$$

(17a)

$${\lambda }_I=l-\frac{H{\alpha }_I-(H-1)}{{\alpha }_I+1}=l-H+\frac{2H-1}{{\alpha }_I+1}$$

(17b)

As a result, the product of ${\lambda }_R$ and ${\lambda }_I$ can be achieved by

$$\begin{array}{l}{\lambda }_R{\lambda }_I={(l\pm H)}^2\mp (l\pm H)(2H-1)(\frac{1}{{\alpha }_R+1}+\frac{1}{{\alpha }_I+1})\\ +{(2H-1)}^2\frac{1}{({\alpha }_R+1)({\alpha }_I+1)}\end{array}$$

(18)where

$$\begin{array}{l}\frac{1}{{\alpha }_R+1}+\frac{1}{{\alpha }_I+1}\\ =2\frac{\left\{\begin{array}{l}{W}^2(-\delta )+W(-\delta )W(-\delta \pm 1)\\ -W^2(2l+\delta )-W(2l+\delta )W(2l+\delta \pm 1)\end{array}\right\}}{\left\{\begin{array}{l}{[W(-\delta )+W(-\delta \pm 1)]}^2\\ -{[W(2l+\delta )+W(2l+\delta \pm 1)]}^2\end{array}\right\}}\end{array}$$

(19a)and

$$\frac{1}{({\alpha }_R+1)({\alpha }_I+1)}=\frac{W^2(-\delta )-W^2(2l+\delta )}{\left\{\begin{array}{l}{[W(-\delta )+W(-\delta \pm 1)]}^2\\ -{[W(2l+\delta )+W(2l+\delta \pm 1)]}^2\end{array}\right\}}$$

(19b)

The frequency estimate is expressed as

$$\widehat{\lambda }=\sqrt{{\lambda }_R{\lambda }_I}$$

(20)

The residual error is

$$\Delta \lambda =\sqrt{{\lambda }_R{\lambda }_I}-{\lambda }_0=\frac{{\lambda }_R{\lambda }_I-{\lambda }_0^2}{\sqrt{{\lambda }_R{\lambda }_I}+{\lambda }_0}=\frac{{\hspace{0.17em}}_{\lambda }^e}{\sqrt{{\lambda }_R{\lambda }_I}+{\lambda }_0}\cong 0$$

(21)where ${\lambda }_0=l\pm H\mp \frac{2H-1}{{\alpha }_0+1}$ and ${\lambda }_e={\lambda }_R{\lambda }_I-{\lambda }_0^2\cong 0$ (see Appendix 1). Equation (21) demonstrates that the error caused by image interference can be completely removed by taking the geometric mean of ${\lambda }_R$ and ${\lambda }_I$ . It should also be emphasized that when $\cos ({\phi }_0+\delta \pi )\ne 0$ or $\sin ({\phi }_0+\delta \pi )\ne 0$ , the values of ${\alpha }_R$ or ${\alpha }_I$ are small and easily disturbed by the noise. In these cases, a so-called $Luo-ar\widehat{e}tes$ phenomena emerge due to the significant errors in the proposed algorithm. ²⁶ In order to avoid the defect and obtain a more robust estimate, a time-shifting technique is recommended. ²⁶ Finally, we obtain

$$\widehat{\lambda }=\left\{\begin{array}{cc}\sqrt{{\lambda }_R{\overline{\lambda }}_I},& \left|\cos [\angle X(l)]\right|\ge \sqrt{2}/2\\ \sqrt{{\overline{\lambda }}_R{\lambda }_I},& \left|\cos [\angle X(l)]\right|<\sqrt{2}/2\end{array}\right.$$

(22)where $\angle X(l)$ denotes the angle of $X(l)$ , the ${\overline{\lambda }}_R$ , or ${\overline{\lambda }}_I$ is calculated from the time-shifted sequence in the same way as ${\lambda }_R$ or ${\lambda }_I$ , respectively.

Performance of the proposed algorithm

This section describes computer simulations performed to make a comparative study between the new algorithm and the following new algorithms: the traditional two-point IpDFT algorithm (Grandke83,¹⁰ Belega09¹¹), the traditional three-point IpDFT algorithm ¹² (Agrez02), the improved three-point IpDFT algorithm recently proposed by Belega²³) (Belega14), the complex spectrum-based algorithms^24,25 (Chen09 and Borkowski14), and the latest two-point algorithm proposed by Luo ²⁶ (Luo16). To make the comparisons simple but without losing generality, all simulations were done with identical parameters: the cosine wave amplitude was A = 1, the number of analyzed samples was N = 256, and the sampling rate was $f_s$ =256 (so that the frequency resolution was 1). In addition, all frequency estimation results are presented in bins.

Analysis of residual errors

In order to clearly demonstrate the advantages of resistance to image component interference, signals were generated with only a few cycles. The normalized frequency ${\lambda }_0$ was varied within the range (1,13) in steps of $1/8$ . For each frequency, the phase $\theta$ scanned the range of $(-\pi ,\pi ]$ in steps of $\pi /72$ . The maximum absolute frequency errors were considered.

The results of the simulations for each algorithm are shown in Figure 1 for both the Hanning window and three-term MSD windows. In Figure 1(a), the estimation errors of Grandke83 and Agrez02 gradually decrease with increasing ${\lambda }_0$ . Similar phenomena also appear in the Belega14 and Luo16 algorithms, but the decreases are more rapid. This implies that image component interference from the negative frequency was more effectively eliminated compared with Grandke83 and Agrez02. It is also shown that the errors of both Chen09 and the proposed algorithm are stable, with oscillatory behavior on the whole range of frequencies, and that error levels are less than order 10⁻⁸ for as few as one or two cycles. The two algorithms are clearly more accurate than the other four. In particular, the proposed algorithm gives excellent results, with a maximum estimation error of only around ${10}^{-9}$ . This demonstrates that Chen09 and the proposed algorithm are strongly resistant to image component interference. For ${\lambda }_0>10$ , the four algorithms (Belega14, Chen09, Luo16, and the proposed algorithm) have very similar performance.

Figure 1.

Maximum estimation error as a function of ${\lambda }_0$ . (a) Hanning window and (b) three-term MSD window.

In Figure 1(b), the results are mostly similar to those in Figure 1(a), except for the following differences. First, the error curves obtained from Belega09, Agrez02, Belega14, and Luo16 decrease more sharply. This is mainly because the three-term MSD window has more rapid side-lobe attenuation than the Hanning window. And second, Borkowski14 and the proposed algorithm exhibit lower errors as well as smaller fluctuations, and the estimation errors for the proposed algorithm are almost independent of ${\lambda }_0$ . The simulation results confirm that the theoretical analysis of the proposed algorithm is robust. In noiseless conditions, Borkowski14 or the proposed algorithm is recommended because of their excellent performance against image component interference.

The dependence of root mean squared error on the signal-to-noise ratio

In order to assess the sensitivity of the new algorithm to random noise, the behavior of the root mean squared error (RMSE) of frequency estimates was analyzed as a function of the signal-to-noise ratio (SNR). For each SNR, 10,000 test signals were generated with the normalized frequency changed randomly in the range [1, 5] and phase varied randomly in the range $(-\pi ,\pi ]$ .

Figure 2(a) shows the results with the Hanning window. The RMSEs of Grandk83, Agrez02, Belega14, and Luo16 all decrease with increasing SNR over a certain range and then are almost independent of SNR. However, Chen09 and proposed algorithm exhibit a continuously decreasing trend over the whole range, and the RMSEs are parallel to the Crame-Rao lower bound (CRLB) over a small range. At low SNR, noise has a greater effect on estimation errors than systematic errors, but as the SNR increases, systematic errors gradually come to dominate. As a result, the RMSEs of the four algorithms (Grandk83, Agrez02, Belega14, and Luo16) become flat once the SNR exceeds a certain value. For Chen09 and the proposed algorithm, since the remaining errors are much smaller than the errors caused by noise, the RMSEs decrease over the whole range of SNR under study. Luo16 and the proposed algorithm exhibit almost the same performance and possess a smaller RMSE than other four algorithms when SNR < 40 dB, and the proposed algorithm has no competitor when SNR > 40 dB. A similar phenomenon also appears in Figure 2(b) with three-term MSD window except that the RMSE level for each algorithm is higher than the corresponding ones in Figure 2(a). This is mainly due to the wider equivalent noise bandwidth (ENBW), which is an index representing the noise properties of the window. The simulation results indicate that the proposed algorithm is more robust than the other algorithms in noisy conditions.

''

Figure 2.

RMSE of different algorithms. (a) Hanning window and (b) three-term MSD window. RMSE: root mean squared error; SNR: signal-to-noise ratio.

Probability density of the errors

To further examine the noise properties of the new algorithm, the probability density of the error (PDE) was studied as a function of normalized frequency. The SNR was set to 0 dB and 20 dB to represent the high- and low-noise levels, respectively. The normalized frequency was scanned from 1 to 11 in steps of 0.1. For each frequency, 300,000 test signals were generated with a random phase. The simulation results for different algorithms at SNR = 20 dB are shown as in Figures 3 to 8. It is interesting to find that each figure shows a saw-tooth boundary concentrated in a small interval around zero. The boundary is convex when ${\lambda }_0$ is close to integer values ( $\delta \approx 0$ ) and concave in the most incoherent sampling condition ( $\delta \approx 0.5$ ) for the two-point-based algorithms. Opposite behavior can be observed in the three-point-based algorithms. The phenomenon can be partly explained by the fact that for the two-point-based algorithms, when ${\lambda }_0$ is close to integer values, the second largest spectral lines are relatively small and can be easily influenced by random noise. In contrast, under the most incoherent sampling condition, since the two largest spectral lines are relatively strong, they have greater immunity to noise. The opposite is true for the three-point-based algorithms, which usually have better performance when $\delta \approx 0$ and are more sensitive to noise when $\delta \approx 0.5$ . Generally, the two-point IpDFT algorithms perform better than the three-point IpDFT algorithms. In addition, it can be observed that the error distribution for each algorithm is more concentrated in the Hanning window compared to the three-term MSD window. This is mainly because the three-term MSD window has a wider main-lobe, resulting in a wider ENBW, and therefore, less immunity to noise (see “The dependence of root mean squared error on the signal-to-noise ratio” section).

Figure 3.

PDE for Granke83 and Belega09 with SNR = 20 dB. (a) Hanning window and (b) three-term MSD window.

Figure 4.

PDE for Agrez02 with SNR = 20 dB. (a) Hanning window and (b) three-term MSD window.

Figure 5.

PDE for Belega14 with SNR = 20 dB. (a) Hanning window and (b) three-term MSD window.

Figure 6.

PDE for Chen09 and Borkowski14 with SNR = 20 dB. (a) Hanning window and (b) three-term MSD window.

Figure 7.

PDE for Luo16 with SNR = 20 dB. (a) Hanning window and (b) three-term MSD window.

Figure 8.

PDE for Proposed with SNR = 20 dB. (a) Hanning window and (b) proposed three-term MSD window.

Simulation results at SNR = 0 dB are presented in Figures 9 to 14. All algorithms have wider spreads, and the saw-tooth boundary is more clear than with SNR = 20 dB. The proposed algorithm and Luo16 have nearly the same performance, and these two algorithms have an overwhelming advantage over the other four in noisy conditions because the error band is narrower and the error distribution is more concentrated. In summary, when the signals are subjected to random noise, the two-point-based IpDFT algorithms have better performance than the three-point-based IpDFT algorithms. For each algorithm, the performance with the Hanning window is more robust than that with the three-term MSD window. If both image frequency interference and random noise are significantly involved, the proposed algorithm is strongly recommended on account of its robustness against random noise and excellent resistance to image frequency interference.

Figure 9.

PDE for Granke83 and Belega09 with SNR = 0 dB. (a) Hanning window and (b) three-term MSD window.

Figure 10.

PDE for Agrez02 with SNR = 20 dB. (a) Hanning window and (b) Agrez02 three-term MSD window.

Figure 11.

PDE for Belega14 with SNR = 0 dB. (a) Hanning window and (b) three-term MSD window.

Figure 12.

PDE for Chen09 and Borkowski14 with SNR = 0 dB. (a) Hanning window and (b) three-term MSD window.

Figure 13.

PDE for Luo16 with SNR = 0 dB. (a) Hanning window and (b) three-term MSD window.

Figure 14.

PDE for proposed algorithm with SNR = 0 dB. (a) Hanning window and (b) three-term MSD window.

Conclusion

In this paper, an improved IpDFT algorithm is proposed that can accurately estimate frequencies when only a few signal cycles are observed. The new estimator is an extension of Luo’s algorithm that exhibits good resistance to image component interference. Systematic errors under noiseless conditions were investigated, and the PDE and the RMSEs for different algorithms were compared under conditions where the signals were disturbed by random additive noise. The theoretical analysis demonstrates that the error caused by image interference can be completely reduced. The simulations confirm that the proposed algorithm not only exhibits strong resistance against image component interference but also demonstrates good immunity to noise. Overall, the results indicate that the Borkowski14 algorithm and the proposed algorithm are acceptable choices in noiseless conditions, and Luo16 and the proposed algorithm are good choices in situations with low SNR. When there is both random noise and significant image frequency interference, the proposed algorithm becomes the best choice.