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

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).


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 . 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 twodimensional (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: Table 1. Analysis of key estimators.

Merits/Demerits Limits of Applicability
IpDFT with a cosine window [20] Reduce the influence of spectral leakage by introducing a cosine window/increase the computation Linear 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 computations 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.

Signal and Frequency Offset Model
For the frequency estimation of complex exponential signals in noise, the received signal can be written as [11] where x(n) = Ae (j(2π f fs 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 ).

CRLB of Frequency Estimation
The Cramér-Rao Lower Bound of frequency estimation is [30] where SNR is the signal-to-noise ratio (SNR). Its definition is SNR = A 2 σ 2 .

DFT and DTFT of Complex Exponential Signals
First, the M-point DFT of x(n) can be obtained by 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 = 3N means that N signal samples are padded with 2N zeros. (3) can be written as Ae jφ e j2π f fs n e −j 2πnk 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 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

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: Then, we obtain their magnitudes: From the definitions of δ and p, it can be deduced that |δ ∓ p| < 3/2. Considering that 0 ) and sin( π M (δ ∓ p)) have the same sign in this range, namely, In this way, (8) can be changed to The ratio of |X( Since M π(δ − p), we obtain (11) can also be expressed as: . (12) According to the transformation formula of a trigonometric function, (12) can be decomposed as Similarly, After adding (13) to (14), we obtain Finally, the estimate of δ can be obtained bŷ The frequency estimate isf 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.

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: where W(k) is the DFT of additive Gaussian white noise w(n). The spectral peak in the noise is Because |W(k)| |X(k)| at a high SNR, |Y k m | = A k m + Re(e −jϕ km 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 where Similarly, the DTFT samples can be expressed as Because 1, the first-order Taylor expansion of the second item on the right side of (21) is By introducing (22) into (21), the frequency offset iŝ The last item is approximately zero if the SNR is high. Thus, after some simplification, (23) becomeŝ Then, we compute the expectation of the estimation error and the MSE of the frequency estimate: and According to deduction in Appendix A, the final MSE expression is obtained as follows:

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.

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 2N-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 = 2N as a compromise between computational complexity and the estimation accuracy.

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 SNR = 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].

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 = 2N 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.

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. 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 SNR = 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 SNR = 10 dB and N = 512.

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 SNR = 0 dB. The results in Figures 11 and 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.

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 log 2 N complex multiplications and Nlog 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.

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.  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.

(A5)
After some algebra, we have To compute the denominator of E (δ − δ) 2 , we return to the expression of A p when δ = 0,