Robust and Efficient Frequency Estimator for Undersampled Waveforms Based on Frequency Offset Recognition

This paper proposes an efficient frequency estimator based on Chinese Remainder Theorem for undersampled waveforms. Due to the emphasis on frequency offset recognition (i.e., frequency shift and compensation) of small-point DFT remainders, compared to estimators using large-point DFT remainders, it can achieve higher noise robustness in low signal-to-noise ratio (SNR) cases and higher accuracy in high SNR cases. Numerical results show that, by incorporating a remainder screening method and the Tsui spectrum corrector, the proposed estimator not only lowers the SNR threshold of detection, but also provides a higher accuracy than the large-point DFT estimator when the DFT size decreases to 1/90 of the latter case.


Introduction
Frequency estimation of undersampled waveforms is widely encountered in radar detection [1] and distance estimation [2] etc. The mainstream approach to solve this problem is Chinese Remainder Theorem (CRT), since it is able to reconstruct a large number from its remainders modulo several moduli. For the CRT-based frequency estimators [3][4][5][6], the remainders are directly acquired from DFT peak bins of multiple channels, which means that the estimation accuracy heavily relies on the DFT frequency resolution. Hence, once the undersampling rates of individual channels are determined, to obtain a high estimation accuracy, the existing estimators [4][5][6] use large-point DFT to enhance the DFT frequency resolutions, which inevitably results in high computational complexity. Therefore, an efficient CRT-based frequency estimator with small-point DFT needs to be developed.
Unlike the CRT-based estimators with large-point DFT, an estimator with small-point DFT tends to incur weak noise robustness and low accuracy. On one hand, as is known, the noise energy is uniformly distributed over all the DFT bins. As the number of DFT bins decreases, the noise contamination on each DFT bin increases, thereby the noise robustness of the estimator is deteriorated. On the other hand, small-point DFT brings a coarse frequency resolution and thus degrades the accuracy.
To overcome these two difficulties, this paper proposes a scheme based on frequency offset recognition, which was ignored by the existed estimators [4][5][6]. The scheme is based on a discovery that, for any reconstruction channel, the underlying frequency offset information can be utilized to highlight the peak DFT bin.
Specifically, in low SNR cases, for any individual undersampled waveform, we use a performance index to recognize the value range of the frequency offset and then implement frequency-shift operation to highlight the peak DFT bin buried in noise. Moreover, we also employ a measure of remainder screening [6,7] to further recognize those channels with reliable DFT spectra. In high SNR cases, we employ the Tsui spectrum corrector [8] to accurately estimate the frequency offset, thereby to overcome the drawback of the coarse DFT frequency resolution. Numerical results demonstrate that, with the above frequency-offset related measures, the proposed estimator can provide better noise robustness and higher accuracy than the existing estimators whose DFT lengths are about 90 times larger than ours.
This paper is organized as follows. Firstly, we present the estimation model and make the error analysis of the existed long-point DFT based estimators. Secondly, we propose a frequency estimator based on frequency shift and compensation, spectrum correction and remainder screening. Lastly, we present the simulation comparison which verifies the proposed estimator's superiority over the existed estimators.

CRT-based Frequency Model
Consider a signal formulated as where f 0 is the high frequency to be determined and a is a non-zero complex amplitude.
To estimate f 0 , L A/D converters with undersampling rates F 1 * F L (assume F 1 < F 2 < Á Á Á < F L ) are used to discrete x(t) in parallel. To meet the requirement of the closed-form CRT [9,10], F 1 * F L need to be integers such that where M is the greatest common divider (GCD) of F 1 * F L and thus integers Γ 1 * Γ L are pairwise coprime. Then, the i-th (i = 1, . . ., L) discrete sequence with length M can be denoted as Therefore, the relationship between f 0 and the individual undersampling rates F 1 * F L is formulated as where n i is the unknown folding integer and η i is the normalized frequency (0 η i < 1). Clearly, Eq (4) constitutes a CRT model, in which F 1 * F L are moduli and η 1 F 1 * η L F L are remainders. In other words, once the normalized remainders η 1 * η L are acquired, the closed-form CRT [9,10] can determine the folding integers n 1 * n L and thus yields the estimate of f 0 .

