Multi-Tone Frequency Estimation Based on the All-Phase Discrete Fourier Transform and Chinese Remainder Theorem

The closed-form robust Chinese Remainder Theorem (CRT) is a powerful approach to achieve single-frequency estimation from noisy undersampled waveforms. However, the difficulty of CRT-based methods’ extension into the multi-tone case lies in the fact it is complicated to explore the mapping relationship between an individual tone and its corresponding remainders. This work deals with this intractable issue by means of decomposing the desired multi-tone estimator into several single-tone estimators. Firstly, high-accuracy harmonic remainders are calculated by applying all-phase Discrete Fourier Transform (apDFT) and spectrum correction operations on the undersampled waveforms. Secondly, the aforementioned mapping relationship is built up by a novel frequency classifier which fully captures the amplitude and phase features of remainders. Finally, the frequencies are estimated one by one through directly applying the closed-form robust CRT into these remainder groups. Due to all the components (including closed-form CRT, the apDFT, the spectrum corrector and the remainder classifier) only involving slight computation complexity, the proposed scheme is of high efficiency and consumes low hardware cost. Moreover, numeral results also show that the proposed method possesses high accuracy.


Introduction
Frequency measurement is a fundamental problem in signal processing, which is widely encountered in instrumentation, digital communication, radar, etc. Numerous frequency estimation methods have been proposed, including the classical discrete Fourier transform (DFT) [1,2], the subspace-based methods [3,4], maximum likelihood [5,6], linear or nonlinear least squares [7,8]. However, as some applications work in a wider band and higher frequency, e.g., the millimeter-wave band in 5G technologies, these methods become impractical, since the realizable sampling rates of the analog to digital converters (ADCs) are limited by the Nyquist theorem. Therefore, investigations on the frequency estimation from undersampled sequences are interesting.
The Chinese Remainder Theorem (CRT) is an effective approach to resolve ambiguity related problems including the undersampling frequency estimation [9][10][11][12][13][14]. The basic idea is to reconstruct a larger number M from its residue set {r l , l = 1, ..., L} modulo multiple moduli {M l , l = 1, ..., L}. Concerning the frequency estimation problem, L ADCs with sub-Nyquist sampling rates {M l , l = 1, 2..., L} are employed to obtain the undersampled waveforms, and then the classical DFT is performed on the acquired waveforms to detect L ambiguous frequencies {r l , l = 1, 2, ..., L}. Moreover, compared to some searching based estimator for undersampled waveforms [9,11], the CRT-based estimator [13,15] can achieve a large reconstruction upper bound, which equals the least common multiple (lcm) of the undersampling rates {M l }.
Recently, some modified versions of the CRT reconstruction algorithms have been presented for computation complexity reduction and robustness enhancement [10,11,13,14,16]. Specifically, the reduced complexity owns to the fact that CRT can work in a closed-form way rather than in a searching-based way. Besides, the robustness enhancement lies in the fact that CRT can also acquire a high reconstruction accuracy when all the remainders are erroneous, as long as these remainder errors do not exceed a quarter of the greatest common divisor (gcd) of all the moduli. Particularly, the remainders of the existing CRT-based estimators are derived from the DFT results of undersampled waveforms. Hence, if we can resort to another high-performance spectral analysis tool to replace DFT, both accuracy and anti-noise robustness will be further enhanced.
As is well known, DFT spectral analysis has two drawbacks: heavy spectrum leakage arising from data truncation and insufficient spectral resolution incurred by picket fence effect. To suppress the spectral leakage, we proposed the all-phase DFT (apDFT) spectral analysis in [17] and pointed out that apDFT spectrum's sidelobe leakage is much slighter than DFT even when dealing with multi-tone waveforms. Moreover, as [18] pointed out, if apDFT is combined with the technique of spectrum correction, the spectral resolution can be improved significantly. This combination actually provides us a good idea to obtain higher-accuracy remainders required by CRT reconstruction.
In recent years, many endeavors have been made to generalize CRT-based estimators to multi-tone undersampled waveform cases [19][20][21]. Generally, these estimators solve this problem through the remainder redundancy coding, which actually pays the cost of heavy computation burden and sacrificing the dynamic estimation range. To break the dilemma, we develop a novel estimator combining closed-form CRT, apDFT and spectrum correction. Different from the existing estimators, our proposed estimator's applicability in the multi-tone case is ensured by building up a mapping relationship between an individual tone and its corresponding remainders. Specifically, this mapping relationship is realized by a novel harmonic-parameter clustering method, which is closely related to apDFT and spectrum correction. With the above considerations, the proposed multi-tone estimator can be converted into multiple single-tone estimators, and thus individual tones can be retrieved one by one. Numerical results show that, our proposed estimator concurrently possesses high accuracy and large dynamic range. Moreover, our proposed estimator is applicable to the case of only 2 undersampling paths, whereas the existing multi-tone estimators cannot apply to this case.
The remaining of this paper is organized as follows. Problem Formulation of CRT-based frequency estimation is given in Section 2. In Section 3, we detail the remainder acquisition and remainder classification based on the apDFT and harmonic-parameter clustering in the proposed method. The numeral results and performance analysis are presented in Section 4. Finally, the conclusions are drawn in Section 5.

