Next Article in Journal
Autonomous Multi-Rotor Aerial Platform for Air Pollution Monitoring
Previous Article in Journal
Laser Curve Extraction of Wheelset Based on Deep Learning Skeleton Extraction Network
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Efficient Frequency Estimator for a Complex Exponential Signal Based on Interpolation of Selectable DTFT Samples

1
Department of Electronic and Information, Zhongyuan University of Technology, Zhengzhou 450007, China
2
Department of Information Engineering, Zhengzhou University, Zhengzhou 450001, China
*
Author to whom correspondence should be addressed.
Sensors 2022, 22(3), 861; https://doi.org/10.3390/s22030861
Submission received: 1 December 2021 / Revised: 13 January 2022 / Accepted: 20 January 2022 / Published: 23 January 2022
(This article belongs to the Section Sensing and Imaging)

Abstract

:
The frequency estimation of complex exponential carrier signals in noise is a critical problem in signal processing. To solve this problem, a new iterative frequency estimator is presented in this paper. By iteratively computing the interpolation of DTFT samples, the proposed algorithm obtains a fine frequency estimate. In addition, its mean square error (MSE) analysis is presented in this paper. By analyzing influences of the selectable parameters on the estimation accuracy of the model, a method for choosing appropriate parameters is discussed, helping to reduce the estimation error of the proposed estimator. Simulation results show that compared with other algorithms with a comparable estimation accuracy, the proposed iterative estimator can obtain a root mean square error (RMSE) that is closer to Cramér–Rao lower bound (CRLB).

1. Introduction

