# Computationally Efficient Robust Chinese Remainder Theorem for Frequency Estimation from Undersampled Waveforms

Authors: Xiaoping Li; Wenjie Wang; Bobin Yao; Bin Yang; Qinye Yin
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We present an efficient algorithm based on the robust Chinese remainder theorem (CRT) to perform single frequency determination from multiple undersampled waveforms. The optimal estimate of common remainder in robust CRT, which plays an important role in the final frequency estimation, is first discussed. To avoid the exhausted searching in the optimal estimation, we then provide an improved algorithm with the same performance but less computation. Besides, the sufficient and necessary condition of the robust estimation was proposed. Numerical examples are also provided to verify the effectiveness of the proposed algorithm and related conclusions.


---

## Body

## 1. Introduction

Frequency estimation from noisy sinusoid is a hot research topic in engineering, which is applied in almost all engineering fields, such as communications, radar, astronomy, medicine, measurements, and geophysical exploration [1–6]. At present, many of the strategies are based on the Nyquist sampling theorem, for instance, discrete Fourier transform (DFT) algorithm [7–9], the least-error squares algorithm [10–13], and equally spaced samples algorithm [14]. However, one of the main drawbacks of these methods is that the sampling frequency is usually proportional to the original frequency, which requires much higher speed analog-to-digital convertor. Recently, frequency estimation from undersampled waveforms has been studied in [15–17]. The basic idea of these methods is to fuse multiple remnant information provided by multiple sub-Nyquist sampled sequences. Then the original frequency can be reconstructed by using the Chinese remainder theorem (CRT) [18–20]. The biggest merit of these methods is that they greatly decrease the rigorous requirement for hardware; that is, a very high frequency can be determined by several lower rate sample devices. In addition, the method is also used in synthetic aperture radar [21, 22] and wireless sensor network [23].The problem of frequency determination in undersampled condition is summarized as follows.Assume that a single frequency signal is modelled as a complex sinusoid form:(1)s0(t)=Aexp(j2πf0t),
where f0=NHz is an unknown frequency we are interested in the signal and it may be high, A is the signal amplitude, and j=-1. In practical applications, noises are unavoidable and thus lead to the received signal having the following expression:
(2)s(t)=Aexp(j2πf0t)+w(t),
where w(t) denotes noise.Let the sampling rates beM1Hz, M2Hz,…,MLHz. Then we can get the corresponding sampled signal as follows:
(3)si(nMi)=Aexp(j2πNnMi)+w(nMi),i=1,2,…L.
If the sampling duration is T seconds, then
(4)si(nMi)=Aexp(j2πNnMi)+w(nMi),0≤n<TMi.TheMi-point DFT of each of the aforementioned sampled signals may provide an estimate of the original frequency. BecauseMiare much smaller than N in practical applications, the Mi-point DFT only provides a folded remainder ri of N modulo Mi. These remainders may be erroneous when the signal is contaminated by the additive noise. So, the above frequency estimation problem is equivalent to determining the original number N from the erroneous remainders r^1, r^2,…,r^L.This problem has been studied in [17] where a fast robust CRT algorithm is proposed. However, the number of searches is high and no closed-form solution is provided. It has been validated that this method has a poorer performance than the improved closed-form robust CRT (ICRT), which is proposed in [24]. The latter method improved the performance because of the remainders’ proper choosing. However, the estimation of common remainder requires searching, which leads to a high computation complexity. Motivated by this, we propose a computational efficient algorithm which has the same performance as the ICRT method but with much less computation.The remaining of this paper is organized as follows. In Section2, we present the CRT and analyse its error sensitivity. In Section 3, we propose the robust CRT frequency estimation algorithm. Besides, we give a sufficient and necessary condition of the robust estimation based on the proposed algorithm. In Section 4, we apply the proposed algorithm for signal in additive noise with low sample rates. We end with concluding remarks in Section 5.


## 2. The Error Sensitivity of CRT

