An Efficient CRT Based Algorithm for Frequency Determination from Undersampled Real Waveform

The Chinese Remainder Theorem (CRT) based frequency estimation has been widely studied during the past two decades. It enables one to estimate frequencies by sub-Nyquist sampling rates, which reduces the cost of hardware in a sensor network. Several studies have been done on the complex waveform; however, few works studied its applications in the real waveform case. Different from the complex waveform, existing CRT methods cannot be straightforwardly applied to handle a real waveform’s spectrum due to the spurious peaks. To tackle the ambiguity problem, in this paper, we propose the first polynomial-time closed-form Robust CRT (RCRT) for the single-tone real waveform, which can be considered as a special case of RCRT for arbitrary two numbers. The time complexity of the proposed algorithm is O(L), where L is the number of samplers. Furthermore, our algorithm also matches the optimal error-tolerance bound.


Introduction
Chinese Remainder Theorem (CRT) is a fundamental number theory result, which shows the reconstruction of a single integer X from its residues modulo multiple co-prime moduli. It has been extensively used in various applications, such as wireless sensor networks [1,2], coding theory [3][4][5][6][7], phase unwrapping [8,9], and frequency estimation from undersampled waveforms [10][11][12][13][14][15][16][17]. In particular, the CRT-based method enables one to estimate frequencies with exponentially smaller sub-Nyquist rates in a distributed setup. This could significantly reduce hardware cost [18,19]. In practice, errors may occur in the spectrum measurement, while CRT is known highly sensitive to residue perturbation [20]. Moreover, in some applications of multiple parameter estimation, we may need to recover multiple real numbers simultaneously. To this end, many works have been proposed during the last two decades to solve the two issues, which can be summarized as follows.
(i) Robustness: On the one hand, to make CRT robust against small errors in residues, Wang et al. introduced a common factor Γ as redundancy to the co-prime moduli {M 1 , M 2 , . . . , M L } in a form {m l = M l Γ|l = 1, 2, . . . , L}. This forms the foundation of the first closed-form Robust CRT (RCRT) for a single real number [20]. RCRT can recover the folding number X/Γ once the error in each residue is upper bounded by Γ/4. Hence, one can ensure the reconstruction error is upper bounded by Γ/4. The Γ/4 error tolerance bound is also proved to be tight in the follow-up works [21]. (ii) Residue Ambiguity: On the other hand, since the observed residues are unordered, there is no clear correspondence between N numbers {X i |i = 1, 2, . . . N} and residues in each residue set R l = {r i,l |i = 1, 2, . . . , N}, l = 1, 2, . . . , L. Here, r i,l denotes the residue of X i modulo m l . Thus, the residue ambiguity makes reconstruction much more complicated for multiple numbers. When N = 2, Xiao et al. proposed a robust generalized CRT, addressing the residue ambiguity by carefully-designed quadratic symmetric polynomials [22]. It is shown that the correspondence between these two numbers and residues can be uniquely determined while the error bound is sacrificed to Γ/8 [23]. As shown in recent works [24,25]; theoretically, one can approach the optimal error bound Γ/4 independent of N when the least common multiple of moduli is sufficiently large. However, to the best of our knowledge, no existing polynomial-time algorithm matches the optimal bound.
In this paper, we focus on CRT based algorithm for frequency determination from the undersampled real waveform. The proposed method can be applied in a sensor network with low power and low transmission rates sensors [26,27] or Synthetic Aperture Radar (SAR) imaging of moving targets [28]. However, in the real waveform scenario, the CRTbased method encounters both the above-mentioned challenges, robustness, and residue ambiguity, simultaneously.
Notably, the real waveform sampling needs less hardware, i.e., only one Analog-to-Digital Converter (ADC) per sampling frequency is required in real waveform sampling rather than two ADCs in complex waveform [29]. However, existing CRT methods for the complex waveform cannot be applied to the real waveform directly due to the existence of the spurious peak [24]. In this paper, we set out to solve these mentioned issues. Our main contributions can be concluded as follows.