Obtaining accurate and fast frequency estimates of complex exponential signals is a classic problem in signal processing. It has extensive applications in navigation [1], power systems [2,3], metabolomics [4], and aerospace communication [5,6]. The frequency offsets in these systems may lead to demodulation failure, especially for highly dynamic applications, such as satellite communication, aviation communication control systems, and missile control systems. Thus, an efficient and simple frequency estimation method is necessary to ensure real-time communication for the above applications.
Frequency estimation can be seen as a search process. A coarse frequency estimate can be obtained by searching the spectral maximum. However, the estimation accuracy of algorithms based on the spectrum calculated by a DFT is limited by the frequency resolution. The obtained estimate must be an integral multiple of the frequency resolution, which deviates from the fact that the true frequency usually lies between spectral lines. Thus, it is necessary to determine a fine frequency estimation to further reduce the residual frequency offset.
Many estimators have been proposed based on the interpolation of DFT and DTFT samples, which is an efficient way to perform fine frequency estimation [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. The interpolations of Jacobsen estimator [7], Candan estimator [8,9], Quinn estimator [10], parabolic interpolation estimator [11,12], and Fang estimator [13] are based on DFT samples. The residual frequency offset can be reduced with only a few additional multiplication operations after obtaining a coarse estimation. However, the RCSTL estimator [14], which uses three DTFT samples to assist with fine frequency estimation, obtains a smaller estimation error than other approaches. Its estimation error still deviates from the CRLB. Therefore, A&M in [16] first introduced the iteration operation into interpolation algorithms and further reduced the offset. IpDFT/e-IpDFT [19] were proposed to reduce the spectral leakage of coarse estimation by iteratively interpolating the weighted signal. The weighted signal is generated by applying different windows on the received signal [20,21,22,23]. Moreover, iterative algorithms can also obtain precise frequency estimates for various signals, such as the real sinusoid signal [26], the damped complex exponential signal [28], and the two-dimensional (2-D) complex exponential signal [29]. The asymptotic variances of some iterative estimators are very close to the CRLB, but an edge effect problem is encountered in which the variance may increase obviously at the edge of the frequency offset range. At the end, we summarize the performance of the key estimators in Table 1:
In fact, even a slight improvement in the estimation accuracy is important for signals in highly dynamic applications. Thus, to further explore the performance potential of estimators based on interpolation and obtain robust estimates, we present a new iterative algorithm in this paper. It includes two stages. In the coarse estimation stage, we first pad M N zeros after the N-point signal samples. Then, a coarse estimate can be obtained by maximizing the magnitude of the M-point FFT, which is an efficient computation method of the DFT. In the fine estimation stage, every iteration, three DTFT samples are calculated to obtain a fine estimate. The general form of the estimator and a mean square error (MSE) analysis are also presented so that the estimator parameters can be adjusted according to the given design requirements. Simulation results reveal that the RMSE of the proposed estimator is closest to the CRLB in most cases. At the same time, this estimator has a relatively simple form and low computational complexity, so it can satisfy the high accuracy and fast estimation speed requirements of the given application environments. It can not only be used in a single complex signal system but can also be used for signal estimation in multi-sinusoid communication systems when introducing the approach proposed in [27].
The rest of the paper is organized as follows. In Section 2, the deduction of the model expression is presented, and the estimation process is also described. In Section 3, a general MSE analysis of the new estimator is performed. The developed parameter design method is discussed in Section 4. In Section 5, a comparison between the new estimator and existing estimators is given by simulation.

2. Signal and Frequency Offset Model

For the frequency estimation of complex exponential signals in noise, the received signal can be written as [11]
y ( n ) = x ( n ) + w ( n ) = A e ( j ( 2 π f f s n + ϕ ) ) + w ( n ) , n = 0 , 1 , , N 1 ,
where x ( n ) = A e ( j ( 2 π f f s n + ϕ ) ) is the useful signal, A is the signal magnitude, ϕ is the initial phase, f s is the sampling rate, N is the number of signal samples, and f is the signal frequency that needs to be estimated. w ( n ) represents complex Gaussian white noise following a normal distribution N ( 0 , σ 2 ) .

2.1. CRLB of Frequency Estimation

The Cramér–Rao Lower Bound of frequency estimation is [30]
C R L B = 3 2 π 2 S N R · N ( N 2 1 ) ,
where S N R is the signal-to-noise ratio (SNR). Its definition is S N R = A 2 σ 2 .

2.2. DFT and DTFT of Complex Exponential Signals

First, the M-point DFT of x ( n ) can be obtained by
X ( k ) = n = 0 M 1 x ( n ) e j 2 π n k M .
It is assumed that M-point data are generated by padding an integral multiple of M N zeros after N signal samples; for example, M = 3 N means that N signal samples are padded with 2 N zeros. (3) can be written as
X ( k ) = n = 0 N 1 A e j ϕ e j 2 π f f s n e j 2 π n k M .
In (3), noise is not considered for simplifying the estimation process. The spectral line can be found by locating the maximum of X ( k ) and is denoted as k m , k m [ 1 , M ] . The corresponding frequency estimation is k m Δ f . Δ f = f s / M is the frequency resolution of the spectrum of the input complex exponential signal. Since the actual spectral peak usually lies between spectral lines f = ( k m ± δ ) Δ f , where δ [ 0.5 , 0.5 ] is the residual frequency offset between the true frequency and the coarse frequency estimate k m Δ f , a smaller frequency resolution is helpful for obtaining a higher frequency estimation accuracy. Thus, frequency resolution has an important impact on the accuracy of frequency estimation. To reduce this influence, it is necessary to perform a second stage to refine the frequency estimation process. This is an efficient method for performing fine estimation using DTFT samples around the DFT peak k m . DTFT samples are usually regarded as interpolations between DFT samples. The DTFT samples around k m can be expressed as
X ( k m + p ) = n = 0 N 1 A e j ϕ e j 2 π f f s n e j 2 π n k M k = k m + p , f = ( k m + δ ) Δ f .
As shown in Figure 1, the two neighboring green lines of k m are the DTFT samples used to assist with the estimation of δ . 1 < p < 1 is the offset between the assistant DTFT samples and the coarse frequency estimate k m Δ f . The meanings of Δ f , δ , k m , and p, and their relationships are also shown in Figure 1.
After some rearrangement, (5) can be simplified as
X ( k m + p ) = A e j ϕ e j π N 1 M ( δ p ) sin ( π ( δ p ) ) sin ( π M ( δ p ) ) .

3. Proposed Algorithm

To ensure the estimation accuracy of our model, we focus on the fine frequency estimation process and propose a new estimator. The steps of the proposed estimator are shown in Algorithm 1. First, we give the expressions of the two neighboring DTFT samples of the spectral peak:
X ( k m ± p ) = A e j ϕ e j 2 π N 1 M ( δ p ) sin ( π N M ( δ p ) ) sin ( π M ( δ p ) ) .
Then, we obtain their magnitudes:
X ( k m ± p ) = A sin ( π N M ( δ p ) ) sin ( π M ( δ p ) ) .
From the definitions of δ and p, it can be deduced that δ p   < 3 / 2 . Considering that 0 < N M < 1 / 2 , we can obtain 3 π / 4 < π N M ( δ p ) < 3 π / 4 . Thus, sin ( π N M ( δ p ) ) and sin ( π M ( δ p ) ) have the same sign in this range, namely, sin ( π N M ( δ p ) ) sin ( π M ( δ p ) ) > 0 . In this way, (8) can be changed to
X ( k m ± p ) = A sin ( π N M ( δ p ) ) sin ( π M ( δ p ) ) .
The ratio of X ( k m + p ) to X ( k m ) is
X ( k m + p ) X ( k m ) = sin ( π N M ( δ p ) ) sin ( π M ( δ p ) ) / sin ( π N M ( δ ) ) sin ( π M ( δ ) ) .
Since M π ( δ p ) , we obtain
X ( k m + p ) X ( k m ) sin ( π N M ( δ p ) ) sin ( π N M ( δ ) ) / δ p δ .
(11) can also be expressed as:
X ( k m + p ) X ( k m ) δ p δ sin ( π N M ( δ p ) ) sin ( π N M ( δ ) ) .
According to the transformation formula of a trigonometric function, (12) can be decomposed as
X ( k m + p ) X ( k m ) δ p δ sin ( π N M δ ) cos ( π N M p ) cos ( π N M δ ) sin ( π N M p ) sin ( π N M δ ) .
Similarly,
X ( k m p ) X ( k m ) δ + p δ sin ( π N M δ ) cos ( π N M p ) + cos ( π N M δ ) sin ( π N M p ) sin ( π N M δ ) .
After adding (13) to (14), we obtain
X ( k m + p ) X ( k m ) δ p δ + X ( k m p ) X ( k m ) δ + p δ 2 cos ( π N M p ) .
Finally, the estimate of δ can be obtained by
δ ^ = p X ( k m + p ) p X ( k m p ) X ( k m + p ) + X ( k m p ) 2 X ( k m ) cos ( π N M p ) .
The frequency estimate is
f ^ = ( k m + δ ^ ) f s M .
From (17), we can obtain a precise estimate of the frequency offset, but the above analysis is based on a signal without noise and some approximations. In fact, the frequency estimation derived from (17) usually deviates from the true frequency due to the influence of noise, so we introduce an iteration operation into the fine estimation process. The procedures of the proposed estimator are shown in Algorithm 1.
Algorithm 1 The proposed iterative estimator.
Input: Signal samples y
Output: Frequency estimate f ^
1:
Coarse estimation: calculate the M-point DFT, and obtain the location of the spectral peak k m as the coarse frequency estimate.
2:
Fine estimation: initialize the iteration variables; δ ^ = 0 , q = 1
3:
for  q < Q do
4:
      Compute the DTFT samples X ( k m + p ) , X ( k m ) , and  X ( k m p ) according to (8)
5:
      Compute δ ^ according to (16)
6:
      Update the iteration variables: k m = k m + δ ^ , q = q + 1
7:
Compute f ^ according to (17)
8:
return  f ^

4. MSE Analysis of the Proposed Algorithm

To perform a MSE analysis for the developed iterative algorithm, we give the expression of the spectral magnitude of a complex exponential signal in noise:
Y ( k ) = X ( k ) + W ( k ) ,
where W ( k ) is the DFT of additive Gaussian white noise w ( n ) . The spectral peak in the noise is
Y ( k m ) = X ( k m ) + W ( k m ) = A sin π N M δ sin π M δ + W ( k m ) e j ϕ e j π N 1 M δ .
Because W ( k ) X ( k ) at a high SNR, | Y k m | = A k m + R e ( e j φ k m W k m ), A k m is the magnitude of X k m , and φ k m is the phase of X k m . Therefore, an approximation of (19) can be expressed as
Y ( k m ) X ( k m ) + Re ( W ( k m ) e j φ k m ) A k m + U k m ,
where A k m =   X ( k m ) = A · sin π N M δ / sin π M δ , U k m = R e ( W ( k m ) e j φ k m ) , φ k m = ϕ + π N 1 M δ . Similarly, the DTFT samples can be expressed as Y ( k m + p ) A k m + p + U k m + p , Y ( k m p ) A k m p + U k m p . To simplify the analysis, we denote Y 0 = Y ( k m ) , Y p = Y ( k m + p ) , Y p = Y ( k m p ) , A 0 = A k m , A p = A k m + p , A p = A k m p , U 0 = U k m , U p = U k m + p , and U p = U k m p . Consequently, (16) becomes
δ ^ = p ( A p + U p ) p ( A p + U p ) ( A p + U p ) + ( A p + U p ) 2 ( A 0 + U 0 ) cos π N M p = p ( A p A p + U p U p ) A p + A p 2 A 0 cos π N M p + U p + U p 2 U 0 cos π N M p = p ( A p A p + U p U p ) A p + A p 2 A 0 cos π N M p · 1 1 + U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p .
Because U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p 1 , the first-order Taylor expansion of the second item on the right side of (21) is
1 1 + U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p 1 U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p .
By introducing (22) into (21), the frequency offset is
δ ^ = p ( A p A p ) A p + A p 2 A 0 cos π N M p + p ( U p U p ) A p + A p 2 A 0 cos π N M p p ( A p A p ) U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p 2 p ( U p U p ) U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p 2 .
The last item is approximately zero if the SNR is high. Thus, after some simplification, (23) becomes
δ ^ = δ + p ( U p U p ) A p + A p 2 A 0 cos π N M p δ · U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p .
Then, we compute the expectation of the estimation error and the MSE of the frequency estimate:
E ( δ ^ δ ) = p ( U p U p ) A p + A p 2 A 0 cos π N M p δ · U p + U p 2 U 0 cos π N M p A p + A p 2 A 0 cos π N M p = ( p δ ) U p ( p + δ ) U p ) 2 U 0 δ cos π N M p A p + A p 2 A 0 cos π N M p ,
and
E ( δ ^ δ ) 2 = E ( p δ ) U p ( p + δ ) U p ) A p + A p 2 A 0 cos π N M p 2 U 0 δ cos π N M p A p + A p 2 A 0 cos π N M p 2 .
According to deduction in Appendix A, the final MSE expression is obtained as follows:
E ( δ ^ δ ) 2 π 2 N M 2 · δ 2 ( δ 2 p 2 ) 2 S N R · p 2 + δ 2 + 2 δ 2 cos 2 π N M p ( 3 δ 2 + p 2 ) sin c 2 π N M p 4 p 2 p sin π N M δ cos π N M p δ cos π N M δ sin π N M p 2 , δ p & δ 0 p 2 4 N · S N R · 1 sin c 2 π N M p sin c π N M p cos π N M p 2 , δ = 0 2 p 2 N · S N R · 1 + cos 2 π N M p 2 sin c 2 π N M p 1 sin c 2 π N M p 2 , δ = p .