CRT Reconstruction Model for the Single-Tone Case
In this subsection, the CRT reconstruction model for the single-tone case is formulated. A narrow-band signal x(t) with a single tone is formulated as where a, θ 0 and f 0 are the amplitude, initial phase and the frequency to be determined, respectively. w(t) is the additive white Gaussian noise with zero mean and variance σ 2 .
The Nyquist sampling theorem requires that the sampling rate F l must be at least as twice as the signal frequency f 0 to avoid the ambiguity. Definitely, the ambiguity occurs in the undersampling case, i.e., F l < f 0 . In this case, the detected frequency r l can be represented as where n l is the folding integer, and r l is the ambiguous frequency which can be acquired by performing the traditional DFT on the undersampled waveforms. Equation (2) is fully in accordance with the model of closed-form robust CRT [13], providing a basis for the undersampling frequency estimation. Guided by the determination procedure in [13], x(t) needs to be discretized with L undersampling rates (also acting as the L moduli of CRT) F 1 , ..., F L f 0 . Typically, denote N be their greatest common divisor (i.e., the largest integer that divides each of them), namely N = gcd{F 1 , · · · , F L }, thus integers Γ 1 , · · · , Γ L valued with constitute a co-prime integer set. As such, the moduli F l can be decided by the coprime integer set {Γ l } and one specified integer N. Then, L undersampled versions of x(t) are generated as Since f 0 F 1 , ..., F L , a simultaneous congruence equation set can be built up as where n 1 , ..., n L are unknown folding integers and ε 1 , ..., ε L are normalized frequencies of L undersampled sequences. As mentioned above, the ambiguous frequencies {ε l F l } (the remainders) can be approximated via the traditional DFT. Then, the frequency estimatef 0 can be acquired by feeding the remainders {ε l F l } and moduli {F l } into the CRT reconstruction algorithm in [13]. The analysis above gives a simple review of frequency estimation in the framework of CRT reconstruction. However, there are some open questions associated with the above method. Firstly, the robust CRT usually requires that the error in each remainder is less than one quarter of the gcd of the moduli. Under this condition, the reconstruction error can be upper bounded by the same range as that of the remainders. For another, the remainder acquisition is always achieved by the traditional N-point DFT, in which the normalized frequency resolution is 1/N. Accordingly, the normalized frequency ε l in (5) can be represented as where k l ∈ {0, 1, ..., N − 1} and δ l is a fractional number ranging in the interval [−0.5, 0.5). Hence, the remainder r l can be rewritten as The traditional DFT allows us to obtain the integer k l , whereas the fractional number δ l (also referring to the frequency offset) tends to be erroneous due to the spectrum leakage and the picket-fence effect in the traditional DFT analysis. Therefore, the errors arising from the DFT analysis tool inevitably give rise to the errors in the reconstruction results. In this sense, it is meaningful to resort to another high-performance spectral analysis tool to replace DFT.

