Next Article in Journal
Coherence Time Evaluation in Indoor Optical Wireless Communication Channels
Previous Article in Journal
Ultrawideband, Wide Scanning Stripline-Fed Tightly Coupled Array Antenna Based on Parallel-Dipole Elements
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

School of Electrical and Information Engineering, Tianjin University, Tianjin 300072, China
*
Author to whom correspondence should be addressed.
Sensors 2020, 20(18), 5066; https://doi.org/10.3390/s20185066
Submission received: 8 July 2020 / Revised: 29 August 2020 / Accepted: 31 August 2020 / Published: 7 September 2020
(This article belongs to the Section Communications)

Abstract

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

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

2. Problem Formulation

2.1. 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
x ( t ) = a e j ( 2 π f 0 t + θ 0 ) + w ( t ) ,
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
f 0 = n l F l + r l , l = 1 , . . . , L
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
Γ l = F l / N , l = 1 , . . . , L ,
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
x l ( n ) = a e j 2 π f 0 F l n + θ 0 + w ( n ) , n = 0 , . . . , N 1 .
Since f 0 F 1 , . . . , F L , a simultaneous congruence equation set can be built up as
f 0 = n 1 F 1 + ε 1 F 1 f 0 = n L F L + ε L F L ,
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 estimate f ^ 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
ε l = ( k l + δ l ) / N , l = 1 , . . . , L ,
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
r l = ε l F l = ( k l + δ l ) F l / N , l = 1 , . . . , L .
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.

2.2. CRT Reconstruction Model for the Multi-Tone Case

For a multi-tone signal x ( t ) formulated as
x ( t ) = m = 1 ρ a m e j ( 2 π f m t + θ m ) + w ( t ) .
Assuming that L undersampling rates are the same as the single-tone case, L discretized versions of x ( t ) can be represented as
x l ( n ) = m = 1 ρ a m e j ( 2 π f m F l n + θ m ) + w ( n ) , l = 1 , . . . , L , n = 0 , . . . , N 1 .
Accordingly, for these tones, ρ simultaneous congruence equation sets can be built up as
f m = n 1 , m F 1 + ε 1 , m F 1 f m = n L , m F L + ε L , m F L , m = 1 , . . . , ρ
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.,
r l , m = ε l , m F l ε l , m = ( k l , m + δ l , m ) / N , l = 1 , . . . , L , m = 1 , . . . , ρ .
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
x l ( n ) = m = 1 ρ 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.

3. Proposed Method