4.1. Discussion on p

Once the ratio M / N has been decided, a smaller p means that the assistant DTFT samples are closer to the DFT peak. Will we obtain a more precise estimate if a smaller p is chosen? How can we choose p? This section provides a discussion about p. We simulate the proposed method versus p for different values of δ to verify the influence of p on the estimation accuracy of the resulting model.
In Figure 2, p has an important influence on the RMSE/CRLB of the algorithm. It decreases as p decreases at the beginning. However, when p is sufficiently small, the improvement in the RMSE/CRLB does not change further. In fact, it does not keep decreasing as p decreases. The minimum RMSE/CRLB varies with δ , so we choose a medium p = 0.3 in Section 5.

4.2. Discussion on M

The spectral frequency resolution directly affects the algorithm’s estimation accuracy. When the sampling frequency is determined, choosing a larger N can reduce the spectral frequency resolution and increase the spectral magnitude. However, this also increases the estimation time and computational complexity of the algorithm. Therefore, we first pad zeros after the signal samples before computing the corresponding DFTs. This approach can decrease the frequency resolution without increasing the observation time. The DFT and DTFT waveforms of samples with and without zero-padding are shown in Figure 3. The figure indicates that the DFT and DTFT waveforms of the signal are extended after the signal samples are padded with N zeros. The two green circles are closer to the peak than the two blue circles. The 2 N -point DFT around the peak, which is closer to the peak, can obtain a smaller RMSE/CRLB in Figure 4. In addition, as M / N increases, the improvement in the MSE decreases. Considering the computational complexity brought by M, we choose M = 2 N as a compromise between computational complexity and the estimation accuracy.