Error Analysis of the Existed Estimators
To acquire η i , the estimators in [4][5][6] implement F i -point DFT (zeros are padded to the end of sequence) on x i (n), i = 1, . . ., L, and calculate η i by searching out the peak DFT bin (assume it falls at k = k i ), i.e., In [4][5][6], f 0 is assumed to be an integer (i.e., integer times of the frequency unit Δf = 1 Hz), which ensures that However, in practice, f 0 also allows to be a real number (i.e., its fractional part is nonzero). In this case, η i − k i /F i 6 ¼ 0, i.e., estimate error inevitably occurs. Furthermore, the DFT size F i = MΓ i is really a bit large, implying that high computational complexity is consumed.

Problems of A Small-point DFT Based Estimator
To reduce the complexity, we attempt to substitute MΓ i -point DFT with M-point DFT. Accordingly, the frequency unit Δf increases from 1 Hz to F i /M = Γ i Hz, indicating that the i-th individual DFT spectrum gets much coarser. To ensure a sufficiently high accuracy, apart from the peak bin index k = k i , the fractional part of the remainder η i F i should also be taken into account. In this case, we have The fractional number δ i in Eq (6) refers to a normalized frequency offset, which was ignored by estimators in [4][5][6]. Clearly, estimation of δ i can make up for the deficiency of coarse resolution of small-point DFT.
The other difficulty is the noise robustness. Specifically, for a same low SNR circumstance, compared to large-point DFT, the true peak bin of small-point DFT is more likely to be buried by heavy noise. As Eq (6) shows, compared to the estimate of δ i , the searching error of the integer peak index k i will result in a larger estimation error of η i . In particular, the searching error of k i might result in an incorrect estimate of the folding number n i , thus leading the CRT reconstruction into failure.

Frequency Shift and Compensation
To improve the estimator's robustness to noise, it is necessary to enhance the magnitude of peak DFT bin of x i (n) as possible.
Combining Eqs (3), (4) with (6), one can deduce x i (n) as Then, the normalized peak DFT bin X i (k i ) can be derived as Eq (8) shows that, the magnitude of peak DFT bin is closely relevant to the frequency offset δ i . Fig 1 gives the curve of |X i (k i )/a| versus δ i .
From Fig 1, one can observe that as |δ i | increases, |X i (k i )/a| monotonously decreases, indicating that the DFT spectral leakage becomes increasingly more serious. As a result, for the low SNR case, it is likely that the peak DFT bin is buried.
To obtain a frequency offset approximating 0 as possilbe, we uniformly divide the value range À Â Á , i.e., the middle region nearby 0. In terms of the property of Fourier transform, x i (n) needs to be applied with a left frequencyshift operation as 2) Case of d i 2 À 1 6 ; 1 6 Â Á : Since δ i ! 0 (meaning that the DFT spectral leakage is not serious), frequency-shift operation is not necessary, i.e., x M i ðnÞ ¼ x i ðnÞ.
However, the problem lies in judging which region δ i falls in. Given the DFT spectra X L i ðkÞ; X M i ðkÞ; X R i ðkÞ of the above shifted sequences, we present a performance index defined by the ratio between the peak DFT bin and its adjacent sub-highest DFT bin, i.e., where c 2 {L, M, R}. Further, the region of δ i can be recognized from the maximum of α L , α M , α R . Here, we present an example to illustrate this recognition criterion. One can observe that, among 3 spectra, |X L (k)| in Fig 2(a) exhibits the smallest spectral leakage, i.e., its peak at k = 6 appears more conspicuous than others. Quantitatively, in terms of Eq (12), we have α L = 13.9805, α M = 1.4995, α R = 2.7479. Thus, we can recognize that d 2 1 6 ; 1 2 Â Á , which is in accord with the fact that δ = 0.4.
Further, in noisy circumstances, while the noise energy uniformly covers all DFT bins, the proposed frequency-shift operation still tends to highlight the peak DFT bin and thus yields a higher noise robustness.

Spectrum Correction and Remainder Screening
The above recognition criterion can provide the reliable peak index k i and the value range of the offset δ i . To improve the accuracy, δ i needs to be further estimated. A lot of spectrum correctors such as Macleod correctors [11], Candan Correctors [12], Phase-difference correctors [13], Tsui correctors [8] are competent to estimate δ i from the DFT result of x i (n). This paper suggests employing the Tsui spectrum corrector proposed in [8] (Its dataflow will be illustrated in the next section), since it currently possesses a high noise robustness and accuracy. Following this, all the CRT remaindersr i ¼ Z i F i ¼ ðk i þ d i ÞF i =M; i ¼ 1; :::; L, can be acquired.
Further, [5,9] claimed that, the necessary condition of successful CRT reconstruction is where r i is the ideal errorless remainder. However, if the SNR is too low, not all L remainders satisfy this error bound condition.
[6] provides a remainder screening method. Specifically, if there are b(L−K)/2c or fewer remainders exceeding the error bound Eq (13), this method can pick out these unrestricted remainders and use the remained remainders to achieve a successful CRT reconstruction. Accordingly, the CRT reconstruction range decreases from f 0 2 ½0; M This procedure of remainder screening will be listed in the next section.