CRT Reconstruction Model for the Multi-Tone Case
For a multi-tone signal x(t) formulated as Assuming that L undersampling rates are the same as the single-tone case, L discretized versions of x(t) can be represented as Accordingly, for these tones, ρ simultaneous congruence equation sets can be built up as where ε l,m F l refers to the required CRT remainder r m,l for the m-th tone at the l-th reconstruction path, i.e., where k l,m ∈ {0, ..., N − 1} and δ l,m ∈ [−0.5, 0.5). For any index l ∈ {1, ..., L}, substituting (10) into (9) yields a m e j(2πε l,m n+θ m ) , n = 0, ..., N − 1.
To emphasize, the above multi-tone remainder acquisition is totally distinct from that of the single-tone case. In the single-tone case, there exists one peak bin of the DFT spectrum at any reconstruction path, consequently, the estimates of remainders can be obtained through collecting all the peak bins. However, in the multi-tone case, the DFT spectrum at any reconstruction path surely contains ρ peak bins and the mapping relationship between the peak bins across different paths is unknown.
We give a simple example to illustrate the unknown relationship in the multi-tone case. We assume the three numbers are {5, 23, 181}, and 181 and two moduli are {7, 9}. In this case, the two remainder sets which can be detected are {2, 5, 6} and {1, 5}, respectively. Considering the second remainder set {1, 5}, we cannot tell which element repeats twice. As for element 5, it is unclear which remainder in the first set corresponds to it. Hence, the difficulty of the CRT algorithm for multiple numbers lies in building up the mapping relationship between each number and the corresponding remainders. Specific to the frequency estimation in the multi-tone case, the difficulty lies in categorizing ρL peak bins across L DFT spectra into ρ remainder classes. In this way, the multi-tone frequency estimators can be decomposed into several single-tone estimators.