4.3. Discussion on Q

We introduce the iteration operation into the frequency estimation and decrease the final estimation error by iteratively updating the estimate of the frequency offset and the DTFT samples used in interpolation. To choose the number of iteration, we make a simulation about the performance of the proposed algorithm versus iteration number Q at different δ when S N R = 0 dB and N = 512 . The RMSE/CRLB is shown in Figure 5. From Figure 5, it can be seen that the RMSE/CRLB decreases obviously after the second iteration, when δ = 0.5 or δ = 0.3 . However, When Q > 2 , the RMSE/CRLB does not decrease as Q increases. Thus, considering the accuracy improvement and the calculation complexity, we choose Q = 2 as the iteration number for simulation in Section 5. This iteration number is also consistent with the research in [24].

5. Experiment

To test the estimation performance of the proposed iterative estimator, we set the simulation parameters as follows: the signal magnitude A = 1 , the sampling frequency f s = 512 kHz, the frequency offset δ < 0.5 , the spectral peak k m = 64 kHz, the true frequency 63.5 kHz < f 0 < 64.5 kHz, the number of iterations Q = 2 [24], and the initial phase ϕ [ 0 , 2 π ] is generated randomly. Moreover, we compare the proposed algorithm with other algorithms including Jacobsen estimator [7], Candan estimator [8,9], A&M estimator [16], Fan estimator [15], and Fang estimator [13], and RCSTL [14]. Every estimation result is obtained from 10,000 simulations.

5.1. Comparison of the CRLB and RMSE of the Proposed Estimator

In this section, we perform simulations of the proposed estimator with different iterations. As shown in Figure 6, the simulation result at the first iteration is consistent with the MSE analysis. The periodic characteristic is also verified. Therefore, the MSE analysis in Section 4 is effective. However, these results provide only the reference standard for the proposed estimator at the first iteration. The RMSE/CRLB at the second iteration can be verified by simulation. In Figure 6, the RMSE/CRLB at the second iteration can approach 1, which means that the proposed estimator has a promising estimation accuracy.
One more aspect that needs to be explained is that the MSE of the proposed estimator in (27) is a periodic function. It can be seen in Figure 6 that between δ = 0.5 and δ = 0.5 , there are two periods for M = 2 N and Q = 1 and one period for M = N and Q = 1 . They are both symmetric about δ = 0 . This means that the RMSE/CLRB is periodical and that the length of its period varies with the ratio of M to N.

5.2. Comparison of RMSEs among the Existing Algorithms

To test the proposed estimator’s RMSE versus the SNR and N, we set δ = 0.2 and simulate seven estimators.
Figure 7, Figure 8, Figure 9 and Figure 10 are figures corresponding to N = 64 , N = 128 , N = 256 , and N = 512 , respectively. Upon comparing these figures, the SNR threshold decreases from 0 dB to 10 dB. The RMSEs of these estimators and the RMSE differences among these estimators decrease as N increases. Detailed figures at S N R = 10 dB are also provided, from which we can see that the RMSEs of A&M, Fan, and the proposed estimator are closest to the CRLB. The RMSE of the proposed estimator is slightly lower than those of the other simulated methods in most cases, although the RMSE difference between the proposed estimator and the Fan estimator is small. In addition, the RMSE of the proposed estimator is only 1.003 times the CRLB when S N R = 10 dB and N = 512 .

5.3. RMSE Comparison among Estimators with a Comparable Accuracy

To further evaluate the RMSE of the proposed estimator, we compare the estimators with the comparable estimation accuracy at S N R = 0 dB. The results in Figure 11 and Figure 12 show that the RMSE of the proposed estimator is closest to the CRLB in most cases. RMSE of the Fan estimator is slightly larger than that of the proposed estimator. Both the two estimators are more stable than the other simulated estimators, showing that the variations in their RMSEs are small throughout the whole δ range. The RMSE of the A&M estimator is also very small, but there exists a kind of edge effect in which the RMSE increases at δ = ± 0.5 . For some values of δ , the corresponding RMSEs may lie slightly below the CRLB and are asymmetric relative to δ = 0 due to the limited number of Monte Carlo trials.

5.4. Computational Complexities of the Estimators Used in the Simulation

Computational complexity is an important aspect to consider when evaluating the performance of estimators. We mainly compare the numbers of complex calculations performed by the various methods, including complex multiplications and complex additions.
Table 2 gives the computational complexities of all simulated estimators. To reduce the computational complexity of the DFT, we choose N = 2 k as the DFT length and compute the FFT instead. An N-point FFT requires N 2 l o g 2 N complex multiplications and N l o g 2 N complex additions. A one-bin N-point DTFT requires N complex multiplications and N 1 complex additions. The additional multiplication and addition operation other than those utilized in the FFT calculation are performed to compute δ ^ . During the two iterations, the proposed estimator uses five magnitudes of DTFTs to compute δ ^ , so it omits several complex multiplications and additions compared with the Fan estimator, which have comparable estimation accuracies as the proposed estimator. Compared with the A&M estimator, the proposed estimator has a lower RMSE and avoids the edge effect when δ approaches to 0.5.