Procedure of the Proposed Estimator
The proposed frequency estimator based on frequency offset recognition is summarized as follows: Step 1 Implement 3 frequency-shift operations on the undersampled sequence x i (n), i = 1, Á Á Á, L, to generate the sequences x L i ðnÞ; x M i ðnÞ; x R i ðnÞ and their DFT results X L i ðkÞ, X M i ðkÞ, X R i ðkÞ. Further, calculate their performance indices a L i ; a M i ; a R i in terms of Eq (12).
Step 2 Use the frequency-offset based criterion to recognize the expected X c i ðkÞ; c 2 fL; M; Rg. Then, employ the Tsui spectrum corrector to obtain its normalized frequency estimate Z c i . Further, use Eq (11) to obtain the compensated result η i . Thus, we can obtain all the remaindersr i ¼ Z i F i ; i ¼ 1; :::; L.
Step 3 Substitute the moduli F 1 ,. . .,F L , the remaindersr 1 ; :::;r L and the specified integer K into the remainder-screening based CRT reconstruction procedure to obtain the final frequency estimatef 0 .
The procedure of the Tsui spectrum corrector addressed in Step 2 is illustrated in Fig 3. The procedure of the remainder-screening based CRT reconstruction is listed in the following: 1. Employr 1 as the reference remainder and use the CRT remaindersr 1 ; :::;r L to calculate the following L − 1 difference remaindersq i;1 : 2. Calculate the remainders ofq i;1 " G i;1 modulo Γ i : where " G i;1 is the modular multiplicative inverse of Γ 1 modulo Γ i and can be calculated in advance.
3. Compute an intermediate integer X 1 using the following formula where Γ ≜ Γ 1 Γ 2 . . .Γ L , and W i, 1 is the modular multiplicative inverse of Γ/(Γ 1 Γ i ) modulo Γ i that can be calculated in advance.

4.
To determine the folding integer n 1 : Construct a set Z ¼ f If there is only one z 0 in Z satisfying 0 X ðz 0 Þ < Q K i¼1 G i , let folding number n 1 ¼Xðz 0 Þ.

5.
To determine other folding integers n j , 2 j L: Exchange the reference remainderr 1 witĥ r j and repeat 1)-4) to acquire the corresponding folding integersn j respectively. It is possible that an individual n j is null in an exchanged operation. Assume that altogether P (P L) times of operations are successful and their exchanged remainder indices are v 1 , . . ., v P . Thus, P frequency estimatesf 0 ðn i Þ ¼n n i F n i þr n i for 1 i P can be obtained.

Numerical Results
In this section, we compare the noise robustness and the accuracy among the large-point-DFT CRT estimator in [5] (LP-DFT estimator) and other 3 small-point-DFT CRT estimators: Tsui correctors with frequency shift (Tsui-FC estimator), Tsui correctors without frequency shift (Tsui estimator), and the proposed estimator. In this way, the improvement contribution of the frequency-shift operation, the spectrum correction and the remainder screening can be quantitatively investigated.
Example 2: SNR varies in a low region [−23dB, −8dB]. For each SNR case and each estimator, 5000 Monte Carlo trials were conducted. We use the detection probability P d to evaluate the noise robustness. Its criterion is as follows: For each trial, iff 0 satisfies jf 0 Àf 0 j < f 0 =1000, then this trial passes; otherwise, it fails. Fig 4 illustrates the P d curves of these 4 estimators.
From Fig 4, the following conclusions can be drawn: 1. The proposed estimator (red circle curve) shows the best noise robustness, and the LP-DFT estimator (black diamond curve), the Tsui-FC estimator (green triangle curve), the Tsui estimator (blue rectangle curve) take the 2nd, 3rd, 4th place respectively, since their P d curves successively appear from left to right.

Conclusions
This paper proposes a small-point DFT based frequency estimator for undersampled waveforms, in which the recognition of the frequency offset δ i plays a critical role in both high SNR cases and low SNR cases. Numerical results verifies the proposed estimators superiority in noise robustness, accuracy, and efficiency over the large-point DFT based estimator, which presents a vast potential for future applications.