All-Phase DFT Based Remainder Acquisition
In the proposed method, we combine the apDFT and spectrum correction to achieve the remainder acquisition. The combination can restrain the spectrum leakage and mitigate the fence effect in the traditional DFT, thereby improving the remainder accuracy.
In the traditional DFT, N samples should be collected for the N-point DFT, as shown in (4) and (9). Unlike the N-length sampling mechanism, a (2N − 1)-length sequence is required in the apDFT, from which a new N-length sequence can be derived for N-point DFT.
the key idea of the apDFT is to derive a new N-length sequence from the sampled (2N − 1)-length sequence. Without loss of generality, the multi-tone case is considered to illustrate the windowed apDFT.
Firstly, the N-length sequence x l (n) in (9) should be expanded to the (2N − 1)-length one, Given a N-length window function w 1 (n), e.g., Hanning window function, N different sequences x l,q (n), q = 0, ..., N − 1 with overlapped each other can be derived from the raw (2N − 1)-length sequence x l (n). This transformation can be achieved by multiplying x l (n) with the shifted window function w 1 (n + q), i.e., where w 1 (n + q) stands for shifting w 1 (n) by q to the left. It is natural to obtain a new N-length sequence x la (n) by simply averaging the corresponding elements in the sequence set {x l,q (n), q = 0, ..., N − 1}. Alternately, it can be done by weighed averaging. Typically, given another N-length window function w 2 (n), the weighted sequences {y l,q (n), q, n = 0, ..., N − 1} can be formulated as follows where the integer r is utilized to perform the N-point cyclic shift of x l,q (n). Furthermore, then the N-length sequence x la (n) can be derived as The analysis above studies N sequences {x l,q (n), q = 0, ..., N − 1} to determine the unique N-length sequence x la (n). Usually, the process is referred to as the all-phase preprocessing, accordingly the signal x la (n) is called the all-phase signal. Similarly, the all-phase DFT X la (k) can be calculated by performing the N-point DFT on x la (n).
Furthermore, for a given sample point, such as x l (0) in (13), all sequences containing x l (0) are derived in (14). On the contrary, the traditional DFT considers one case q = 0 only. Hence, the DFT spectrum based on the all-phase signal x la (n) can estimate the frequency components with smaller spectrum leakage that arises from the data truncation.
From (14), it is possible to notice that x(0) exits all the possible points in the N-length sequence, i.e., all the possible phase, so it is referred to as the all-phase prepossessing. This also leads to the phase invariance in all-phase DFT.
In order to illustrate the superiority clearly, consider a multi-source exponential signal x(n) = ∑ ρ l=0 e j(w l n+φ l ) , where N = 64, w l = β l · 2π N . Assume that x(n) consists of three frequency components with β 1 = 12, β 2 = 22.2, β 3 = 28.4 respectively. The initial phases of the three components are set as φ 1 = 100 • , φ 2 = 50 • , φ 3 = 80 • . The amplitude spectra and phase spectra for N-point windowed DFT and double-windowed all-phase DFT are shown in Figures 1 and 2  From Figure 1, we observe that the amplitude spectrum for windowed all-phase DFT has clearer peak bins and smaller side bins, verifying that the all-phase DFT can effectively restrain spectrum leakage.
From Figure 2, the phase estimates can be directly obtained around the corresponding peak bins in the phase spectrum of apDFT in Figure 2b. Especially, the phase spectrum of apDFT does not change with k, which differs from that of traditional DFT in Figure 2a.
Through the analysis above, all-phase DFT outperforms the traditional DFT since that the preprocessing on the (2N − 1)-length sequence can significantly restrain the spectral leakage arising from the data truncation. Moreover, to reduce the errors incurred by picket effect, it is vital to adopt some correction methods to obtain accurate harmonic parameters from finite spectrum lines.
Ref. [22] pointed out that, if an exponential sequence is implemented with the Hanning windowed DFT, the ratio between the peak DFT bin |X l (k l,m )| and its sub-peak neighbor contains the information of the frequency offset. Specifically, this amplitude ratio v can be calculated as v = |X l (k l,m )| max{|X l (k l,m − 1)|, |X l (k l,m + 1)|} .
In order to apply the ratio-based spectrum correction to the all-phase DFT, the amplitude ratio v a in our paper is specified as v a = |X a,l (k l,m )| max{|X a,l (k l,m − 1)|, |X a,l (k l,m + 1)|} .
On basis of [22], we can deduce the frequency offset estimateδ l,m aŝ Accordingly, the normalized frequency estimate is derived aŝ Then, the corrected amplitude is obtained by substitutingε l,m into the exponential term in (12) a l,m = πδ l,m (1 −δ 2 l,m ) · |X l,m (k l,m )|/ sin(πδ l,m ), As mentioned earlier, we can directly extract the initial phase information from the phase spectrum of all-phase DFT, since that the all-phase DFT adopted in the proposed method has the excellent performance of the initial phase invariance. That means only the amplitude and frequency need be corrected via (20) and (21).

Harmonic Parameter Clustering Based Remainder Classification
As mentioned in Section 2.2, the main difficulty in multi-tone case lies in developing the mapping relationship between {r l,m } and { f m }. Therefore, we present a harmonic parameter clustering to solve this problem.
It can be inferred from (12) that, due to multi-path undersamplings, the resultant normalized frequencies ε 1,m , ..., ε L,m of the m-th tone surely differ with each other. Nevertheless, they share a common harmonic parameter pair {a m , θ m }, which provides the basis for the remainder classification.
Note that the harmonic-parameter triple {ε l,m ,â l,m ,θ l,m } actually is bound together for the m-th tone at the l-th reconstruction path. To construct patterns convenient for further remainder clustering, the following vector quantities need to be built up as: z l,m =â l,m e jθ l,m , l = 1, ..., L, m = 1, ..., ρ.
For any individual pattern of the m-th tone at the l-th path, we might as well select the pattern z 1,m of the 1-st path as the reference. Then, its clustering indicators c l,m can be determined by finding the closest distance among ρ patternsẑ l,1 , ...,ẑ l,ρ , i.e., In this way, altogether ρL indicators c l,m , m = 1, ..., ρ, l = 1, ..., L, can be collected, besides that the reference indicators c 1,m = 1, m = 1, ..., ρ.
Particularly, for the two-path case (i.e., L = 2), two indicator-involved remainders r 1,m , r 2,c 2,m can be determined via the above harmonic-parameter clustering operations. Following this, feeding the moduli F 1 , F 2 and r 1,m , r 2,c 2,m into the procedure of the closed-form robust CRT in [13] yields the final estimatef m , m = 1, ..., ρ.