6. Conclusions

In this paper, we present a new and simple iterative frequency estimator that is used for the fine estimation of complex exponential signals. In this algorithm, the magnitudes of DTFT samples are used to calculate the frequency offset to reduce the number of required complex calculations, and iterative calculations are introduced into the estimation process to further improve the estimation accuracy of the resulting model. It is verified by simulation that the developed algorithm achieves a higher estimation accuracy compared with other algorithms. Moreover, it has a lower computational complexity than the algorithm that has a comparable estimation accuracy. In addition, the general form of the proposed estimator is also presented in this paper to help select appropriate parameters according to given application requirements.

Author Contributions

Conceptualization, M.W.; methodology, M.W.; software, B.L. and J.S.; validation, B.L. and J.S.; formal analysis, M.W.; investigation, J.S.; resources, B.L.; data curation, J.S.; writing—original draft preparation, M.W.; writing—review and editing, A.Z.; visualization, B.L.; supervision, L.Q.; project administration, L.Q.; funding acquisition, A.Z. All authors have read and agreed to the published version of the manuscript.

Funding

The National Natural Science Foundation of China (No. 61901537).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data available on request due to restrictions. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Deduction of the MSE Expression

The deduction process of the final expression of the MSE is given in this section. To obtain E ( δ ^ δ ) 2 , we need to compute the numerator and denominator of the MSE:
( δ ^ δ ) 2 n o m i = ( p δ ) 2 U p 2 + ( p + δ ) 2 U p 2 + 4 U 0 2 δ 2 cos 2 π N M p 2 ( p 2 δ 2 ) U p U p + 4 U p U + p ( p δ ) δ cos π N M p + 4 U p U p ( p + δ ) δ cos π N M p ,
and
( δ ^ δ ) 2 de n o m i = A p + A p 2 A 0 cos π N M p 2 .
Because U p is the real part of W p at all frequencies and w ( n ) N ( 0 , σ 2 ) , E ( U p 2 ) = E ( U p 2 ) = E ( U 0 2 ) = N σ 2 / 2 . To compute E ( U p U + p ) , we go back to the DTFT expression:
U p = Re ( W p e j φ p ) = Re ( n = 0 N 1 w ( n ) e j ( 2 π f p n + φ p ) ) .
The expectation of U p U p is
E U p U p = E n = 0 N 1 w r ( n ) cos ( 2 π f p n + φ p ) + w i ( n ) sin ( 2 π f p n + φ p ) · m = 0 N 1 w r ( m ) cos ( 2 π f p m + φ p ) + w i ( m ) sin ( 2 π f p m + φ p ) ,
where φ p = φ + π N 1 M ( δ p ) and φ p = φ + π N 1 M ( δ + p ) are the phases of X p and X p , respectively. f p = N M ( δ p ) and f p = N M ( δ + p ) are the frequencies of X p and X p , respectively. According to the independence characteristic of Gaussian white noise, E ( w r ( n ) w i ( m ) ) = 0 , and E ( w i ( n ) w r ( m ) ) = 0 . If n m , E ( w r ( n ) w r ( m ) ) = 0 , and E ( w i ( n ) w i ( m ) ) = 0 . Therefore, only autocorrelation items ( n = m ) are retained in the expectation, which can be written as
E U p U p = E [ n = 0 N 1 w r 2 ( n ) cos ( 2 π f p n + φ p ) cos ( 2 π f p n + φ p ) + w i 2 ( n ) sin ( 2 π f p n + φ p ) sin ( 2 π f p n + φ p ) ] = N σ 2 / 2 n = 0 N 1 cos ( 2 π f p n + φ p ) cos ( 2 π f p n + φ p ) + sin ( 2 π f p n + φ p ) sin ( 2 π f p n + φ p ) = N σ 2 / 2 n = 0 N 1 cos ( 2 π ( f p f p ) n + ( φ p φ p ) ) .
After some algebra, we have
E U p U p = N σ 2 4 e j ( φ p φ p ) 1 e j ( 2 π ( f p f p ) n ) 1 e j ( 2 π ( f p f p ) n / M ) + e j ( φ p φ p ) 1 e j ( 2 π ( f p f p ) n ) 1 e j ( 2 π ( f p f p ) n / M ) .
By substituting φ p φ p = π N 1 M p , f p f p = 2 N M p into (A6), we obtain
E U p U p = N σ 2 2 sin ( 2 π N p M ) sin ( 2 π p M ) N σ 2 2 sin ( 2 π N p M ) 2 π p M .
In the same way, E U 0 U p = E U 0 U p = N σ 2 / 2 sin ( π N p / M ) sin ( π p / M ) N σ 2 / 2 sin ( π N p / M ) π p / M . Then, we introduce the above results into (A1) and obtain
E ( δ ^ δ ) 2 n o m i = N σ 2 p 2 + δ 2 + 2 δ 2 cos 2 π N M p ( 3 δ 2 + p 2 ) sin 2 π N M p 2 π M p .
To compute the denominator of E ( δ ^ δ ) 2 , we return to the expression of A p
A p = X ( k m + p ) = A · sin π N M ( δ p ) sin π M ( δ p ) A · sin π N M ( δ p ) π M ( δ p ) .
Similarly, A p A · sin π N M ( δ + p ) π M ( δ + p ) , A 0 A · sin π N M δ π M δ .
By introducing A p , A p , and A 0 into (A2), we obtain
A p + A p 2 A 0 cos π N M p A · sin π N M ( δ p ) π M ( δ p ) + sin π N M ( δ + p ) π M ( δ + p ) 2 sin π N M δ cos π N M p π M δ 2 A M π · p 2 sin π N M δ cos π N M p δ ( δ 2 p 2 ) p δ cos π N M δ sin π N M p δ ( δ 2 p 2 ) .
Then,
E ( δ ^ δ ) 2 de n o m i = A p + A p 2 A 0 cos π N M p 2 = 2 A M p π 2 · p sin π N M δ cos π N M p δ 2 ( δ 2 p 2 ) 2 δ cos π N M δ sin π N M p 2 δ 2 ( δ 2 p 2 ) 2 .
Since there are two cases in which the denominator equals zero, which leads to a meaningless result, we further analyze these two cases. When δ = p ,
A p + A p 2 A 0 cos π N M p 2 A M p π · p sin π N M δ cos π N M p δ cos π N M δ sin π N M p δ ( δ 2 p 2 ) 2 A M p π · δ sin π N M ( δ p ) ( δ p ) sin π N M δ cos π N M p δ ( δ 2 p 2 ) 2 A N p δ + p · sin c π N M ( δ p ) sin c π N M δ cos π N M p A N · 1 sin c π N M δ cos π N M p ;
when δ = 0 ,
A p + A p 2 A 0 cos π N M p 2 A M π · p 2 sin π N M δ cos π N M p p δ cos π N M δ sin π N M p δ ( δ 2 p 2 ) 2 A M π · sin π N M δ cos π N M p δ + cos π N M δ sin π N M p p 2 A N · sin c π N M p cos π N M p .