In this section, we first briefly describe the extended CRT where all the moduli have a common factor larger than 1. We then briefly describe the error sensitivity of CRT.For convenience, we denote by〈x〉M the remainder of x modulo M.Let0<M1<M2<⋯<ML be the moduli, and let r1, r2,…,rL be the remainders of the positive integer N; that is,
(5)N≡rimodMi,i=1,2,…,L,
where 0≤ri<Mi. It is not hard to see that the original number N can be uniquely reconstructed from its remainders if and only if N is less than the least common multiple of the remainders. In the case when all the moduli are coprime, the CRT has a simple formula. In this paper, we consider the case when any paired moduli Mi have the same greatest common remainder (gcd) M>1 and the remaining integers factorized by M are coprime integers Γi; that is,
(6)ri=〈N〉MΓi,i=1,2,…,L.
It follows from (6) that
(7)〈ri〉M=〈N〉M,i=1,2,…,L.
One can see from (7) that all of the remainder ri moduli M have the same value, which is called the common remainder, denoted by rc.Letqi=(ri-rc)/M, and let N0=(N-rc)/M. Then we obtain from (5) that
(8)N0≡qimodΓi.
According to the CRT, if and only if 0≤N0<Γ=Γ1Γ2⋯ΓL, N0 can be uniquely reconstructed as
(9)N0=〈∑i=1L‍γ¯iγiqi〉Γ,
where γi=Γ/Γi, γ¯i is the modular multiplicative inverse of γi modulo Γi. Therefore, N can be uniquely reconstructed by
(10)N=MN0+rc.Now, we consider the case when the remainders have errors.Let theith erroneous remainder be
(11)r^i=ri+Δri,
where Δri denotes the error. In order to reconstruct the original number N from its erroneous remainders r^1, r^2,…,r^L, we need to determine rc and N0 according to (10). Suppose that the estimate of the common remainder is r^c. Then qi can be estimated by
(12)q^i=[r^i-r^cM]=qi+[Δri-ΔrcM],
where Δrc=r^c-rc and [·] stands for the rounding integer.According to (12), if the error of the common remainder Δrc is constrained by
(13)-M2+Δri≤Δrc<M2+Δri,
then q^i=qi. Consequently, N0 can be accurately reconstructed. It follows from (10) that the estimate of the original number is
(14)N^=MN0+r^c.
Obviously, the estimate error is Δrc.Now, we consider condition (13), which leads to the correct estimate of N0. It can be rewritten as
(15)-M2+Δri+rc≤r^c<M2+Δri+rc.Note that the fact thatr^c does not satisfy (15) will lead to the error estimate of qi by (12). As γ¯iγi is usually a large integer, this may cause a large error of N0. For example, let N0=51, Γ1=11, Γ2=13, and Γ3=17. Then we have q1=7, q2=12, and q3=0. Suppose that the contaminated data are q^1=6, q^2=12, and q^3=0. Then N0 can be reconstructed as 2661. However, in the case when all errors are the same, the error of reconstruction is reduced significantly. Suppose that the estimate of qi is q^i=qi+Δq for i=1,2,…,L. Then we have N^0=N0+Δq, which means that the estimate error of N0 is Δq. Therefore, the estimate of the common remainder r^c is significant to the estimation.


## 3. Robust CRT Frequency Estimation Algorithm

In this section, we give the optimal estimate of the common remainder firstly; then we propose the efficient frequency estimation algorithm. Before introducing our method, we introduce a kind of circular distance function.Definition 1.
For real numbersx and y, the circular distance of x to y for cycle M is defined as
(16)d(x,y∣M)=Δx-y-k0M,
where k0=[(x-y)/M] and M∈ℝ∖{0}.Putting(17)r^ic=〈r^i〉M,
then we can conclude that r^ic may be different from each other due to errors Δri. The optimal estimate r^c is the point which minimizes the summation of the deviation between r^ic and variable x, where x∈[0,M). Based on the circular distance defined above, we can obtain the summation of deviation as follows:
(18)D(x)=∑i=1L‍d2(r^ic,x∣M).
Then, the optimal estimate r^c can be determined by
(19)r^c=argmin0≤x<MD(x).
Next, we give the optimal estimate in the following theorem.Theorem 2.
The optimal estimater^c in (19) can be determined by
(20)r^c=argminx∈Ω∑i=1L‍d2(r^ic,x∣M),
where Ω={〈(1/L)∑i=1Lr^i-(M/L)k〉M,k=0,1,…,L-1}.Proof.
According to the definition of the circular distance, we have(21)d(r^ic,x∣M)=d(r^i,x∣M).
Letfi(x)=d2(r^i,x∣M). One can easily see that fi(x) is continuous in field of real number. Now, we consider the differentiability of fi(x) at xt which satisfies
(22)d(r^i,xt∣M)=-M2.
Forx∈(xt-ϵ,xt), where ϵ is a small positive number, we have d(r^i,x∣M)=xt-M/2-x. It follows that the left-sided derivative of fi(x) at xt is M.
Forx∈(xt,xt+ϵ), we obtain d(r^i,xM)=xt+M/2-x. It follows that the right-sided derivative of fi(x) at xt is -M. Thus, fi(x) is not differentiable at xt.
Letf(x)=∑i=1Ld2(r^i,x∣M). Then we can obtain
(23)f+′(xt)-f-′(xt)=-2kM<0,
where k is the number of remainders satisfying (22). Thus, f(x) does not have a local minimum value at xt. Note that f(x) is differentiable at all x∈ℝ∖{xt}, so the minimum point must satisfy f′(x)=0, that is, r^c in (19), such that
(24)(∑i=1L‍d2(r^i,x∣M))′=0.
After equating the result, we obtain
(25)x=1L∑i=1L‍r^i-ML∑i=1L‍ki.
Since ∑i=1Lki∈ℤ and r^c∈[0,M), it follows that r^c belongs to
(26)Ω={〈1L∑i=1L‍r^i-MLk〉M,k=0,1,…,L-1}.
Hence, the optimal estimate r^c in (19) can be simplified as
(27)r^c=argminx∈ΩD(x).
Combining (26) and (27), we can draw the conclusion.Theorem2 tells us that the estimate of common remainder belongs to the set Ω. It is not difficult to find that the cardinality of the set equals the number of the remainders. So, the number of searches is reduced sharply when compared with other searching methods. After optimal estimate r^c is determined, the original frequency can be reconstructed. For convenience, we conclude the proposed efficient algorithm in Algorithm 1.Algorithm 1:Summary of the proposed algorithm.
Step  1.  From contaminated remainders r^i, calculate L0:
L0=argmink=0,1,…,L-1∑i=1Ld2(r^i,1L∑i=1Lr^i-MLk∣M)   (*)
Step  2.   Calculate r^c:
r^c=〈1L∑i=1Lr^i-MLL0〉M.         (**)
Step  3.  Calculate q^i  (i=1,2,…,L) by (12).
Step  4.  Calculate N^0:
N^0=〈∑i=1Lγ-iγiq^i〉M.        (***)
Step  5.  Calculate N^:
N^=MN^0+r^c.       (****)Based on the proposed algorithm, we have a sufficient and necessary condition of the robust estimate as shown in Algorithm1.Theorem 3.
LetΔr¯=(1/L)∑i=1L‍Δri. Assume that |Δr¯|<M/2 and |Δri-Δr¯|<M/2. N is a real number in the range [M,M(Γ-1)). Then
(28)N^-N=Δr¯
holds if and only if
(29)L0=〈∑i=1L‍qi〉M.Proof.
We proof the sufficiency.
From (**), we have(30)r^c=〈1L∑i=1L‍r^i-MLŁ0〉M=〈rc+Δr¯〉M.
Since |Δr¯|<M/2 and 0≤rc<M, it follows that
(31)-M2≤rc+Δr¯<3M2.
We make this into three cases as follows.
Case A. (-M/2≤rc+Δr¯<0).
According to (30), we obtain r^c=rc+Δr¯+M. It follows from (12) that q^i=qi-1 for i=1,2,…,L. From (9), we have
(32)N^0=〈∑i=1L‍γ¯iγi(qi-1)〉M=N0-1.
Here, we use the conclusion that 〈∑i=1L‍γ¯iγi〉M=1. Consequently,
(33)N^=M(N0-1)+r^c=N+Δr¯.
Case B. (0≤rc+Δr¯<M).
According to (30), we obtain r^c=rc+Δr¯. It follows from (12) that q^i=qi. By using CRT, we have N^0=N0. Consequently, N^=MN0+r^c=N+Δr¯.
Case C. (M≤rc+Δr¯<3M/2).
According to (30), we obtain r^c=rc+Δr¯-M. It follows from (12) that q^i=qi+1. From (9), we have N^0=N0+1. Consequently, N^=N+Δr¯. This completes the proof of the sufficiency.
We then proof the necessity.
Combination of (10) and (****) yields
(34)N^-N=M(N^0-N0)+(r^c-rc)=Δr¯.
This leads to
(35)|M(N^0-N0)+(r^c-rc)|<M2.
Since -M<r^c-rc<M, it follows that
(36)|N^0-N0|<32.
Note that both of N^0 and N0 are integers; we have three possible cases: N^0-N0=-1, N^0=N0, and N^0-N0=1.
Case A. (N^0-N0=-1).
By (34), we have r^c-rc-M=Δr¯. That is, r^c=rc+Δr¯+M.
Case B. (N^0=N0).
By (34), we have r^c=rc+Δr¯.
Case C. (N^0-N0=1).
By (34), we have r^c-rc+M=Δr¯. That is, r^c=rc+Δr¯-M.
Therefore, we have(37)r^c=〈rc+Δr¯〉M.
Combing (**) and (37), we can draw the conclusion.


## 4. Simulation Results

In the simulations, the coprime numbers fromΓ1 to Γ7 are 5, 7, 11, 13, 19, 23, and 29, respectively. Two possibilities of M are considered: M=100 and M=200. The frequencies which are to be estimated is taken integers randomly and uniformly distributed in the range (MΓ/Γ1,MΓ) for each M. The additive noise w(t) in (4) is additive white Gauss noise. The sampling frequencies areMΓi, which are much smaller than the unknown frequencies. The observation of the time duration is 1 s, and the number of the trials is 10000 for each signal-to-noise ratio (SNR).We compare the proposed method with the ICRT algorithm which is the most effective one of the existing methods. To assess the performance, we take the trail-fail rate (TFR) as the measurement. In each trial, if the estimate of the frequency is(38)|f^0-f0|<M2,
then the trial passed; otherwise, the trail failed.Figure1 gives the curve of the TFR on different SNR. It shows that the probabilities of robust estimation are almost equal when the two methods have the same M. That is to say, the two methods have nearly the same performance. The figure also shows that the performance of the estimation is better when M=200 than that when M=100, which is in accordance with the theory. From (13) we can conclude that the larger M can lead to the larger error tolerance.Figure 1
TFR versus SNR.Figure2 presents the curve of the total running time over the number of remainders. The two methods are carried out under the same operational conditions. One can see that our proposed algorithm has nearly equal time no matter L is, while the ICRT algorithm has essentially linear running time. This is because the number of searches for the common remainder is L for the proposed algorithm. However, the number of searches is M/ϵ for the ICRT method, where ϵ is the searching step. In this simulation, the searching step is ϵ=1, and M=100.Figure 2
Total running time versus remainders number.Figure3 gives the robust estimation probability of the common remainder for searching methods when ϵ is 10, 1, and 0.1, respectively. It shows that the estimate of the common remainder can achieve an approximation of the optimal estimate whenϵis small enough. Clearly, decreasing the searching step can improve the performance at the cost of increasing the computation amount.Figure 3
Robust estimation probability of common remainder versus SNR.


## 5.  Conclusions

In this paper, we have proposed an efficient algorithm based on the robust CRT to perform single frequency estimation from multiple undersampled waveforms. We transformed the problem into the robust CRT firstly. Then we presented the optimal estimate of the common remainder and the original frequency. Based on the proposed algorithm, we presented the sufficient and necessary condition of the robust estimation. We finally applied the efficient algorithm to estimate the original frequency from multiple undersampled waveforms. Compared with the improved robust CRT algorithm, it has nearly the same performance but less computation. After using the undersampled waveforms, the sampling rate reduced sharply.



---

*Source: 413431-2013-12-03.xml*