Determination Procedure of Multi-Tone Frequency Based on apDFT Analysis and CRT
The proposed method for multi-tone frequency estimation can be summarized as follows (Algorithm 1).

Simulation Results and Discussion
In this section, the simulation is carried in MATLAB R2016b, with an Intel Core i5 2.60 GHz. To emphasize, the existing CRT-based multi-tone frequency determination approaches [19][20][21] cannot apply to the case that the reconstruction path number L is smaller than the component number ρ. Therefore, this section will first verify the feasibility of the case L < ρ, meaning that there are more frequency components to be estimated than the data acquisition paths. As Ref. [22] pointed out, the Hanning Window is adopted to guarantee the applicability of spectrum correction.

Procedure Demonstration
For the case of L = 2 and ρ = 4, the corresponding frequencies, amplitudes and initial phases of the four components are listed in Table 1. To validate the proposed method in high-frequency scenarios, four tones with high frequencies (up to the GHz level) are considered. The input parameters are set as follows: the co-prime integers Γ 1 = 3301, Γ 2 = 3307, the gcd N = 512. In terms of (3), two ADC sampling rates F 1 = NΓ 1 = 1.690112 × 10 6 samples/s and F 2 = NΓ 2 = 1.693184 × 10 6 samples/s, much lower than the signal frequencies f 1 , ..., f 4 listed in Table 1. The noise w(t) in is additive white Gaussian noise with mean zero and variance σ 2 n . In this paper, the SNR is defined as where σ 2 s indicates the signal mean power. In this simulation, the SNR is specifies as 13 dB. Through undersampling and all-phase preprocessing (Hanning double-window is adopted) in Step 1, Step 2, two N-length all-phase sequences x 1a (n), x 2a (n) can be generated. Following Step 3, two apDFT spectra X 1a (k), X 2a (k) are acquired and illustrated in Figure 3

Following
Step 4, through applying ratio-based spectrum correction on these peak DFT bins, the ρL = 8 harmonic parameter triples are determined. With that, the vector quantities and remainder information {ẑ l,m ,ε l,m }, l = 1, 2, m = 1, ..., 4 can be generated, which are listed in Table 2. As a comparison, the uncorrected vector quantities (generated from the DFT spectra directly) and these corrected vector quantitiesẑ l,m are illustrated in Figures 4 and 5, respectively. It becomes clear that the ρL = 8 vector quantities in Figure 5 can be intuitively divided into four sets, which correspond to the four tones, respectively, whereas the uncorrected patterns in Figure 4 appear chaotic.  Figure 5 depicts. Further, in terms of r l,m =ε l,c l,m F l , l = 1, 2, ρ = 4 remainder sets {r 1,m ,r 2,m } are acquired and listed in Table 3. Step 6, successively substituting these 4 remainder sets, the gcd N = 512 and the moduli {F 1 , F 2 } into the closed-form CRT reconstruction procedure [13] yields the final frequency estimatesf 1 , ...,f 4 listed in Table 4, which reflect that the errors are almost negligible.
In this simulation, we present the step-by-step frequency determination procedure based on the given parameters. For the case of more frequency components than the data acquisition paths, the effectiveness of the proposed method has been verified. In addition, the results also verify the applicability to the high-frequency scenario. we conduct the experiment to verify the superiority of the apDFT over the traditional DFT in the remainder acquisition. Based on the same multi-tone signal given in Section 4.1, the estimation accuracy of the proposed approach across different tones versus SNRs ranging from 6 dB to 23 dB is studied. Herein, the significance of SNR is identical with that in Section 4.1. The root-mean-square error (RMSE) where the superscript i refers to the i-th trial and Q denotes the number of Monte Carlo tests. For each SNR case, 500 Monte Carlo trials were conducted. The RMSE curves of four tones are shown in Figure 6. To ensure the comparability, another experiment was conducted under the same condition except that the remainder acquisition was achieved by the traditional DFT. Similarly, the RMSE results of four tones are shown in Figure 7.  As Figures 6 and 7 depict that, for any tone, there is a SNR threshold under which the estimation error abruptly increases. For the region where the SNR is above the threshold, the relative errors corresponding to any tone are below the level 100/10 9 × 100% = 0.00001%. The result matches perfectly with that of the single trial in Section 4.1.
As a comparison, for any tone, the threshold value of the proposed method (utilizing the all-phase DFT) in Figure 6 is smaller than that utilizing the traditional DFT in Figure 7. The simulation confirms that the apDFT can result in the improvement of anti-noise robustness compared with the traditional DFT. The improvement in anti-noise is due to the property of restraining spectral leakage of apDFT, which can detect remainders with higher accuracy. Essentially, it can be attributed to the all phase preprocessing mechanism.

Performance Analysis on Different Number of Data Acquisition Path
In this simulation, some simulation results are presented to investigate the estimation performance in different data acquisition number. For simplicity and effectiveness, the overall RMSE of four tones is used to assess the frequency estimation performance. Herein, the method utilizing the traditional DFT (as mentioned in Section 4.2) is also considered as a reference. Figure 8 illustrates the RMSE results under the channel number L = 2, 3, 4 based on the model parameterization in Section 4.1. For each channel number, the apDFT-based method has a better performance. This corresponds to the results in Figures 6 and 7. For another, as to the proposed method, the more data acquisition paths, the smaller the value of threshold is. To some extent, the robustness can be enhanced by increasing channel numbers, whereas the enhancement cannot be ensured when the peak bins are contaminated by the heavy noise. In addition, it is interesting to observe that, in the low SNR region, the error magnitude associated with more data acquisition paths like the channel number L = 4 is significantly higher. This occurs when the clustering results are invalid due to heavy noise contamination, and then more cumulative errors are generated with the increasing channel number.

Analysis of Computation Complexity
In the proposed method, the multi-tone frequency estimation is achieved by incorporating the apDFT, the harmonic-parameter clustering and the closed-form CRT reconstruction algorithm. Note that the CRT reconstruction algorithm adopted works in a closed-form way. Furthermore, the spectrum correction and the harmonic-parameter clustering only involve in simple algebra calculations. Therefore, the computation complexity of the proposed estimator mainly depends on the N-point apDFT in L data acquisition paths, i.e., O(N/2log 2 N). Compared with the existing CRT-based estimators, the proposed method works well even when the data acquisition path number equals two and no searching step is required. Obviously, no heavy computation burden is required in the proposed method.

Conclusions
In this paper, we propose a multi-tone frequency estimator from undersampled waveforms, which incorporates the closed-form CRT, the apDFT spectral analysis, the spectrum correction, and the harmonic-parameter clustering. This organic technique combination allows that the multi-tone estimator is decomposed into multiple single-tone estimators. Thus these tones can be recovered one by one. Two remarkable merits should be emphasized.
On one hand, different from the existing CRT-based multi-tone estimators with remainder redundancy coding, the proposed estimator can still work well even if the reconstruction path number L = 2, which actually greatly reduces the hardware cost. On the other hand, due to the utilization of apDFT, which yields excellent suppressing effect of spectral leakage and noise, our proposed estimator can acquire a higher anti-noise robustness than the DFT case.
The above two merits give the proposed estimator with vast potentials in high-frequency measurement related applications such as radar, future beyond 5G communications.