References

  1. Cheung, K.; Lee, C.; Jun, W.; Lightsey, G. Single-Satellite Doppler Localization with Law of Cosines (LOC). In Proceedings of the 2019 IEEE Aerospace Conference, Big Sky, MT, USA, 2–9 March 2019; pp. 1–12. [Google Scholar] [CrossRef]
  2. Zhang, J.; Tang, L.; Mingotti, A.; Peretto, L.; Wen, H. Analysis of White Noise on Power Frequency Estimation by DFT-Based Frequency Shifting and Filtering Algorithm. IEEE Trans. Instrum. Meas. 2020, 69, 4125–4133. [Google Scholar] [CrossRef]
  3. Sun, J.; Aboutanios, E.; Smith, D.B.; Fletcher, J.E. Robust Frequency, Phase, and Amplitude Estimation in Power Systems Considering Harmonics. IEEE Trans. Power Deliv. 2020, 35, 1158–1168. [Google Scholar] [CrossRef]
  4. Aboutanios, E.; Kopsinis, Y.; Rubtsov, D. Instantaneous frequency based spectral analysis of nuclear magnetic spectroscopy data for metabolomics. In Proceedings of the 2010 3rd International Congress on Image and Signal Processing, Yantai, China, 16–18 October 2010; Volume 7, pp. 3428–3432. [Google Scholar] [CrossRef]
  5. Divsalar, D.; Net, M.S.; Cheung, K.M. Adaptive Sweeping Carrier Acquisition and Tracking for Dynamic Links with High Uplink Doppler. In Proceedings of the 2020 IEEE Aerospace Conference, Big Sky, MT, USA, 7–14 March 2020; pp. 1–14. [Google Scholar] [CrossRef]
  6. Wang, J.; Jiang, C.; Kuang, L.; Guo, S. Iterative Doppler Frequency Offset Estimation in Low SNR Satellite Communications. In Proceedings of the 2020 International Wireless Communications and Mobile Computing (IWCMC), Limassol, Cyprus, 15–19 June 2020; pp. 36–41. [Google Scholar] [CrossRef]
  7. Jacobsen, E.; Kootsookos, P. Fast, Accurate Frequency Estimators [DSP Tips & Tricks]. IEEE Signal Process. Mag. 2007, 24, 123–125. [Google Scholar]
  8. Candan, C. A Method For Fine Resolution Frequency Estimation From Three DFT Samples. IEEE Signal Process. Lett. 2011, 18, 351–354. [Google Scholar] [CrossRef]
  9. Candan, C. Analysis and Further Improvement of Fine Resolution Frequency Estimation Method From Three DFT Samples. IEEE Signal Process. Lett. 2013, 20, 913–916. [Google Scholar] [CrossRef]
  10. Quinn, B.G. Estimation of frequency, amplitude, and phase from the DFT of a time series. IEEE Trans. Signal Process. 1997, 45, 814–817. [Google Scholar] [CrossRef]
  11. Djukanović, S. Sinusoid frequency estimator with parabolic interpolation of periodogram peak. In Proceedings of the 2017 40th International Conference on Telecommunications and Signal Processing (TSP), Barcelona, Spain, 5–7 July 2017; pp. 470–473. [Google Scholar] [CrossRef]
  12. Djukanović, S.; Popović, T.; Mitrović, A. Precise sinusoid frequency estimation based on parabolic interpolation. In Proceedings of the 2016 24th Telecommunications Forum (TELFOR), Belgrade, Serbia, 22–23 November 2016; pp. 1–4. [Google Scholar] [CrossRef]
  13. Fang, L.; Duan, D.; Yang, L. A new DFT-based frequency estimator for single-tone complex sinusoidal signals. In Proceedings of the MILCOM 2012—2012 IEEE Military Communications Conference, Orlando, FL, USA, 29 October–1 November 2012; pp. 1–6. [Google Scholar] [CrossRef]
  14. Yang, C.; Wei, G. A Noniterative Frequency Estimator with Rational Combination of Three Spectrum Lines. IEEE Trans. Signal Process. 2011, 59, 5065–5070. [Google Scholar] [CrossRef]
  15. Fan, L.; Qi, G.; Xing, J.; Jin, J.; Liu, J.; Wang, Z. Accurate Frequency Estimator of Sinusoid Based on Interpolation of FFT and DTFT. IEEE Access 2020, 8, 44373–44380. [Google Scholar] [CrossRef]
  16. Aboutanios, E.; Mulgrew, B. Iterative frequency estimation by interpolation on Fourier coefficients. IEEE Trans. Signal Process. 2005, 53, 1237–1242. [Google Scholar] [CrossRef]
  17. Chen, X.; Zhi, W.; Xiao, Y. Carrier Frequency Estimation based on Frequency Domain and Wavelet Ridge. Int. J. Perform. Eng. 2019, 15, 2107. [Google Scholar]
  18. Liu, W.; Bian, X.; Deng, Z.; Mo, J.; Jia, B. A Novel Carrier Loop Algorithm Based on Maximum Likelihood Estimation (MLE) and Kalman Filter (KF) for Weak TC-OFDM Signals. Sensors 2018, 18, 2256. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Belega, D.; Petri, D. Effect of Noise and Harmonics on Sine-wave Frequency Estimation by Interpolated DFT Algorithms Based on Few Observed Cycles. Signal Process. 2017, 140, 207–218. [Google Scholar] [CrossRef]
  20. Belega, D.; Petri, D. Influence of the Noise on DFT-based Sine-wave Frequency and Amplitude Estimators. Measurement 2019, 137, 527–534. [Google Scholar] [CrossRef]
  21. Belega, D.; Petri, D.; Dallet, D. Accurate Frequency Estimation of a Noisy Sine-wave by Means of an Interpolated Discrete-Time Fourier Transform Algorithm. Measurement 2017, 116, 685–691. [Google Scholar] [CrossRef]
  22. Belega, D.; Petri, D. Accuracy Analysis of an Enhanced Frequency-Domain Linear Least-Squares Algorithm. In Proceedings of the 2019 IEEE International Instrumentation and Measurement Technology Conference (I2MTC), Auckland, New Zealand, 20–23 May 2019; pp. 1–5. [Google Scholar] [CrossRef]
  23. Belega, D.; Petri, D. Accuracy of DTFT-based Sine-wave Amplitude Estimators. J. Phys. Conf. Ser. 2018, 1065, 052030. [Google Scholar] [CrossRef]
  24. Serbes, A. Fast and Efficient Sinusoidal Frequency Estimation by Using the DFT Coefficients. IEEE Trans. Commun. 2019, 67, 2333–2342. [Google Scholar] [CrossRef]
  25. Belega, D.; Petri, D. Effect of Frequency Uncertainty on the Sine-wave Amplitude Estimator Returned by Linear Least-Squares Fitting. In Proceedings of the 2020 IEEE International Instrumentation and Measurement Technology Conference (I2MTC), Dubrovnik, Croatia, 25–28 May 2020; pp. 1–5. [Google Scholar] [CrossRef]
  26. Djukanović, S. An Accurate Method for Frequency Estimation of a Real Sinusoid. IEEE Signal Process. Lett. 2016, 23, 915–918. [Google Scholar] [CrossRef]
  27. Djukanović, S.; Popović-Bugarin, V. Efficient and Accurate Detection and Frequency Estimation of Multiple Sinusoids. IEEE Access 2019, 7, 1118–1125. [Google Scholar] [CrossRef]
  28. Aboutanios, E.; Ye, S. Efficient Iterative Estimation of the Parameters of a Damped Complex Exponential in Noise. IEEE Signal Process. Lett. 2014, 21, 975–979. [Google Scholar] [CrossRef]
  29. Ye, S.; Aboutanios, E. Two dimensional frequency estimation by interpolation on Fourier coefficients. In Proceedings of the 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, Japan, 25–30 March 2012; pp. 3353–3356. [Google Scholar] [CrossRef]
  30. Rife, D.; Boorstyn, R. Single tone parameter estimation from discrete-time observations. IEEE Trans. Inf. Theory 1974, 20, 591–598. [Google Scholar] [CrossRef]