•
We present the first polynomial-time closed-form RCRT for frequency determination from undersampled single-tone real waveform, which provides a feasible and efficient solution. Moreover, the proposed method fixes the gap in the CRT-based method for frequency determination for the real waveform case. • By fully utilizing the prior knowledge of the real waveform, we reach the optimal error tolerance bound, i.e., Γ/4, which is twice better than the best-known robust generalized CRT proposed in [23].
The remaining content is organized as follows. In Section 2, we give an overview of the problem formulation. Section 3 details our closed-form reconstruction for the real waveform. In Section 4, we present some simulation results to support the theory. In Section 5, we discuss and interpret the simulation results. The conclusion is drawn in Section 6.

Problem Formulation
We first describe the frequency estimation model from the undersampled real waveforms.

Signal Model and Sampling
A sinusoidal waveform is defined as where A denotes the amplitude, X represents the frequency. Sampling x(t) with L ADCs at frequency rates of {m l |l = 1, 2, . . . , L} [24,30], where max l m l < 2X, i.e., the sampling rates are below the Nyquist rate, we have Applying the m l -point Discrete Fourier Transform (DFT) to x m l [u] [31,32], we obtain Here, δ( * ) is the the Kronecker delta function, i.e., δ(k − X m l ) equals 1 when k = X m l or 0 otherwise, where k represents a frequency bin and X m l denote the residue of X modulo m l . Clearly, the locations of the spectrum peaks correspond to the residues X m l and −X m l , which leads to two symmetric peaks over the frequency spectrum domain in the noiseless case. Thus, one can recover the frequency X with sampling rates (moduli) m l and the locations of the spectrum peaks (residues) X m l via CRT.

Noise Model and RCRT Procedure
In the following, we further consider the noisy case and review RCRT. Still, let X i ∈ {X 1 , X 2 } represent the real number to be recovered, where X 1 > 0 and X 2 = −X 1 . The moduli are in a form {m l = M l Γ|l = 1, 2, . . . L}, where {M l } are pairwise co-prime. r i,l = X i + ∆ i,l m l denotes the erroneous residue of X i modulo m l , where ∆ i,l represents the underlying error such that |∆ i,l | < Γ/4. Moreover, r c i = X i Γ denotes the common residue of X i .r c i,l = X i + ∆ i,l Γ = r c i + ∆ i,l Γ denotes the erroneous common residue. In practical cases,r c i,l is calculated fromr i,l based on the number theory, i.e., r c i,l = X i + ∆ i,l Γ = X i + ∆ i,l M l Γ Γ = r i,l Γ , which ensures thatr c i,l andr i,l share the same ∆ i,l . For clarity, all the notations are listed in Table 1. In the following, we aim to estimate the real number X 1 with known erroneous residuesr i,l and moduli m l .

Notations Explanation
Since X i = X i /Γ Γ + r c i , we recover X i by estimating the folding number X i /Γ and common residue r c i successively. We adopt the reconstruction framework proposed in [24], which consists of three steps: (i) Estimate the folding number: Based on the fact that X i = k i m l + r i,l = k i M l Γ + r i,l , where k i ∈ Z, we have X i /Γ = k i M l + r i,l /Γ = k i M l + (r i,l − r i,l Γ )/Γ. Clearly, r i,l Γ = X i M l Γ Γ = X i Γ = r c i . Thus, we have X i /Γ = k i M l + (r i,l − r c i )/Γ. By taking the modulo arithmetic, one can obtain From (4), the folding number X i /Γ is estimated by the equation below via CRT [24], where q i denotes the estimation of X i /Γ . (ii) Estimate the common residues: Calculate ∑ L l=1r c i,l /L as the estimation of the common residue r c i . (iii) Estimate the number: Based on X i = X i /Γ Γ + r c i , X i is reconstructed bŷ whereX i represents the estimation of X i .

Key Issues in Real Waveform
(i) Robustness: Trivially estimating the folding number by (5) may lead to large errors due to the ambiguity of (r i,l −r c i,l )/Γ. In other words, since r c i ∈ [0, Γ) and |∆ i,l | < Γ/4, (r i,l −r c i,l )/Γ must satisfy one of the three subcases below based on (4) [33], If r c i + ∆ i,l 1 and r c i + ∆ i,l 2 fall into different subcases in (7), where l 1 , l 2 ∈ {1, 2, . . . , L}, simply aggregating them via CRT will bring unpredictable reconstruction errors. Thus, we need to unify (r i,l −r c i,l )/Γ such that all of them fall into one subcase in (7) to ensure robustness. This can be achieved by sortingr c i,l , wherer c i,l = r c i + ∆ i,l Γ , in the order such that the corresponding ∆ i,l are in an ascending order for each i. However, the above operation is only implementable when |∆ i,l | < Γ/8 [34], while it still remains open in the generic setup |∆ i,l | < Γ/4. (ii) Residue Ambiguity: Due to the loss of the correspondence between X i andr i,l , we cannot clusterr i,l corresponding to X i to calculate q i from (5) for each i.

Robust Reconstruction for Frequency Estimation of Single-Tone Real Waveform
This section presents the polynomial-time RCRT-based frequency estimation for a noisy single-tone real waveform. Before proceeding, the following notations are introduced.
We first define a metric to represent the minimum circular distance betweenr c i,l 1 and r c i,l 2 on the circle of length Γ, i.e., ,r c i,l 2 ) represents the interval whose length is maximal. As shown in Figure 1, when i = 1, the maximum interval is I(r c 1,1 ,r c 1,3 ); when i = 2, I(r c 2,1 ,r c 2,3 ) is the maximum one.

The Order of Residues
Now, we consider the first key issue stated in Section 2.3, i.e., sortingr c i,l in the order such that the corresponding errors ∆ i,l are in ascending order for each i, wherẽ r c i,l = r c i + ∆ i,l Γ . According to [21], sorting is equivalent to finding a cutting point ξ on the circle of length Γ and stretching it to a real axis. If ξ ∈ max l I( r c i,l 1 , r c i,l 2 ) for each i, the shifted common residuesr c i,l on the real axis are sorted in ascending order of ∆ i,l . For example, in Figure 1, Γ/2 is not in the maximum intervals, i.e., I(r c 1,1 ,r c 1,3 ) and I(r c 2,1 ,r c 2,3 ). Then, cutting the circle at Γ/2 leads tor c 2,3 <r c 2,2 <r c 2,1 andr c 1,1 <r c 1,2 <r c 1,3 sorted in ascending order of ∆ i,l , i.e., ∆ 2,3 < ∆ 2,2 < 0 < ∆ 2,1 and ∆ 1, shown in Figure 2a. On the contrary, if we cut the circle at 0, where 0 ∈ I(r c 1,1 ,r c 1,3 ), ∆ 1,1 breaks the ascending order, as shown in Figure 2b. Here,r c 1,2 <r c 1,3 <r c 1,1 , but ∆ 1,2 , ∆ 1,3 , ∆ 1,1 are in non-ascending order. However, the key remaining problem is how to find the proper cutting point without the correspondence betweenr c i,l and X i , i.e., max l I(r c i,l 1 ,r c i,l 2 ) is unknown. To this end, ref. [21] sacrifices the error bound to Γ/8. Nonetheless, we reach the error bound Γ/4 by using the symmetry of residues, i.e., That is to say, r c 1 and r c 2 are axially symmetric about line α shown in Figure 1.
,r c i,l 2 ) cannot contain 0 and Γ/2 simultaneously. Based on the symmetry, max l I(r c 1,l 1 ,r c 1,l 2 ) ∪ max l I(r c 2,l 1 ,r c 2,l 2 ) cannot contain both 0 and Γ/2. Thus, either 0 or Γ/2 is the cutting point. To figure out the cutting point, we state Lemma 1, which is proved in Appendix A, that once 0 ∈ max l I(r c Thus, the unknown max l I(r c i,l 1 ,r c i,l 2 ) problem is converted to distance comparison, i.e., min il d(0,r c i,l ) and min il d(Γ/2,r c i,l ). Before giving Lemma 1, we define Operation 1 and 2 corresponding to the cutting point is 0 and Γ/2, respectively.

Residue Ambiguity
Withr c i,l sorted in ascending order of ∆ i,l , (r i,l −r c i,l )/Γ fall into one subcase in (7) [24]. Now, we discuss the second key issue, i.e., residue ambiguity. If we can divider i,l into two sets corresponding to X 1 and X 2 , respectively, the folding number X i /Γ can be estimated based on (5) [34]: However, the correspondence betweenr i,l and X i is unknown. To determine the correspondence, Li et al. proposed a scheme for positive numbers, which cannot be directly applied to the real waveform since X 2 < 0 [23]. To solve this issue, we form a quadratic equation by the prior condition X 2 = −X 1 . First, we consider the two residues {(r 1,l −r c 1,l )/Γ, (r 2,l −r c 2,l )/Γ} as a pair for each l. By multiplying each pair, we can reconstruct q 1 q 2 via CRT based on (11) [22], i.e., Then, with X 2 = −X 1 , it can be proved that either q 2 = −q 1 or q 2 = −q 1 − 1 holds, which is stated in Lemma 2 and proved in Appendix B. Hence, we can form a quadratic equation in one unknown by replacing q 2 with −q 1 or −q 1 − 1 in (12) based on Lemma 2. In a nutshell, the residue ambiguity is addressed by solving one of the two quadratic equations below via CRT, corresponding to q 2 = −q 1 and q 2 = −q 1 − 1, respectively.
Lemma 2. If |∆ i,l | < Γ/4, {q 1 , q 2 } must fall into one of the following two cases: • q 2 = −q 1 − 1, when Operation 1 is the appropriate operation and performed onr c i,l , where when Operation 2 is the appropriate operation and performed onr c i,l . If r c 1 ∈ [0, Γ/2),

Simulation Results
In this section, we first present some simulations to verify our proposed theory. Then the simulation results are shown to demonstrate the performance of the proposed method compared with that of the robust generalized CRT [23] and searching−based algorithm [29].
In the following, we first consider the estimation error versus the error upper bound for our proposed theory, i.e., Theorem 1. To begin with, the simulation setup is given as follows. The moduli are m l = {11 × 80, 13 × 80, 17 × 80}, where the greatest common factor Γ = 80 and the maximal error level τ ∈ {1, 2, . . . , 25}. Based on Theorem 1, τ needs to be bounded by Γ/4 = 20 to ensure robustness.
For a trial, one unknown real number X is chosen randomly, which belongs to the dynamic range [0, 3880), where the negative duplicate is −X. Moreover, 10,000 trials are implemented for each τ. Figure 3 shows the mean absolute error E τ between the estimateX and the true number X for each error bound. The mean absolute error E τ is defined as below, where E trials denotes the mean of all the trials,X and X are the estimate and true number in a trial, respectively. Clearly, E τ is less than τ when τ ≤ Γ/4 = 20, which matches well with our conclusion. Once τ exceeds the error bound, the reconstruction error increases rapidly.
In Figure 4, we present the curve of the probability of failure P e versus the error bound τ, where P e = P(|X − X| > τ).
One can see that when τ ≤ Γ/4, the probability of failure is zero while non−zero when τ exceeds the bound. In a word, if τ < Γ 4 , the reconstruction error is linearly bounded by τ, the probability of which is 1.  Next, we compare the performance of the proposed algorithm with that of the robust generalized CRT for two numbers in [23] and the searching−based algorithm in [29]. We consider the real sinusoidal waveform case and select L = 3 sampling rates (moduli) in a form m[1 : L] = Γ × {11, 13, 17}, which share a greatest common factor Γ. We test different sampling rates where Γ = {40, 80}. The unknown frequency X is randomly selected from the range [0, 48.5 × Γ). Each noise ∆ i,l is assumed to be some independent uniform noise within (−τ, τ), where τ varies from 1 to 25.
We repeat 5000 trials for each selection of Γ and τ. On the on hand, the root mean square error (RMSE) is investigated, where Figure 5a,c show that our method outperforms the best known robust generalized CRT, where the maximal error tolerance is improved from Γ/8 to Γ/4. Morever, our method performs as well as the best searching−based method when the maximal error level is less than Γ/4. On the other hand, we compare the test fail rate (TFR). We say that the test fails when As shown in Figure 5b,d, if τ ≤ Γ/4, the estimation error is bounded by Γ/4, the probability of which is one. Once the maximal error level exceeds Γ/4, the reconstruction error is almost unpredictable. In a word, our method outperforms the robust generalized CRT while slightly worse than the searching−based algorithm when τ > Γ/4. However, it's worth pointing out that our method provides a closed−form solution that cannot be realized by the searching−based method. Then, we consider the real running time consumption, where the computing equipment is Lenovo xiaoxin Pro 13. The real running time of our method that runs for 125,000 times is about 7.96 s, while the robust generalized CRT proposed in [23] requires about 86.05 s since the algorithm involves a lot of loops. The searching-based method proposed in [29] sightly outperforms our method, which only needs about 5.75 s.

Discussion
The experiment results in Figure 5 suggest a clear improvement in the error bound from Γ/8 to Γ/4 compared with the method proposed in [23]. The reason why we can improve the error bound is that we fully utilized the prior knowledge of the real waveform, i.e., symmetry. For the real waveform, the real peak and the corresponding spurious peak are symmetric at about 0 points in the spectrum. Thus, the frequency determination from undersampled single-tone real waveform can be formulated as RCRT for two numbers {X 1 , X 2 }, where these two numbers are in a form X 1 = −X 2 . Based on this symmetry, the corresponding error-free common residues {r c 1 , r c 2 } are symmetric on the circle of length Γ. The geometric property of symmetry ensures that even if the error bound is improved to Γ/4, we can still shift the erroneous common residues correctly to obtain a robust reconstruction. In addition, our algorithm is also highly efficient according to the real running time and the theoretical analysis. We use the prior condition of the real waveform to form a quadratic equation in one unknown to determine the folding numbers, which realizes the high efficiency of the algorithm.
In summary, our proposed method provides a feasible solution for the frequency determination from the undersampled single-tone real waveform. In addition, we complete the study of CRT-based frequency determination from undersampled waveform, which shows that the optimal error tolerance bound can be achieved in the real waveform case. The limitation of our proposed method is that since it is based on the prior knowledge of the real waveform and the prior condition is invalid; it cannot handle the complex waveform. In addition, this algorithm cannot deal with the case of multiple frequency estimation from undersampled real waveforms. We will investigate these problems in our future studies.

Conclusions
We proposed the first polynomial-time RCRT-based frequency estimation for a noisy single-tone real waveform, which matches the optimal error bound. The proposed method can be applied in SAR imaging of moving targets or sensor networks where the sampling rate may be lower than the Nyquist rate of the input signal. The time complexity of the proposed method is linear to the number of samplers. Moreover, the proposed method can estimate the frequency from the real waveform by sub-Nyquist rates, which reduces the cost and system size, especially in sensor networks that require noticeable sensors. We believe the method can be further extended to the multiple frequencies case.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Proof of Lemma 1
Proof. For brevity, we only consider the case when r c 1 ≤ r c 2 , i.e., r c 1 is on the right half of the circle Γ. For the case r c 1 ≥ r c 2 , which is obviously symmetric with r c 1 ≤ r c If 0 ∈ max l I(r c 1,l 1 ,r c 1,l 2 ), there must exist at least oner c 1,k to the left of the point 0, and anotherr c 1,l to the right of the point 0, where k, l ∈ {1, 2, . . . , L}. Since r c 1 is on the right half of the circle Γ whiler c 1,k = r c 1 + ∆ 1,k Γ is to the left of the point 0 and |∆ 1,k | < Γ/4, we have r c 1 ∈ [0, Γ/4). Based on the fact that r c 2 = Γ − r c 1 , r c 2 ∈ (3Γ/4, Γ]. The same conclusion can be derived if 0 ∈ max l I(r c 2,l 1 ,r c 2,l 2 ). Since r c 1 ∈ [0, Γ/4) andr c 1,k is to the left of the point 0, we obtain min Now, we calculate min il d( Γ 2 ,r c i,l ). When i = 1, since r c 1 ∈ [0, Γ/4) and |∆ i, In conclusion, when 0 ∈ max l I(r c Since Cases (1) and (2) cannot hold simultaneously, we have Γ 2 ∈ max l I(r c i,l 1 ,r c i,l 2 ) for each i. Thus, cutting the circle at Γ/2, i.e., applying Operation 2 onr c i,l leads tor c i,l sorted in ascending order of ∆ i,l .
(1) min d(0,r c i,l ) ≥ min d( Γ 2 ,r c i,l ) (2) min d(0,r c i,l ) < min d( Γ 2 ,r c i,l ) Case (1): Since Operation 1 is the appropriate operation based on Lemma 1, we havê r c i,l =r c i,l based on (9). Therefore, replacingr c i,l withr c i,l , (11) is equal to Then, we discuss r c i + ∆ i,l , which determines q i in (A5). With r c i ∈ [0, Γ) and |∆ i,l | < Γ 4 , r c i + ∆ i,l must fall into one of the three subcases: ) are residues of ( X 1 Γ , X 2 Γ ) modulo M l , where Therefore, when Operation 1 is the appropriate operation, q 2 = −q 1 − 1 and q 1 = X 1 Γ . Case (2): In this case, Operation 2 is the appropriate operation based on Lemma 1. Then, based on (10),r c i,l =r c i,l whenr c i,l ∈ [0, Γ/2); otherwise,r c i,l =r c i,l − Γ. Thus, replacinĝ r c i,l withr c i,l , (11) is equal to Likewise, we discuss r c i + ∆ i,l in the following to figure out q i , where r c i + ∆ i,l must fall into one of the five subcases: are residues of X i Γ modulo M l based on (7). So q i = X i Γ . Subcase (c): Likewise, since r c i + ∆ i,l ∈ [ 3Γ 4 , Γ), we haver c i,l = r c i + ∆ i,l Γ ∈ [ 3Γ 4 , Γ). Therefore, based on (A6), we have q i ≡˜r i,l −r c i,l Γ + 1. From (7), one can obtainr i,l −r c i,l Γ are residues of X i Γ modulo M l , which leads to q i = X i Γ + 1.
. Therefore, based on (A6), we have q i ≡˜r i,l −r c i,l Γ + 1. From (7), one can obtainr i,l −r c i,l Γ are residues of . Therefore, based on (A6), we have q i ≡˜r i,l −r c i,l Γ . From (7), one can obtainr i,l −r c i,l Γ are residues of X i Γ + 1 modulo M l , which results in q i = X i Γ + 1. Since |∆ i,l | < Γ/4, only subcases (b) and (d) or subcases (c) and (e) can occur simultaneously. Clearly, subcases (b) and (d) result in the same q i , i.e., q i = X i Γ . Similarly, subcases c) and e) lead to the same q i , where q i = X i Γ + 1. Based on the discussions above, we start to consider the residue pair {r c 1 , r c 2 }. Since r c 2 = Γ − r c 1 , {r c 1 , r c 2 } must fall into one the four subcases: are two scenarios for the reconstruction error: (i) Operation 1 is the appropriate operation and applied; (ii) Operation 2 is the appropriate operation and implemented. When Operation 1 is applied: In this case,r c 1,l =r c 1,l based on (9). Then, we discuss r c 1,l = r c 1 + ∆ 1,l Γ , where r c 1 + ∆ 1,l must fall into one of the three subcases: (a) r c , a contradiction to the fact that Operation 1 is the appropriate operation. Likewise, when r c 1 + ∆ 1,l falls into subcase (c), we have min il d(0,r c i,l ) < Γ 4 < min il d( Γ 2 ,r c i,l ), which contradicts the fact that Operation 1 is the appropriate operation. Therefore, only subcase a) can happen, i.e., r c 1 + ∆ 1,l ∈ (0, Γ). Thus,r c 1,l = r c 1 + ∆ 1,l Γ = r c 1 + ∆ 1,l .