3.1. 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 ( 2 N 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 ( 2 N 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 ( 2 N 1 ) -length one,
x l ( n ) = x l ( n ) , N + 1 n N 1 0 , else .
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 ( 2 N 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.,
x l , q ( n ) = x l ( n ) w 1 ( n + q ) , 0 q N 1
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 l a ( 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
y l , q ( n ) = r = r = + x l , q ( n + r N ) w 2 ( q ) , 0 q N 1 ,
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 l a ( n ) can be derived as
x l a ( n ) = q = 0 N 1 y l , q ( n ) , 0 n N 1 0 , else .
The analysis above studies N sequences { x l , q ( n ) , q = 0 , . . . , N 1 } to determine the unique N-length sequence x l a ( n ) . Usually, the process is referred to as the all-phase preprocessing, accordingly the signal x l a ( n ) is called the all-phase signal. Similarly, the all-phase DFT X l a ( k ) can be calculated by performing the N-point DFT on x l a ( 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 l a ( 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 Figure 1 and Figure 2, respectively.
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 ( 2 N 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 as
δ ^ l , m = 2 v a 1 + v a , if | X l ( k l , m + 1 ) | > | X l ( k l , m 1 ) | 2 v a 1 + v a , else .
Accordingly, the normalized frequency estimate is derived as
ε ^ l , m = ( k l , m + δ ^ l , m ) / N ,
Then, the corrected amplitude is obtained by substituting ε ^ l , m into the exponential term in (12)
a ^ l , m = π δ ^ l , m ( 1 δ ^ l , m 2 ) · | 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).

3.2. 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 , a ^ 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 = a ^ 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 z ^ l , 1 , . . . , z ^ l , ρ , i.e.,
c l , m = arg min m = 1 , . . . , ρ z ^ 1 , m z ^ l , m .
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 estimate f ^ m , m = 1 , . . . , ρ .

3.3. 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).
Algorithm 1: The optimized discrete spectral analysis (ODSA).
  • Input: L coprime integers Γ 1 , , Γ L , the constant N;
  • Output: Frequency estimates f ^ m , m = 1 , . . . , ρ ;
  • Acquire ( 2 N 1 ) -length sequences x l ( n ) , n = N + 1 , . . . , N 1 by undersampling the signal x ( t ) at the sampling rates F l = Γ l N , l = 1 , . . . , L respectively.
  • At any l-th path, preprocess the acquired sequences x l ( n ) according to (14)–(16), from which the N-length all-phase sequence x l a ( n ) , n = 0 , . . . , N 1 can be obtained.
  • Implement N-point DFT on sequences x 1 a ( n ) , . . . , x L a ( n ) to acquire the all-phase DFT results X l a ( k ) , . . . , X L a ( k ) , k = 0 , . . . , N 1 . For any l-th path (l=1,...,L), record the peak DFT bin indices k l , m , m = 1 , . . . , , ρ .
  • For each DFT bin, implement Spectrum Correction through (18)–(21) to acquire ρ L harmonic parameter triples { a ^ l , m , θ ^ l , m , ε ^ l , m } , from which ρ L vector quantities z ^ l , m can be constructed in terms of (22).
  • For any individual z ^ 1 , m of the m-th tone at the 1st path ( m = 1 , . . . , ρ ), find out the peak indicator c l , m amongst { z ^ l , 1 , . . . , z ^ l , ρ } at the l-th path ( l = 2 , . . . , L ) according to the principle in (23).
  • For m = 1 , . . . , ρ , substituting the remainder set { r ˜ 1 , m , . . . r ˜ l , m } , the gcd N, the moduli { F 1 , . . . , F L } into the CRT reconstruction procedure [13] yields the final frequency estimate f ^ m of the m-th tone.

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

4.1. 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 σ n 2 . In this paper, the SNR is defined as
SNR = 10 lg σ s 2 σ n 2
where σ s 2 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 1 a ( n ) , x 2 a ( n ) can be generated.
Following Step 3, two apDFT spectra X 1 a ( k ) , X 2 a ( k ) are acquired and illustrated in Figure 3, from which two sets of peak indices are recorded as { k 1 , 1 , . . . , k 1 , 4 } = { 92 , 178 , 261 , 347 } , { k 2 , 1 , . . . , k 2 , 4 } = { 244 , 278 , 402 , 432 } .
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 { z ^ 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 z ^ l , m are illustrated in Figure 4 and Figure 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.
Following Step 5, in terms of (23), we can acquire two sets of indicators { c 1 , 1 , . . . , c 1 , 4 } = { 1 , 2 , 3 , 4 } , { c 2 , 1 , . . . , c 2 , 4 } = { 3 , 4 , 2 , 1 } . In other words, 8 vector quantities are clustered into ρ = 4 pairs as { z ^ 1 , 1 , z ^ 2 , 3 } , { z ^ 1 , 2 , z ^ 2 , 4 } , { z ^ 1 , 3 , z ^ 2 , 2 } , { z ^ 1 , 4 , z ^ 2 , 1 } , as 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.
Following 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 estimates f ^ 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.

4.2. Performance Analysis on the All-Phase DFT and Traditional DFT

By following the above steps, the multiple frequency estimates can be obtained in the two-path case. Since that the traditional DFT is the most commonly used technique for remainder acquisition, 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) is defined as RMSE = 1 / Q i = 1 Q ( f ^ m i f m ) 2 , 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 Figure 6 and Figure 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.

4.3. 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 Figure 6 and Figure 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.

4.4. 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 / 2 log 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.

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

Author Contributions

Conceptualization, X.H.; methodology, X.H. and L.C.; software, X.H.; validation, X.H. and W.L.; formal analysis, W.L.; investigation, X.H.; resources, X.H.; writing—original draft preparation, C.L.; writing—review and editing, W.L. and X.H.; funding acquisition, W.L. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China under Grant 61671012.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Belega, D.; Petri, D. Frequency estimation by two- or three-point interpolated Fourier algorithms based on cosine windows. Signal Process. 2015, 117, 115–125. [Google Scholar] [CrossRef]
  2. Candan, C. Fine resolution frequency estimation from three DFT samples: Case of windowed data. Signal Process. 2015, 114, 245–250. [Google Scholar] [CrossRef]
  3. Schmidt, R. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, 34, 276–280. [Google Scholar] [CrossRef] [Green Version]
  4. Besson, O.; Stoica, P. Analysis of MUSIC and ESPRIT frequency estimates for sinusoidal signals with lowpass envelopes. IEEE Trans. Signal Process. 1996, 44, 2359–2364. [Google Scholar] [CrossRef]
  5. Bykhovsky, D.; Cohen, A. Electrical Network Frequency (ENF) Maximum-Likelihood Estimation Via a Multitone Harmonic Model. IEEE Trans. Inf. Forensics Secur. 2013, 8, 744–753. [Google Scholar] [CrossRef]
  6. Klapuri, A.P. Multiple fundamental frequency estimation based on harmonicity and spectral smoothness. EEE Trans. Audio Speech Lang. Process. 2003, 11, 804–816. [Google Scholar] [CrossRef] [Green Version]
  7. McKilliam, R.G.; Quinn, B.G.; Clarkson, I.V.L.; Moran, B. Frequency Estimation by Phase Unwrapping. IEEE Trans. Signal Process. 2010, 58, 2953–2963. [Google Scholar] [CrossRef]
  8. Christensen, M.G.; Jensen, S.H. New Results on Perceptual Distortion Minimization and Nonlinear Least-Squares Frequency Estimation. IEEE Trans. Audio Speech Lang. Process. 2011, 19, 2239–2244. [Google Scholar] [CrossRef] [Green Version]
  9. Xia, X.G. On Estimation of Multiple Frequencies in Undersampled Complex Valued Waveforms. IEEE Trans. Signal Process. 1999, 47, 3417–3419. [Google Scholar]
  10. Wang, W.; Li, X.; Wang, W.; Xia, X.G. Maximum Likelihood Estimation based Robust Chinese Remainder Theorem for Real Numbers and its Fast Algorithm. IEEE Trans. Signal Process. 2015, 63, 3317–3331. [Google Scholar] [CrossRef]
  11. Li, X.; Liang, H.; Xia, X.G. A Robust Chinese Remainder Theorem with its Applications in Frequency Estimation from Undersampled Waveforms. IEEE Trans. Signal Process. 2009, 57, 4314–4322. [Google Scholar]
  12. Vaidyanathan, P.P.; Pal, P. Sparse Sensing with Co-prime Samplers and Arrays. IEEE Trans. Signal Process. 2011, 59, 573–586. [Google Scholar] [CrossRef]
  13. Wang, W.; Xia, X.G. A Closed-form Robust Chinese Remainder Theorem and its Performance Analysis. IEEE Trans. Signal Process. 2010, 58, 5655–5666. [Google Scholar] [CrossRef] [Green Version]
  14. Xiao, L.; Xia, X.G. A new robust Chinese remainder theorem with improved performance in frequency estimation from undersampled waveforms. Signal Process. 2015, 117, 242–246. [Google Scholar] [CrossRef]
  15. Huang, X.; Wang, H.; Qin, H.; Nie, W. Frequency Estimator Based on Spectrum Correction and Remainder Sifting for Undersampled Real-Valued Waveforms. IEEE Access 2019, 7, 25980–25988. [Google Scholar] [CrossRef]
  16. Li, G.; Xu, J.; Peng, Y.N.; Xia, X.G. An Efficient Implementation of a Robust Phase-unwrapping Algorithm. IEEE Signal Process Lett. 2007, 14, 393–396. [Google Scholar] [CrossRef]
  17. Huang, X.; Wang, Z.; Ren, L.; Zeng, Y.; Ruan, X. A novel high-accuracy digitalized measuring phase method. In Proceedings of the 2008 9th International Conference on Signal Processing, Beijing, China, 26–29 October 2008; pp. 120–123. [Google Scholar]
  18. Huang, X.; Cui, H.; Wang, Z.; Zeng, Y. Mechanical fault diagnosis based on all-phase FFT parameters estimation. In Proceedings of the IEEE 10th International Conference on Signal Processing Proceedings, Beijing, China, 24–28 October 2010; pp. 176–179. [Google Scholar]
  19. Xiao, H.; Huang, Y.; Ye, Y.; Xiao, G. Robustness in Chinese Remainder Theorem for Multiple Numbers and Remainder Coding. IEEE Trans. Signal Process. 2018, 66, 4347–4361. [Google Scholar] [CrossRef]
  20. Xiao, H.; Xiao, G. Notes on CRT-based robust frequency estimation. Signal Process. 2017, 133, 13–17. [Google Scholar] [CrossRef]
  21. Xiao, H.; Xiao, G. On Solving Ambiguity Resolution With Robust Chinese Remainder Theorem for Multiple Numbers. IEEE Trans. Veh. Technol. 2019, 68, 5179–5184. [Google Scholar] [CrossRef] [Green Version]
  22. Zhang, F.; Geng, Z.; Yuan, W. The algorithm of interpolating windowed FFT for harmonic analysis of electric power system. IEEE Trans. Power Deliv. 2001, 16, 160–164. [Google Scholar] [CrossRef]
Figure 1. (a) Traditional 64-point DFT amplitude spectrum; (b) All-phase 64-point DFT amplitude spectrum.
Figure 1. (a) Traditional 64-point DFT amplitude spectrum; (b) All-phase 64-point DFT amplitude spectrum.
Sensors 20 05066 g001
Figure 2. (a) Traditional 64-point DFT phase spectrum; (b) All-phase 64-point DFT phase spectrum.
Figure 2. (a) Traditional 64-point DFT phase spectrum; (b) All-phase 64-point DFT phase spectrum.
Sensors 20 05066 g002
Figure 3. The upper pane shows the windowed apDFT amplitude spectrum of the undersampled sequence x1a(n). The bottom pane shows that of the undersampled sequence x2a(n).
Figure 3. The upper pane shows the windowed apDFT amplitude spectrum of the undersampled sequence x1a(n). The bottom pane shows that of the undersampled sequence x2a(n).
Sensors 20 05066 g003
Figure 4. Clustering of uncorrected vector quantities.
Figure 4. Clustering of uncorrected vector quantities.
Sensors 20 05066 g004
Figure 5. Clustering of corrected vector quantities z ^ l , m .
Figure 5. Clustering of corrected vector quantities z ^ l , m .
Sensors 20 05066 g005
Figure 6. RMSE results of four tones utilizing the all-phase DFT.
Figure 6. RMSE results of four tones utilizing the all-phase DFT.
Sensors 20 05066 g006
Figure 7. RMSE results of four tones utilizing the traditional DFT.
Figure 7. RMSE results of four tones utilizing the traditional DFT.
Sensors 20 05066 g007
Figure 8. RMSE results of frequency estimation under different number of data acquisition path.
Figure 8. RMSE results of frequency estimation under different number of data acquisition path.
Sensors 20 05066 g008
Table 1. Harmonic parameters of the multi-tone signal.
Table 1. Harmonic parameters of the multi-tone signal.
Tones f m (Hz) A m θ m
s 1 ( t ) 1.3 × 10 9 1 20
s 2 ( t ) 1.4 × 10 9 1 50
s 3 ( t ) 0.9 × 10 9 1 80
s 4 ( t ) 1.12 × 10 9 1 110
Table 2. The vector quantities and normalized frequencies of two undersampled sequences in Step 4.
Table 2. The vector quantities and normalized frequencies of two undersampled sequences in Step 4.
m = 1 m = 2 m = 3 m = 4
l = 1 z ^ 1 , m 0.99 e j 22 . 5 0.98 e j 51 . 2 1.03 e j 81 . 9 1.00 e j 107 . 1
ε ^ 1 , m 0.17980.34750.50900.6780
l = 2 z ^ 2 , m 1.00 e j 107 . 7 1.01 e j 83 . 5 1.01 e j 18 . 2 1.01 e j 49 . 2
ε ^ 2 , m 0.47570.54290.78420.8446
Table 3. The classified CRT remainders and normalized frequencies in Step 5.
Table 3. The classified CRT remainders and normalized frequencies in Step 5.
m = 1 m = 2 m = 3 m = 4
ε ^ 1 , c 1 , m 0.1798 0.3475 0.5090 0.6780
ε ^ 2 , c 2 , m 0.7842 0.8446 0.5429 0.4757
r ˜ 1 , m 0.3038 × 10 6 0.5873 × 10 6 0.8604 × 10 6 1.1459 × 10 6
r ˜ 2 , m 0.8054 × 10 6 0.9193 × 10 6 1.3279 × 10 6 1.4300 × 10 6
Table 4. Final results of frequency recovery.
Table 4. Final results of frequency recovery.
m = 1 m = 2 m = 3 m = 4
f m (Hz) 1.3 × 10 9 1.4 × 10 9 0.9 × 10 9 1.12 × 10 9
f ^ m (Hz) 1.2999 × 10 9 1.3999 × 10 9 0.8999 × 10 9 1.1200 × 10 9

Share and Cite

MDPI and ACS Style

Huang, X.; Cao, L.; Lu, W. Multi-Tone Frequency Estimation Based on the All-Phase Discrete Fourier Transform and Chinese Remainder Theorem. Sensors 2020, 20, 5066. https://doi.org/10.3390/s20185066

AMA Style

Huang X, Cao L, Lu W. Multi-Tone Frequency Estimation Based on the All-Phase Discrete Fourier Transform and Chinese Remainder Theorem. Sensors. 2020; 20(18):5066. https://doi.org/10.3390/s20185066

Chicago/Turabian Style

Huang, Xiangdong, Lu Cao, and Wei Lu. 2020. "Multi-Tone Frequency Estimation Based on the All-Phase Discrete Fourier Transform and Chinese Remainder Theorem" Sensors 20, no. 18: 5066. https://doi.org/10.3390/s20185066

APA Style

Huang, X., Cao, L., & Lu, W. (2020). Multi-Tone Frequency Estimation Based on the All-Phase Discrete Fourier Transform and Chinese Remainder Theorem. Sensors, 20(18), 5066. https://doi.org/10.3390/s20185066

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