Figure 1. Magnitudes of the DTFT and DFT of a complex exponential signal.
Figure 1. Magnitudes of the DTFT and DFT of a complex exponential signal.
Sensors 22 00861 g001
Figure 2. The RMSE/CRLB of the proposed algorithm versus p for different values of δ when the S N R = 0 dB and N = 512 .
Figure 2. The RMSE/CRLB of the proposed algorithm versus p for different values of δ when the S N R = 0 dB and N = 512 .
Sensors 22 00861 g002
Figure 3. DFTs and DTFTs of complex exponential signals at M = N and M = 2 N .
Figure 3. DFTs and DTFTs of complex exponential signals at M = N and M = 2 N .
Sensors 22 00861 g003
Figure 4. The RMSE/CRLB of the proposed estimator versus δ at different values of M when S N R = 0 dB and N = 512 .
Figure 4. The RMSE/CRLB of the proposed estimator versus δ at different values of M when S N R = 0 dB and N = 512 .
Sensors 22 00861 g004
Figure 5. The RMSE/CRLB of the proposed estimator versus Q at different values of δ when S N R = 0 dB and N = 512 .
Figure 5. The RMSE/CRLB of the proposed estimator versus Q at different values of δ when S N R = 0 dB and N = 512 .
Sensors 22 00861 g005
Figure 6. The RMSE/CRLB of the proposed estimator versus δ under different settings of M and Q when S N R = 0 dB and N = 512 .
Figure 6. The RMSE/CRLB of the proposed estimator versus δ under different settings of M and Q when S N R = 0 dB and N = 512 .
Sensors 22 00861 g006
Figure 7. RMSE comparison among frequency estimates versus SNR at N = 64 and δ = 0.2 .
Figure 7. RMSE comparison among frequency estimates versus SNR at N = 64 and δ = 0.2 .
Sensors 22 00861 g007
Figure 8. RMSE comparison among frequency estimates versus SNR at N = 128 and δ = 0.2 .
Figure 8. RMSE comparison among frequency estimates versus SNR at N = 128 and δ = 0.2 .
Sensors 22 00861 g008
Figure 9. RMSE comparison among frequency estimates versus SNR at N = 256 and δ = 0.2 .
Figure 9. RMSE comparison among frequency estimates versus SNR at N = 256 and δ = 0.2 .
Sensors 22 00861 g009
Figure 10. RMSE comparison among frequency estimates versus SNR at N = 512 and δ = 0.2 .
Figure 10. RMSE comparison among frequency estimates versus SNR at N = 512 and δ = 0.2 .
Sensors 22 00861 g010
Figure 11. RMSE comparison among the frequency estimates versus δ at N = 256 and S N R = 0 dB.
Figure 11. RMSE comparison among the frequency estimates versus δ at N = 256 and S N R = 0 dB.
Sensors 22 00861 g011
Figure 12. RMSE comparison among the frequency estimates versus δ at N = 512 and S N R = 0 dB.
Figure 12. RMSE comparison among the frequency estimates versus δ at N = 512 and S N R = 0 dB.
Sensors 22 00861 g012
Table 1. Analysis of key estimators.
Table 1. Analysis of key estimators.
EstimatorsMerits/DemeritsLimits of Applicability
IpDFT with a cosine window [20]Reduce the influence of spectral leakage by introducing a cosine window/increase the computationLinear computation based on all DFT samples
IpDTFT with a cosine window [21]Reduce the influence of spectral leakage by introducing a cosine window/increase the computation, extra two DTFT computationsAlmost two times CRLB
e-FLLs/FLLs (frequency-domain linear least-squares) [22,23]Robust to harmonics/matrix computation for window, matrix inversion for linear least-squaresAmplitude and phase estimation
QSE/HAQSE (hybrid A&M and q-shift estimator) [24]DFT-based, easily realizable, within ±0.003 dB CRLB/at least extra four DTFT computations and several complex computations for the estimateEdge effect
Aboutanios and Mulgrew (A&M) [25,26]DFT-based, easily realizable, 1.0147 times the asymptotic CRLB/extra four DTFT computations and several complex computations for the estimateEdge effect
Parabolic interpolation [27]DFT-based, easily realizable, within ±0.526 dB CRLB/extra three DTFT computations and several complex computations for the estimateLimited accuracy
Table 2. Computational complexities of the estimators used in the simulation.
Table 2. Computational complexities of the estimators used in the simulation.
EstimatorsComplex MultiplicationsComplex Additions
Jacobsen (N points) N 2 l o g 2 N + 1 N l o g 2 N + 3
Candan (N points) N 2 l o g 2 N + 1 N l o g 2 N + 3
A&M (N points, 2 iterations) N 2 l o g 2 N + 4 N + 2 N l o g 2 N + 4 N
Fan (M points, 2 iterations) M 2 l o g 2 M + 5 N + 8 M l o g 2 M + 5 N + 1
Fang (M points) M 2 l o g 2 M M l o g 2 M
RCSTL (M points) M 2 l o g 2 M + 1 M l o g 2 M
Proposed (M points, 2 iterations) M 2 l o g 2 M + 5 N M l o g 2 M + 5 ( N 1 )
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wei, M.; Zhang, A.; Qi, L.; Li, B.; Sun, J. An Efficient Frequency Estimator for a Complex Exponential Signal Based on Interpolation of Selectable DTFT Samples. Sensors 2022, 22, 861. https://doi.org/10.3390/s22030861

AMA Style

Wei M, Zhang A, Qi L, Li B, Sun J. An Efficient Frequency Estimator for a Complex Exponential Signal Based on Interpolation of Selectable DTFT Samples. Sensors. 2022; 22(3):861. https://doi.org/10.3390/s22030861

Chicago/Turabian Style

Wei, Miaomiao, Aihua Zhang, Lin Qi, Bicao Li, and Jun Sun. 2022. "An Efficient Frequency Estimator for a Complex Exponential Signal Based on Interpolation of Selectable DTFT Samples" Sensors 22, no. 3: 861. https://doi.org/10.3390/s22030861

APA Style

Wei, M., Zhang, A., Qi, L., Li, B., & Sun, J. (2022). An Efficient Frequency Estimator for a Complex Exponential Signal Based on Interpolation of Selectable DTFT Samples. Sensors, 22(3), 861. https://doi.org/10.3390/s22030861

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop