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

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:

Table 1

Analysis of key estimators.

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	Almost 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-squares	Amplitude 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 estimate	Edge 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 estimate	Edge effect
Parabolic interpolation [27]	DFT-based, easily realizable, within ±0.526 dB CRLB/extra three DTFT computations and several complex computations for the estimate	Limited accuracy

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 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] where is the useful signal, A is the signal magnitude, is the initial phase, is the sampling rate, N is the number of signal samples, and f is the signal frequency that needs to be estimated. represents complex Gaussian white noise following a normal distribution .

2.1. CRLB of Frequency Estimation

The Cramér–Rao Lower Bound of frequency estimation is [30] where is the signal-to-noise ratio (SNR). Its definition is .

2.2. DFT and DTFT of Complex Exponential Signals

First, the M-point DFT of can be obtained by

It is assumed that M-point data are generated by padding an integral multiple of zeros after N signal samples; for example, means that N signal samples are padded with zeros. (3) can be written as

In (3), noise is not considered for simplifying the estimation process. The spectral line can be found by locating the maximum of and is denoted as , . The corresponding frequency estimation is . is the frequency resolution of the spectrum of the input complex exponential signal. Since the actual spectral peak usually lies between spectral lines , where is the residual frequency offset between the true frequency and the coarse frequency estimate , 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 . DTFT samples are usually regarded as interpolations between DFT samples. The DTFT samples around can be expressed as

As shown in Figure 1, the two neighboring green lines of are the DTFT samples used to assist with the estimation of . is the offset between the assistant DTFT samples and the coarse frequency estimate . The meanings of , , , and p, and their relationships are also shown in Figure 1.

Figure 1

Magnitudes of the DTFT and DFT of a complex exponential signal.

After some rearrangement, (5) can be simplified as

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:

Then, we obtain their magnitudes:

From the definitions of and p, it can be deduced that . Considering that , we can obtain . Thus, and have the same sign in this range, namely, . In this way, (8) can be changed to

The ratio of to is

Since , we obtain

(11) can also be expressed as:

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 by

The frequency estimate is

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.

Input: Signal samples y

Output: Frequency estimate

Coarse estimation: calculate the M-point DFT, and obtain the location of the spectral peak as the coarse frequency estimate.

Fine estimation: initialize the iteration variables;

for  do

Compute the DTFT samples , , and  according to (8)

Compute according to (16)

Update the iteration variables:

Compute according to (17)

return

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: where is the DFT of additive Gaussian white noise . The spectral peak in the noise is

Because at a high SNR, ), is the magnitude of , and is the phase of . Therefore, an approximation of (19) can be expressed as where , , . Similarly, the DTFT samples can be expressed as , . To simplify the analysis, we denote , , , , , , , , and . Consequently, (16) becomes

Because , the first-order Taylor expansion of the second item on the right side of (21) is

By introducing (22) into (21), the frequency offset is

The last item is approximately zero if the SNR is high. Thus, after some simplification, (23) becomes

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:

4.1. Discussion on p

Once the ratio 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 in Section 5.

Figure 2

The RMSE/CRLB of the proposed algorithm versus p for different values of when the dB and .

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 -point DFT around the peak, which is closer to the peak, can obtain a smaller RMSE/CRLB in Figure 4. In addition, as increases, the improvement in the MSE decreases. Considering the computational complexity brought by M, we choose as a compromise between computational complexity and the estimation accuracy.

Figure 3

DFTs and DTFTs of complex exponential signals at and .

Figure 4

The RMSE/CRLB of the proposed estimator versus at different values of M when dB and .

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 dB and . 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 or . However, When , the RMSE/CRLB does not decrease as Q increases. Thus, considering the accuracy improvement and the calculation complexity, we choose as the iteration number for simulation in Section 5. This iteration number is also consistent with the research in [24].

Figure 5

The RMSE/CRLB of the proposed estimator versus Q at different values of when dB and .

5. Experiment

To test the estimation performance of the proposed iterative estimator, we set the simulation parameters as follows: the signal magnitude , the sampling frequency kHz, the frequency offset , the spectral peak = 64 kHz, the true frequency kHz < < 64.5 kHz, the number of iterations [24], and the initial phase 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.

Figure 6

The RMSE/CRLB of the proposed estimator versus under different settings of M and Q when dB and .

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 and , there are two periods for and and one period for and . They are both symmetric about . 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 and simulate seven estimators.

Figure 7, Figure 8, Figure 9 and Figure 10 are figures corresponding to , , , and , respectively. Upon comparing these figures, the SNR threshold decreases from 0 dB to dB. The RMSEs of these estimators and the RMSE differences among these estimators decrease as N increases. Detailed figures at 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 times the CRLB when dB and .

Figure 7

RMSE comparison among frequency estimates versus SNR at and .

Figure 8

RMSE comparison among frequency estimates versus SNR at and .

Figure 9

RMSE comparison among frequency estimates versus SNR at and .

Figure 10

RMSE comparison among frequency estimates versus SNR at and .

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 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 . For some values of , the corresponding RMSEs may lie slightly below the CRLB and are asymmetric relative to due to the limited number of Monte Carlo trials.

Figure 11

RMSE comparison among the frequency estimates versus at and dB.

Figure 12

RMSE comparison among the frequency estimates versus at and dB.

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 as the DFT length and compute the FFT instead. An N-point FFT requires complex multiplications and complex additions. A one-bin N-point DTFT requires N complex multiplications and 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.

Table 2

Computational complexities of the estimators used in the simulation.

Estimators	Complex Multiplications	Complex Additions
Jacobsen (N points)	N2log_2N+1	Nlog_2N+3
Candan (N points)	N2log_2N+1	Nlog_2N+3
A&M (N points, 2 iterations)	N2log_2N+4N+2	Nlog_2N+4N
Fan (M points, 2 iterations)	M2log_2M+5N+8	Mlog_2M+5N+1
Fang (M points)	M2log_2M	Mlog_2M
RCSTL (M points)	M2log_2M+1	Mlog_2M
Proposed (M points, 2 iterations)	M2log_2M+5N	Mlog_2M+5(N−1)

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.