# Digital Instantaneous Frequency Measurement of a Real Sinusoid Based on Three Sub-Nyquist Sampling Channels

Authors: Yu Su; Defu Jiang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Multichannel sub-Nyquist sampling is an efficient technique to break through the limitation of the Nyquist sampling theorem for the wideband digital instantaneous frequency measurement (DIFM) receiver. The significant challenge is calculating the folding frequency and solving the ambiguity quickly and accurately. Usually, the researchers adopt a fast Fourier transform (FFT) to calculate the folding frequency and a Chinese remainder theorem (CRT) to solve the ambiguity caused by undersampling. However, these algorithms have the drawback of long response time. In this paper, we use a frequency deduction algorithm based on the total least-squares estimation to measure the folding frequency accurately within a few clocks. We propose a frequency-band division algorithm by dividing the measurable frequency range into a few sub-bands to solve the frequency ambiguity. This deblurring algorithm is more efficient and faster than CRT. We analyse the performance of our algorithms by numerical simulations and functional hardware simulations and compare them with FFT and CRT. We also conduct experiments on a hardware platform to confirm the effectiveness of our algorithms. The simulation results show that our algorithms have better performance than FFT and CRT in accuracy and response time. Board-level measurement results verify that the proposed DIFM technique can capture the frequency from 0 to 5,000 MHz with a maximum of 0.8 MHz root-mean-squared error in less than 200 ns.


---

## Body

## 1. Introduction

Instantaneous frequency measurement (IFM) receiver plays an important role in military and civilian applications, such as electronic countermeasures (ECM), electronic counter-countermeasures, radar warning systems, and collision avoidance radars [1]. Usually, the IFM receiver needs to have characteristics including wide instantaneous bandwidth (IBW), low latency, high dynamic ranges, high-frequency measurement accuracy, and low hardware complexity.Traditional IFM receivers use analog technology to estimate frequency using multiple delay line channels and correlators [2, 3]. The analog delay line has a fixed time constant and is difficult to change in wideband applications. To improve IBW, multiple IFM channels must be used, which increase hardware complexity and cost. Furthermore, analog components are susceptible to environmental factors such as temperature, electromagnetic noise, and interference.In recent decades, with the development of high-speed analog to digital converter (ADC), field-programmable gate array (FPGA), and photonic devices, many researchers have developed various digital instantaneous frequency measurement (DIFM) technologies [4–14]. DFT algorithm is a simple algorithm widely used in frequency estimation [4–8]. In general, researchers focus on developing a variety of interpolation algorithms [4–6] or approximate models [7, 8] to mitigate the peak location errors when the sampling frequency is inconsistent with the signal frequency. However, to achieve satisfactory accuracy, the DFT technique must sample enough points, which results in longer processing time. Besides, multiple sampling channels must be used for large bandwidth. This will increase hardware complexity and cost. Monobit ADC sampling [9, 10] and photonic technology [11–13] can achieve a large bandwidth of tens of GHz owing to their ultrahigh sampling rate. However, the drawbacks are obvious, with errors of up to several hundred megahertz. Meanwhile, the low instantaneous dynamic range of monobit ADC also limits its application range. In [14], a practical DIFM technique is proposed, which adopts high-speed ADC and FPGA. It can measure the frequency with high accuracy in less than 180 ns. However, its IBW cannot exceed half of the sampling rate, so it cannot meet the engineering requirements. In general, all the DIFM techniques described above are based on the Nyquist sampling theorem. These technologies cannot deliver superior performance in terms of IBW, response time, precision, and hardware complexity at the same time.The multichannel sub-Nyquist sampling technique is an efficient way to break through the limitation of the Nyquist sampling theorem, thus achieving an ultrawide IBW [15]. However, the main challenge is to accurately calculate the frequency fk (k is the channel number, k = 1, 2, 3, …, K, K is the number of channels) folding into the first Nyquist region and solve the frequency ambiguity in a short time. Usually, a fast Fourier transform (FFT) [16, 17] algorithm is used to analyse the frequency to obtain high accuracy. However, the FFT algorithm has obvious shortcomings in processing time. Chinese remainder theorem (CRT), as an effective method to solve the ambiguity, has been widely used in frequency determination [16, 18–20]. In [16], a searching-based robust CRT is proposed for frequency estimation. Compared with conventional CRT, it reduces the number of searches and thus the amount of computation. In [21], a closed-form robust CRT theorem is proposed. It does not need to search for numbers and has better computational performance than [16]. However, this algorithm needs to perform integer and modulus calculations, which are not easy for the real-time hardware platform. In [22, 23], the authors proposed optimized algorithms to achieve better root-mean-squared error (RMSE) performance. But, both of them are computation intensive. Furthermore, the CRT model built for frequency estimation in [16, 18–20] ignores that the remainders may be folding frequency fk or image frequency fsk-fk (fsk is the sampling rate of the kth channel). In general, most researchers analyse CRT theoretically for optimal performance, with little attention paid to engineering applications. So, CRT is not suitable for the IFM application due to huge computation.In this paper, we implement a DIFM receiver on a hardware platform by combining frequency deduction and frequency-band division deblurring algorithms. The frequency deduction algorithm can measure the folding frequency accurately within a few clocks when combining with the total least-squares (TLS) estimation algorithm. We establish a relationship model between carrier frequencyfc and fk for the first time. Based on this model, we propose a measurable frequency-band division algorithm to solve the problem frequency ambiguity caused by undersampling. This deblurring algorithm does not require that all the sampling rates are coprime as CRT. So, it is convenient for engineering. The performances of frequency deduction and frequency-band division deblurring algorithms are compared with FFT and CRT [16, 21], respectively, by numerical simulations and functional hardware simulations. The simulation results show that compared with FFT, the frequency deduction algorithm can measure the frequency faster with the same accuracy. Also, the proposed frequency-band division deblurring algorithm is simpler, more accurate, faster, and more suitable for IFM application in a limited frequency range. We have conducted experiments on a hardware platform, and the results show that the DIFM technique which is based on the proposed algorithms can measure the carrier frequency from 0 to 5,000 MHz with a maximum of 0.8 MHz RMSE in less than 200 ns. Besides, the proposed deblurring algorithm can also be applied in phase unwrapping [24], DOA estimation [25], and SAR [26] as CRT.The remainder of this paper is organized as follows. In Section2, we will briefly describe the sub-Nyquist sampling and then present details of the frequency deduction and deblurring algorithm. In Section 3, the numerical results obtained by processing simulated and real data are presented. We end with concluding remarks in Section 4.


## 2. DIFM Algorithm

### 2.1. Sub-Nyquist Sampling

When the input signal is sampled by an ADC with sample frequencyfsk, we can obtain the frequency folding into the first Nyquist zone. As shown in Figure 1, the folding frequency varies with the frequency of the input signal periodically. This can be mathematically expressed as(1)fk=fc−dfsk,dfsk<fc≤d+12fsk,d+1fsk−fc,d+12fsk<fc≤d+1fsk,where d is a nonnegative integer.Figure 1
Output spectrum of three samplers (fs1 < fs2 < fs3).For each channel, iffc∈dfsk,d+1/2fsk, fk is the remainder frequency as described in CRT; otherwise, fk is the image frequency. Obviously, fk corresponds to multiple fc, so the carrier frequency cannot be uniquely determined by one sampling channel. So, the DIFM problem based on multiple channel sub-Nyquist sampling can be converted into two parts, folding frequency measurement and deblurring.


### 2.2. Frequency Deduction Based on the TLS Estimation Algorithm

Assume a monochromatic sinusoidal signal of thekth channel which is defined as(2)skt=Acos2πfkt+φ+ωt,where A is the signal amplitude, φ is the initial phase, ωt∼ℕ0,σ2, and ℕ0,σ2 denotes a Gaussian distribution with a zero mean and variance σ2. The signal-to-noise ratio (SNR) of the input signal is defined as SNR=A2/2σ2.If the sample period isTsk Tsk=1/fsk, then(3)skn=Acos2πfknTsk+φ+ωn,n=1,2,…,N,where N is the length of the signal.When considering the idea condition without noise, i.e.,ωn=0, it is easy to prove that(4)skn+1+skn−1=2skn⋅P,n=2,3,…,N−1,where P=cos2πfkTsk is a constant corresponding to fk one by one, P≤1.Therefore, the calculation offk is converted to the calculation of P. For any value of n, as long as there is no additional noise in the input signal, P can be accurately calculated as follows:(5)P=skn+1+skn−12skn.When the input signal contains noise termωn, then (4) is invalid and P calculated by (5) deviates from the true value. When the noise term is included, the estimated value of P can be expressed as(6)P^=Acos2πfkn+1Tsk+φ+Acos2πfkn−1Tsk+φ+ωn+1+ωn−12Acos2πfknTsk+φ+2ωn=2Acos2πfknTsk+φcos2πfkTsk+ωn+1+ωn−12Acos2πfknTsk+φ+2ωn.The errorΔP can be expressed as(7)ΔP=P−P^=ωn+1+ωn−1−2Pωn2Acos2πfknTsk+2ωn.Dividing the numerator and denominator of the right side of (7) by 2σ, we obtain(8)ΔP=ωn+1/2σ+ωn−1/2σ−2Pωn/2σA/σcos2πfknTsk+ωn/σ=ΗΤ,where Η=ωn+1/2σ+ωn−1/2σ−2Pωn/2σandΤ=A/σcos2πfknTsk+ωn/σ denote the numerator and denominator of the fraction, respectively.Sinceωn∼ℕ0,σ2, we get to know(9)Η∼ℕ0,0.5+P2,Τ∼ℕAσcos2πfknTsk,1.It is obvious that when the expectation ofΗ and the variance of Τ are constant, the variance of ΔP is proportional to the variance of Η and inversely proportional to the expectation of Τ. So, we conclude that the RMSE value of P^ is proportional to P and inversely proportional to 2SNRcos2πfknTsk.The estimation ofP using (6) is based on an assumption that the denominator skn is nonzero. However, in practice, due to the noise and quantization level error, we cannot ensure skn is nonzero. Also, if skn is close to zero, the estimation result will be significantly affected by noise.To solve the problem and estimateP correctly, we construct a linear mathematical model with L sampling points as(10)Y=PX,where Y=y1y2⋯yLT,X=x1x2⋯xLT, yn=skn+1+skn−1,andxn=2skn.The problem of estimating the value ofP is converted to fit an optimal straight line that passes through the origin and minimizes the total errors of L points. Usually, the least-squares (LS) method is adopted to estimate P. The LS method involves finding P^LS such that(11)Y−P^LSX2=ΔY2minimum,where ΔY is an L vector representing the errors in Y and ⋅2 is the l2 norm given by(12)ΔY2=ΔYTΔY=∑n=1LΔ2yn1/2.Then, the optimal fitting forP using the LS method is(13)P^LLS=XTX−1XTY=∑n=1Lxnyn∑n=1Lxn2=UV,where U and V are defined as(14)U=1N∑n=1Lynxn,V=1N∑n=1Lxn2.However, since bothY and X are contaminated by noise, the conventional LS method only considers the error of Y, and we adopt the TLS method to estimate P. The error of the TLS method takes into account both the Y errors and the X errors. So, the TLS method involves finding P^LTLS such that(15)ΔYΔXFminimum,where ⋅F denotes the Frobenius norm given by(16)ΔYΔXF=ΔYTΔY+ΔXTΔX=∑n=1LΔ2yn+Δ2xn1/2.The direct calculation of (14) is complicated. In [27], the authors have shown in this case of (10) that the error in TLS can be equivalent to perpendicular distance from a point to a line as follows:(17)enTLS=yn−P^LTLSxn1+P^LTLS2.The sum of the squared errors ofL points is(18)fP^LTLS=∑n=1LenTLS2=11+P^LTLS2∑n=1Lyn−P^LTLSxn2=11+P^LTLS2∑n=1Lyn2−2P^LTLS1+P^LTLS2∑n=1Lynxn+P^LTLS21+P^LTLS2∑n=1Lxn2=N⋅W1+P^LTLS2−U⋅2P^LTLS1+P^LTLS2+V⋅P^LTLS21+P^LTLS2,where W is defined as(19)W=1N∑n=1Lyn2.Settingf′P^LTLS=0 (′ represents differential operation), we obtain(20)U⋅PLTLS2+V−W⋅P^LTLS−U=0.We get the value ofP^LTLS that makes fP^LTLS the minimal value as(21)P^LTLS=4U2+V−W2+W−V2U.Obviously,U and V can hardly be zero especially for large L. Through analysis of (13) and (21), we can see both LS and TLS methods avoid the problem of zero denominators and fit the line well. However, in the LS problem, it is the vertical distances that are critical, while in the TLS problem, it is the perpendicular distances that are important. The geometric interpretation of the LS problem and the TLS problem is depicted in Figure 2. So, the TLS method is better than LS with respect to both Y errors and X errors in the curve fitting.Figure 2
LS and TLS.In addition, we note that estimatingP from (6) is a special case for LS and TLS methods in the case of L = 1.We can get the folding frequency as follows:(22)f^k=fsk2πcos−1P^.Coordinated rotation digital computer (CORDIC) algorithm or look-up-table can be used to calculatecos−1P^. In this paper, we select look-up-table to compute the folding frequency for less response time at the expense of RAM resources. To further reduce the error of f^k, we have done moving average of S points at the expense of increasing S clock response time.The frequency error is defined asΔfk=f^k−fk, which is influenced by P, 2SNRcos2πfknTsk, L, SNR, S, and the quantization level of look-up-table. For purposes of analysis, we consider Δfk being uniformly distributed over a certain range.The proposed method for estimating folded frequency is simple, fast, and accurate. However, it is important to realize that the derivation is only valid for the monochromatic signal since equations (10)–(22) do not hold in the case of multiple sinusoidal signals.


### 2.3. Frequency-Band Division Deblurring

Based on the Nyquist sampling theorem, only the frequency underfs/2 can be exactly recovered without aliasing. Ambiguity occurs when the frequency is extended to more than fs/2. To reconstruct the unambiguous frequency, we employ a multichannel sampling technique. (1) can also be expressed as(23)fc=mkfsk+bkfk,0≤fk<fsk2,mk=0,1,2,…,Mk+1,bk=±1,where mk is an unknown integer, bk = 1 when fk is the remainder frequency, bk = –1 when fk is the image frequency, and Mk is the integer part of fup/fsk, where fup is the upper bound of measurable frequency. Usually, Mk will not be too large as the limitation of full power bandwidth (FPBW) of ADC. In [28], the authors proposed and demonstrated a theorem that gives the relationship between the sampling rates and the maximum unambiguous frequency. The maximum unambiguous frequency determined by K channels is derived and can be expressed as(24)fmax=minI,JLCMfs1G1,…,fsIG1+LCMfs1G2,…,fsJG22,I+J=K,I,J=1,…,K,fs1,fs2,…,fsK=fs1G1,…,fsIG1∪fs1G2,…,fsJG2,where LCM denotes least common multiple.IfK = 2, it is easy to find that fmax is less than the maximum sampling rate of ADC and cannot meet the application requirements. In this paper, we define K = 3. In this case, fmax = {LCM (fs1, fs2) + fs3}/2 (fs1 < fs2 < fs3). Proper setting of sampling rates of three channels can guarantee that the value of fmax is larger than the FPBW of ADC. Therefore, all the frequencies that fall into the FPBW can be uniquely reconstructed.From (23), the reconstruction of fc is transformed into the calculation of fk, mk, and bk. The calculation of fk has been introduced in Section 2.2. In this section, we study the fast algorithm for determining mk and bk from f^k. We divide the measurable frequency range into Q sub-bands as shown in Figure 1. Each frequency sub-band corresponds uniquely to a group of m1,m2,m3,b1,b2,b3. The values of mk and bk can be calculated in advance as(25)mk=fcfsk,rk<fsk2,fcfsk+1,rk≥fsk2,bk=1,rk<fsk2,−1,rk≥fsk2,where ⌊⋅⌋ denotes the flooring operation and rk is the remainder of fc modulo fsk.Table1 presents an example for fs1 = 1500, fs2 = 1600, and fs3 = 1700. Usually, Q will not be too large since Mk is small. So, we do not need to search for too many numbers. For example, if fc ∈ [0, 2fs3), Q = 12. So, the frequency range can be divided into 12 sub-bands. We just have to judge the sub-band which the carrier frequency falls into.Table 1
Parameters corresponding to frequency sub-bands.
Number
f
c
Parameters
m
1
m
2
m
3
b
1
b
2
b
3
1
[0,fs1/2)
0
0
0
1
1
1
2
[fs1/2, fs2/2)
0
0
1
1
1
−1
3
[fs2/2, fs3/2)
0
1
1
1
−1
−1
4
[fs3/2, fs1)
1
1
1
−1
−1
−1
5
[fs1, fs2)
1
1
1
−1
−1
1
6
[fs2, fs3)
1
1
1
−1
1
1
7
[fs3, 3fs1/2)
1
1
1
1
1
1
8
[3fs1/2, 3fs2/2)
1
1
2
1
1
−1
9
[3fs2/2, 3fs3/2)
1
2
2
1
−1
−1
...
...
...
...
...
...
...
...The optimal sub-band is determined by calculating the error of each sub-band. We list all the values ofm1,m2,m3,b1,b2,b3 and carry out subtraction of m1fs1+b1f^1, m2fs2+b2f^2, and m3fs3+b3f^3. The sum of error between arbitrary two channels can be described as(26)ferrori=m1ifs1+b1if^1−m2ifs2+b2if^2+m1ifs1+b1if^1−m3ifs3+b3if^3+m2ifs2+b2if^2−m3ifs3+b3if^3,i=1,2,…,Q,where ferrori is the total error of three channels and i is the sub-band index.AllQ data are compared by the pipeline technique to get the minimum value. This process can be expressed as follows:(1)
DivideQ data into Q/2 groups and do a pairwise comparison among all Q data.(2)
Choose smallerQ/2 from the results of (1), and continue the next comparison. This processing procedure will continue until minimum ferrorl is achieved.(3)
Register the values ofm1l,m2l,m3l,b1l,b2l,b3l that give the least value of ferrorl.When the right sub-band is determined,fc is calculated using the average of three channels as(27)f^c=m1lfs1+b1lf^1+m2lfs2+b2lf^2+m3lfs3+b3lf^33.For all carrier frequencies, we can obtain(28)fc=m1lfs1+b1lf1=m2lfs2+b2lf2=m3lfs3+b3lf3.Therefore, the error of carrier frequency measurement is(29)Δfc=f^c−fc=m1lfs1+b1lf^1+m2lfs2+b2lf^2+m3lfs3+b3lf^33−fc=Δf1+Δf2+Δf33.Assume thatΔfk is uniformly distributed between −τ,τ, τ is the error bound. The variance of Δfk is σ2=τ2/3. Only if the sub-band is determined correctly, we can ensure that Δfc≤τ.It should be noted that the proposed deblurring method is still valid in the case of multiple sinusoidal signals. However, the computation complexity isT3 times that of the case of a single sinusoidal signal. Here, T is the number of signals. For example, if T = 2, then there are two folded frequencies in each channel. Therefore, there are 8 combinations of the three channels. Since the proposed deblurring method requires only simple addition and search operations, the computational complexity is acceptable.


## 2.1. Sub-Nyquist Sampling

When the input signal is sampled by an ADC with sample frequencyfsk, we can obtain the frequency folding into the first Nyquist zone. As shown in Figure 1, the folding frequency varies with the frequency of the input signal periodically. This can be mathematically expressed as(1)fk=fc−dfsk,dfsk<fc≤d+12fsk,d+1fsk−fc,d+12fsk<fc≤d+1fsk,where d is a nonnegative integer.Figure 1
Output spectrum of three samplers (fs1 < fs2 < fs3).For each channel, iffc∈dfsk,d+1/2fsk, fk is the remainder frequency as described in CRT; otherwise, fk is the image frequency. Obviously, fk corresponds to multiple fc, so the carrier frequency cannot be uniquely determined by one sampling channel. So, the DIFM problem based on multiple channel sub-Nyquist sampling can be converted into two parts, folding frequency measurement and deblurring.


## 2.2. Frequency Deduction Based on the TLS Estimation Algorithm

Assume a monochromatic sinusoidal signal of thekth channel which is defined as(2)skt=Acos2πfkt+φ+ωt,where A is the signal amplitude, φ is the initial phase, ωt∼ℕ0,σ2, and ℕ0,σ2 denotes a Gaussian distribution with a zero mean and variance σ2. The signal-to-noise ratio (SNR) of the input signal is defined as SNR=A2/2σ2.If the sample period isTsk Tsk=1/fsk, then(3)skn=Acos2πfknTsk+φ+ωn,n=1,2,…,N,where N is the length of the signal.When considering the idea condition without noise, i.e.,ωn=0, it is easy to prove that(4)skn+1+skn−1=2skn⋅P,n=2,3,…,N−1,where P=cos2πfkTsk is a constant corresponding to fk one by one, P≤1.Therefore, the calculation offk is converted to the calculation of P. For any value of n, as long as there is no additional noise in the input signal, P can be accurately calculated as follows:(5)P=skn+1+skn−12skn.When the input signal contains noise termωn, then (4) is invalid and P calculated by (5) deviates from the true value. When the noise term is included, the estimated value of P can be expressed as(6)P^=Acos2πfkn+1Tsk+φ+Acos2πfkn−1Tsk+φ+ωn+1+ωn−12Acos2πfknTsk+φ+2ωn=2Acos2πfknTsk+φcos2πfkTsk+ωn+1+ωn−12Acos2πfknTsk+φ+2ωn.The errorΔP can be expressed as(7)ΔP=P−P^=ωn+1+ωn−1−2Pωn2Acos2πfknTsk+2ωn.Dividing the numerator and denominator of the right side of (7) by 2σ, we obtain(8)ΔP=ωn+1/2σ+ωn−1/2σ−2Pωn/2σA/σcos2πfknTsk+ωn/σ=ΗΤ,where Η=ωn+1/2σ+ωn−1/2σ−2Pωn/2σandΤ=A/σcos2πfknTsk+ωn/σ denote the numerator and denominator of the fraction, respectively.Sinceωn∼ℕ0,σ2, we get to know(9)Η∼ℕ0,0.5+P2,Τ∼ℕAσcos2πfknTsk,1.It is obvious that when the expectation ofΗ and the variance of Τ are constant, the variance of ΔP is proportional to the variance of Η and inversely proportional to the expectation of Τ. So, we conclude that the RMSE value of P^ is proportional to P and inversely proportional to 2SNRcos2πfknTsk.The estimation ofP using (6) is based on an assumption that the denominator skn is nonzero. However, in practice, due to the noise and quantization level error, we cannot ensure skn is nonzero. Also, if skn is close to zero, the estimation result will be significantly affected by noise.To solve the problem and estimateP correctly, we construct a linear mathematical model with L sampling points as(10)Y=PX,where Y=y1y2⋯yLT,X=x1x2⋯xLT, yn=skn+1+skn−1,andxn=2skn.The problem of estimating the value ofP is converted to fit an optimal straight line that passes through the origin and minimizes the total errors of L points. Usually, the least-squares (LS) method is adopted to estimate P. The LS method involves finding P^LS such that(11)Y−P^LSX2=ΔY2minimum,where ΔY is an L vector representing the errors in Y and ⋅2 is the l2 norm given by(12)ΔY2=ΔYTΔY=∑n=1LΔ2yn1/2.Then, the optimal fitting forP using the LS method is(13)P^LLS=XTX−1XTY=∑n=1Lxnyn∑n=1Lxn2=UV,where U and V are defined as(14)U=1N∑n=1Lynxn,V=1N∑n=1Lxn2.However, since bothY and X are contaminated by noise, the conventional LS method only considers the error of Y, and we adopt the TLS method to estimate P. The error of the TLS method takes into account both the Y errors and the X errors. So, the TLS method involves finding P^LTLS such that(15)ΔYΔXFminimum,where ⋅F denotes the Frobenius norm given by(16)ΔYΔXF=ΔYTΔY+ΔXTΔX=∑n=1LΔ2yn+Δ2xn1/2.The direct calculation of (14) is complicated. In [27], the authors have shown in this case of (10) that the error in TLS can be equivalent to perpendicular distance from a point to a line as follows:(17)enTLS=yn−P^LTLSxn1+P^LTLS2.The sum of the squared errors ofL points is(18)fP^LTLS=∑n=1LenTLS2=11+P^LTLS2∑n=1Lyn−P^LTLSxn2=11+P^LTLS2∑n=1Lyn2−2P^LTLS1+P^LTLS2∑n=1Lynxn+P^LTLS21+P^LTLS2∑n=1Lxn2=N⋅W1+P^LTLS2−U⋅2P^LTLS1+P^LTLS2+V⋅P^LTLS21+P^LTLS2,where W is defined as(19)W=1N∑n=1Lyn2.Settingf′P^LTLS=0 (′ represents differential operation), we obtain(20)U⋅PLTLS2+V−W⋅P^LTLS−U=0.We get the value ofP^LTLS that makes fP^LTLS the minimal value as(21)P^LTLS=4U2+V−W2+W−V2U.Obviously,U and V can hardly be zero especially for large L. Through analysis of (13) and (21), we can see both LS and TLS methods avoid the problem of zero denominators and fit the line well. However, in the LS problem, it is the vertical distances that are critical, while in the TLS problem, it is the perpendicular distances that are important. The geometric interpretation of the LS problem and the TLS problem is depicted in Figure 2. So, the TLS method is better than LS with respect to both Y errors and X errors in the curve fitting.Figure 2
LS and TLS.In addition, we note that estimatingP from (6) is a special case for LS and TLS methods in the case of L = 1.We can get the folding frequency as follows:(22)f^k=fsk2πcos−1P^.Coordinated rotation digital computer (CORDIC) algorithm or look-up-table can be used to calculatecos−1P^. In this paper, we select look-up-table to compute the folding frequency for less response time at the expense of RAM resources. To further reduce the error of f^k, we have done moving average of S points at the expense of increasing S clock response time.The frequency error is defined asΔfk=f^k−fk, which is influenced by P, 2SNRcos2πfknTsk, L, SNR, S, and the quantization level of look-up-table. For purposes of analysis, we consider Δfk being uniformly distributed over a certain range.The proposed method for estimating folded frequency is simple, fast, and accurate. However, it is important to realize that the derivation is only valid for the monochromatic signal since equations (10)–(22) do not hold in the case of multiple sinusoidal signals.


## 2.3. Frequency-Band Division Deblurring

Based on the Nyquist sampling theorem, only the frequency underfs/2 can be exactly recovered without aliasing. Ambiguity occurs when the frequency is extended to more than fs/2. To reconstruct the unambiguous frequency, we employ a multichannel sampling technique. (1) can also be expressed as(23)fc=mkfsk+bkfk,0≤fk<fsk2,mk=0,1,2,…,Mk+1,bk=±1,where mk is an unknown integer, bk = 1 when fk is the remainder frequency, bk = –1 when fk is the image frequency, and Mk is the integer part of fup/fsk, where fup is the upper bound of measurable frequency. Usually, Mk will not be too large as the limitation of full power bandwidth (FPBW) of ADC. In [28], the authors proposed and demonstrated a theorem that gives the relationship between the sampling rates and the maximum unambiguous frequency. The maximum unambiguous frequency determined by K channels is derived and can be expressed as(24)fmax=minI,JLCMfs1G1,…,fsIG1+LCMfs1G2,…,fsJG22,I+J=K,I,J=1,…,K,fs1,fs2,…,fsK=fs1G1,…,fsIG1∪fs1G2,…,fsJG2,where LCM denotes least common multiple.IfK = 2, it is easy to find that fmax is less than the maximum sampling rate of ADC and cannot meet the application requirements. In this paper, we define K = 3. In this case, fmax = {LCM (fs1, fs2) + fs3}/2 (fs1 < fs2 < fs3). Proper setting of sampling rates of three channels can guarantee that the value of fmax is larger than the FPBW of ADC. Therefore, all the frequencies that fall into the FPBW can be uniquely reconstructed.From (23), the reconstruction of fc is transformed into the calculation of fk, mk, and bk. The calculation of fk has been introduced in Section 2.2. In this section, we study the fast algorithm for determining mk and bk from f^k. We divide the measurable frequency range into Q sub-bands as shown in Figure 1. Each frequency sub-band corresponds uniquely to a group of m1,m2,m3,b1,b2,b3. The values of mk and bk can be calculated in advance as(25)mk=fcfsk,rk<fsk2,fcfsk+1,rk≥fsk2,bk=1,rk<fsk2,−1,rk≥fsk2,where ⌊⋅⌋ denotes the flooring operation and rk is the remainder of fc modulo fsk.Table1 presents an example for fs1 = 1500, fs2 = 1600, and fs3 = 1700. Usually, Q will not be too large since Mk is small. So, we do not need to search for too many numbers. For example, if fc ∈ [0, 2fs3), Q = 12. So, the frequency range can be divided into 12 sub-bands. We just have to judge the sub-band which the carrier frequency falls into.Table 1
Parameters corresponding to frequency sub-bands.
Number
f
c
Parameters
m
1
m
2
m
3
b
1
b
2
b
3
1
[0,fs1/2)
0
0
0
1
1
1
2
[fs1/2, fs2/2)
0
0
1
1
1
−1
3
[fs2/2, fs3/2)
0
1
1
1
−1
−1
4
[fs3/2, fs1)
1
1
1
−1
−1
−1
5
[fs1, fs2)
1
1
1
−1
−1
1
6
[fs2, fs3)
1
1
1
−1
1
1
7
[fs3, 3fs1/2)
1
1
1
1
1
1
8
[3fs1/2, 3fs2/2)
1
1
2
1
1
−1
9
[3fs2/2, 3fs3/2)
1
2
2
1
−1
−1
...
...
...
...
...
...
...
...The optimal sub-band is determined by calculating the error of each sub-band. We list all the values ofm1,m2,m3,b1,b2,b3 and carry out subtraction of m1fs1+b1f^1, m2fs2+b2f^2, and m3fs3+b3f^3. The sum of error between arbitrary two channels can be described as(26)ferrori=m1ifs1+b1if^1−m2ifs2+b2if^2+m1ifs1+b1if^1−m3ifs3+b3if^3+m2ifs2+b2if^2−m3ifs3+b3if^3,i=1,2,…,Q,where ferrori is the total error of three channels and i is the sub-band index.AllQ data are compared by the pipeline technique to get the minimum value. This process can be expressed as follows:(1)
DivideQ data into Q/2 groups and do a pairwise comparison among all Q data.(2)
Choose smallerQ/2 from the results of (1), and continue the next comparison. This processing procedure will continue until minimum ferrorl is achieved.(3)
Register the values ofm1l,m2l,m3l,b1l,b2l,b3l that give the least value of ferrorl.When the right sub-band is determined,fc is calculated using the average of three channels as(27)f^c=m1lfs1+b1lf^1+m2lfs2+b2lf^2+m3lfs3+b3lf^33.For all carrier frequencies, we can obtain(28)fc=m1lfs1+b1lf1=m2lfs2+b2lf2=m3lfs3+b3lf3.Therefore, the error of carrier frequency measurement is(29)Δfc=f^c−fc=m1lfs1+b1lf^1+m2lfs2+b2lf^2+m3lfs3+b3lf^33−fc=Δf1+Δf2+Δf33.Assume thatΔfk is uniformly distributed between −τ,τ, τ is the error bound. The variance of Δfk is σ2=τ2/3. Only if the sub-band is determined correctly, we can ensure that Δfc≤τ.It should be noted that the proposed deblurring method is still valid in the case of multiple sinusoidal signals. However, the computation complexity isT3 times that of the case of a single sinusoidal signal. Here, T is the number of signals. For example, if T = 2, then there are two folded frequencies in each channel. Therefore, there are 8 combinations of the three channels. Since the proposed deblurring method requires only simple addition and search operations, the computational complexity is acceptable.


## 3. Results and Discussion

### 3.1. Simulation

#### 3.1.1. Simulation for Frequency Deduction

The performance of the frequency deduction algorithm is verified by a single channel frequency simulation. The configuration of simulation isfs = 1600 MHz and pulse width = 20 us.Figures3(a)–3(c) show the change rule of RMSE value of P^ to the input frequency, time index, and SNR in different values of L. It can be seen that the RMSE value of P^ is proportional to P and inversely proportional to cos2πfknTs and SNR, which validates our theory. A few sampling points for both LS and TLS methods can significantly weaken the influence caused by P, cos2πfknTs, and SNR. Also, the TLS method is better than the LS method which also verifies our theory. Figure 3(d) analyses the influence of sampling points and the quantitative level to our algorithm under different S. It can be seen that a few points of moving average can significantly decrease the RMSE value of f^k. Also, a 12 ∗ 12 bit look-up-table is a good choice for a trade-off between accuracy and resources.Figure 3
Performance analysis for frequency deduction: (a) RMSE ofP^ versus input frequency with SNR = 20 dB; (b) RMSE of P^ versus n with SNR = 20 dB and fk = 50 MHz; (c) RMSE of P^ versus sample SNR with fk = 50 MHz; (d) RMS error of f^k versus size of look-up-table with SNR = 20 dB, TLS L = 8, and fk = 50 MHz.
(a)
(b)
(c)
(d)Functional hardware simulations have been conducted to compare the performance of the frequency deduction algorithm and the FFT algorithm. The input frequency is an arbitrary value between 0 and 800 MHz. The input data are divided into 8 phases, and the operation clock isfclock = fs/8 = 200 MHz. The latency is simulated by Quartus II and Modelsim. The RMSE is simulated by Matlab using 12,000 trials. The results are shown in Table 2. We can conclude that the frequency deduction algorithm performs better performance in accuracy and response time.Table 2
Implementation results for FFT and frequency deduction algorithms.
SNR (dB)
L
Comparison between FFT and frequency deduction
FFT
Frequency deduction
N_FFT
RMSE (MHz)
Latency (ns)
X
RMSE (MHz)
Latency (ns)
10
8
1024
15.773
10790
4
10.239
105
2048
15.735
21080
8
8.477
125
16
1024
5.269
10795
4
6.831
110
2048
5.181
21085
8
6.237
130
20
8
1024
9.500
10790
4
2.262
105
2048
9.522
21080
8
1.670
125
16
1024
3.161
10795
4
1.219
110
2048
3.039
21085
8
0.967
130
30
8
1024
8.665
10790
4
0.676
105
2048
8.689
21080
8
0.487
125
16
1024
2.847
10795
4
0.339
110
2048
2.739
21085
8
0.251
130Table3 compares the computation complexity for frequency deduction and FFT with the parameters N_FFT = 1024 and L = 16. For the FFT method, its complexity can be found easily from [6]. The N_FFT-point FFT requires to calculate 2 × N_FFT × log2N_FFT real multiplications, 2 × N_FFT × log2N_FFT real additions, and N_FFT-point peak search. The complexity of frequency deduction is derived from (8) to (11) in this paper, which needs 3L + 4 multiplications, 3(L − 1) + 2 additions, 1 division, and some other operations. It is obvious that the proposed method requires far fewer resources than FFT.Table 3
Computation complexity comparison of the FFT and frequency deduction.
Method
MULT.
ADD.
DIV.
Search numbers
Square root
Look-up-table
FFT
20480
20480
0
1024
0
0
Frequency deduction
52
47
1
0
1
1
MULT. represents multiplication operation, ADD. represents addition operation, and DIV. represents division operation.


#### 3.1.2. Performance Comparison between Frequency-Band Division and CRT

CRT described in [16, 21] applies only to the complex exponential signal, i.e., bk = 1. In the real signal case, the value of bk must be estimated in advance. In [29], a method to estimate bk by employing a triggering circuit to detect the waveform’s zero crossing point ‘O’ is proposed. The method needs extra hardware. For simplicity, we assume that the value of bk is estimated correctly for CRT. The configuration of simulation is K = 3, fs1=MΓ1=1500MHz, fs2=MΓ2=1600MHz, fs3=MΓ3=1700MHz, and fc is a random number and is uniformly distributed between 0 to 5,000 MHz. The folding frequency f^k is obtained by f^k=fk+Δfk, and fk is the accurate folding frequency. 12,000 trials are implemented for different SNR. The SNR is defined as SNR=20logM/σ which is the same as [21]. So, we can derive that(30)τ=3M⋅10−SNR/20,where M = 100 MHz is the greatest common divisor of fs1, fs2, and fs3.From the analysis of Section2.2, we can see that accurate judgment of mk is quite important. In Figure 4(a), we simulated the probability for judging mk correctly in different SNR. The result shows that the proposed algorithm can judge the folding integers mk more accurately than CRT. Since both robust CRT and closed-form CRT need to calculate mk one by one, the calculation process is complicated especially in low SNR, and it is difficult to guarantee the calculation accuracy of each mk. However, the proposed algorithm judges mk by comparing the errors of each sub-band, so it is less sensitive to noise.Figure 4
Performance comparison between robust CRT, closed-form CRT, and the proposed algorithm for the real signal: (a) probability of correction versus SNR; (b) RMSE versus SNR; (c) mean error versus SNR.
(a)
(b)
(c)Figures4(b) and 4(c) show the curves of RMSE value and mean error versus SNR. The theoretical RMSE value of CRT is given by Wang and Xia [21] as(31)fc,RMSE=σK.Since our algorithm is less sensitive to noise, it has better RMSE and mean error performance at low SNR. Additionally, we can see that when the SNR is larger than the bound which from the error boundτ=M/4, RMSE and mean error decrease linearly with SNR due to the robustness of CRT. In this case, the bound is about 16.8 dB, indicated by the black line in figures. However, the bound of our proposed algorithm is about 11 dB, which is lower than the CRT. So, it shows better robustness.Table4 compares the computation complexity between the frequency-band division algorithm and robust CRT and closed-form CRT. Through careful analysis (17), (18), and (28) in [16], we find that the robust CRT method needs to calculate 4Γ3(K − 1) +3 (Γ2 + Γ1) multiplications, 4Γ3 (K − 1) +3 (Γ2+Γ1) divisions, 4Γ3(K − 1) +4 (Γ2 + Γ1) additions, and (K − 1) Γ1+Γ2+⋯+ΓK searches. Similarly, the closed-form CRT in [21] needs to calculate 4K multiplications, 3K divisions, 2K + 1 additions, K rounding of integers, and 2K modular operations. Since the values of mkfsk can be registered in advance, there is no need for multiplication in our algorithm. The frequency-band division deblurring method needs only 6Q additions, 3Q operations for finding absolute values, and Q searches. Therefore, we can conclude that the deblurring algorithm proposed is faster and more effective, so it is more suitable for the application of DIFM.Table 4
Computation complexity comparison of the CRT and frequency-band division.
Method
MULT.
DIV.
Mod.
Round.
ADD.
ABS.
Search numbers
Robust CRT [16]
205
205
0
0
236
0
63
Closed-form CRT [21]
8
6
4
2
5
0
0
Frequency-band division
0
0
0
0
108
54
18
Mod. represents a modular operation, Round. represents rounding integer operation, and ABS. represents absolute value operation.


### 3.2. Experiment

Experiments have also been conducted on a digital radio frequency memory (DRFM) module (Figure5(a)) to verify the practicability of the DIFM technique. This module is originally designed for ECM and can be reused in our experiments with minor modifications. Signal source 1 is used to generate sampling clocks of different frequencies for the ADC. Signal source 2 is used to generate a monochromatic pulsed signal with a frequency from 0 to 5,000 MHz and a step size of 100 MHz. After being converted to the differential signal through Anaren B0430J50100A00 Balun, the input signal is injected into National Semiconductor ADC12D1800 ADC. The ADC can digitize the input signal with a sampling clock of up to 1,800 MHz and a data bit width of 12 bits. In the experiment, we set the sampling rates to be the same as in the simulation. The sampling data are transmitted to the Altera Stratix4SGX360 FPGA. Due to the limitation of experimental conditions, we only sample the input signal of one channel at a time and calculate folding frequency. Then, a deblurring operation is completed by combining the three channel results. All the algorithms were completed in FPGA. Quartus II 13.1 firmware fitting and timing analysis software shows a maximum clock speed of 200 MHz and FPGA logic resource usage of less than 5% for the 4SGX360 device. We have also used 12,000 trials for implementation. Figure 5(b) shows the experimental results of the mean deviation, maximum absolute error, and RMS error. It demonstrates that the DIFM technique proposed in this paper can measure the frequency from 0 to 5,000 MHz with high accuracy. Also, it can be observed that RMS errors of more than 90% frequency points are less than 0.5 MHz. The latency of this implementation is simulated as 38 FPGA clocks, which at the given clock rate yields(32)TDIFM=38∗fclk=190ns.Figure 5
Experiment platform and results. (a) DIFM hardware for the implementation of DIFM. (b) Experiment results of DIFM.
(a)
(b)


## 3.1. Simulation

### 3.1.1. Simulation for Frequency Deduction

The performance of the frequency deduction algorithm is verified by a single channel frequency simulation. The configuration of simulation isfs = 1600 MHz and pulse width = 20 us.Figures3(a)–3(c) show the change rule of RMSE value of P^ to the input frequency, time index, and SNR in different values of L. It can be seen that the RMSE value of P^ is proportional to P and inversely proportional to cos2πfknTs and SNR, which validates our theory. A few sampling points for both LS and TLS methods can significantly weaken the influence caused by P, cos2πfknTs, and SNR. Also, the TLS method is better than the LS method which also verifies our theory. Figure 3(d) analyses the influence of sampling points and the quantitative level to our algorithm under different S. It can be seen that a few points of moving average can significantly decrease the RMSE value of f^k. Also, a 12 ∗ 12 bit look-up-table is a good choice for a trade-off between accuracy and resources.Figure 3
Performance analysis for frequency deduction: (a) RMSE ofP^ versus input frequency with SNR = 20 dB; (b) RMSE of P^ versus n with SNR = 20 dB and fk = 50 MHz; (c) RMSE of P^ versus sample SNR with fk = 50 MHz; (d) RMS error of f^k versus size of look-up-table with SNR = 20 dB, TLS L = 8, and fk = 50 MHz.
(a)
(b)
(c)
(d)Functional hardware simulations have been conducted to compare the performance of the frequency deduction algorithm and the FFT algorithm. The input frequency is an arbitrary value between 0 and 800 MHz. The input data are divided into 8 phases, and the operation clock isfclock = fs/8 = 200 MHz. The latency is simulated by Quartus II and Modelsim. The RMSE is simulated by Matlab using 12,000 trials. The results are shown in Table 2. We can conclude that the frequency deduction algorithm performs better performance in accuracy and response time.Table 2
Implementation results for FFT and frequency deduction algorithms.
SNR (dB)
L
Comparison between FFT and frequency deduction
FFT
Frequency deduction
N_FFT
RMSE (MHz)
Latency (ns)
X
RMSE (MHz)
Latency (ns)
10
8
1024
15.773
10790
4
10.239
105
2048
15.735
21080
8
8.477
125
16
1024
5.269
10795
4
6.831
110
2048
5.181
21085
8
6.237
130
20
8
1024
9.500
10790
4
2.262
105
2048
9.522
21080
8
1.670
125
16
1024
3.161
10795
4
1.219
110
2048
3.039
21085
8
0.967
130
30
8
1024
8.665
10790
4
0.676
105
2048
8.689
21080
8
0.487
125
16
1024
2.847
10795
4
0.339
110
2048
2.739
21085
8
0.251
130Table3 compares the computation complexity for frequency deduction and FFT with the parameters N_FFT = 1024 and L = 16. For the FFT method, its complexity can be found easily from [6]. The N_FFT-point FFT requires to calculate 2 × N_FFT × log2N_FFT real multiplications, 2 × N_FFT × log2N_FFT real additions, and N_FFT-point peak search. The complexity of frequency deduction is derived from (8) to (11) in this paper, which needs 3L + 4 multiplications, 3(L − 1) + 2 additions, 1 division, and some other operations. It is obvious that the proposed method requires far fewer resources than FFT.Table 3
Computation complexity comparison of the FFT and frequency deduction.
Method
MULT.
ADD.
DIV.
Search numbers
Square root
Look-up-table
FFT
20480
20480
0
1024
0
0
Frequency deduction
52
47
1
0
1
1
MULT. represents multiplication operation, ADD. represents addition operation, and DIV. represents division operation.


### 3.1.2. Performance Comparison between Frequency-Band Division and CRT

CRT described in [16, 21] applies only to the complex exponential signal, i.e., bk = 1. In the real signal case, the value of bk must be estimated in advance. In [29], a method to estimate bk by employing a triggering circuit to detect the waveform’s zero crossing point ‘O’ is proposed. The method needs extra hardware. For simplicity, we assume that the value of bk is estimated correctly for CRT. The configuration of simulation is K = 3, fs1=MΓ1=1500MHz, fs2=MΓ2=1600MHz, fs3=MΓ3=1700MHz, and fc is a random number and is uniformly distributed between 0 to 5,000 MHz. The folding frequency f^k is obtained by f^k=fk+Δfk, and fk is the accurate folding frequency. 12,000 trials are implemented for different SNR. The SNR is defined as SNR=20logM/σ which is the same as [21]. So, we can derive that(30)τ=3M⋅10−SNR/20,where M = 100 MHz is the greatest common divisor of fs1, fs2, and fs3.From the analysis of Section2.2, we can see that accurate judgment of mk is quite important. In Figure 4(a), we simulated the probability for judging mk correctly in different SNR. The result shows that the proposed algorithm can judge the folding integers mk more accurately than CRT. Since both robust CRT and closed-form CRT need to calculate mk one by one, the calculation process is complicated especially in low SNR, and it is difficult to guarantee the calculation accuracy of each mk. However, the proposed algorithm judges mk by comparing the errors of each sub-band, so it is less sensitive to noise.Figure 4
Performance comparison between robust CRT, closed-form CRT, and the proposed algorithm for the real signal: (a) probability of correction versus SNR; (b) RMSE versus SNR; (c) mean error versus SNR.
(a)
(b)
(c)Figures4(b) and 4(c) show the curves of RMSE value and mean error versus SNR. The theoretical RMSE value of CRT is given by Wang and Xia [21] as(31)fc,RMSE=σK.Since our algorithm is less sensitive to noise, it has better RMSE and mean error performance at low SNR. Additionally, we can see that when the SNR is larger than the bound which from the error boundτ=M/4, RMSE and mean error decrease linearly with SNR due to the robustness of CRT. In this case, the bound is about 16.8 dB, indicated by the black line in figures. However, the bound of our proposed algorithm is about 11 dB, which is lower than the CRT. So, it shows better robustness.Table4 compares the computation complexity between the frequency-band division algorithm and robust CRT and closed-form CRT. Through careful analysis (17), (18), and (28) in [16], we find that the robust CRT method needs to calculate 4Γ3(K − 1) +3 (Γ2 + Γ1) multiplications, 4Γ3 (K − 1) +3 (Γ2+Γ1) divisions, 4Γ3(K − 1) +4 (Γ2 + Γ1) additions, and (K − 1) Γ1+Γ2+⋯+ΓK searches. Similarly, the closed-form CRT in [21] needs to calculate 4K multiplications, 3K divisions, 2K + 1 additions, K rounding of integers, and 2K modular operations. Since the values of mkfsk can be registered in advance, there is no need for multiplication in our algorithm. The frequency-band division deblurring method needs only 6Q additions, 3Q operations for finding absolute values, and Q searches. Therefore, we can conclude that the deblurring algorithm proposed is faster and more effective, so it is more suitable for the application of DIFM.Table 4
Computation complexity comparison of the CRT and frequency-band division.
Method
MULT.
DIV.
Mod.
Round.
ADD.
ABS.
Search numbers
Robust CRT [16]
205
205
0
0
236
0
63
Closed-form CRT [21]
8
6
4
2
5
0
0
Frequency-band division
0
0
0
0
108
54
18
Mod. represents a modular operation, Round. represents rounding integer operation, and ABS. represents absolute value operation.


## 3.1.1. Simulation for Frequency Deduction

The performance of the frequency deduction algorithm is verified by a single channel frequency simulation. The configuration of simulation isfs = 1600 MHz and pulse width = 20 us.Figures3(a)–3(c) show the change rule of RMSE value of P^ to the input frequency, time index, and SNR in different values of L. It can be seen that the RMSE value of P^ is proportional to P and inversely proportional to cos2πfknTs and SNR, which validates our theory. A few sampling points for both LS and TLS methods can significantly weaken the influence caused by P, cos2πfknTs, and SNR. Also, the TLS method is better than the LS method which also verifies our theory. Figure 3(d) analyses the influence of sampling points and the quantitative level to our algorithm under different S. It can be seen that a few points of moving average can significantly decrease the RMSE value of f^k. Also, a 12 ∗ 12 bit look-up-table is a good choice for a trade-off between accuracy and resources.Figure 3
Performance analysis for frequency deduction: (a) RMSE ofP^ versus input frequency with SNR = 20 dB; (b) RMSE of P^ versus n with SNR = 20 dB and fk = 50 MHz; (c) RMSE of P^ versus sample SNR with fk = 50 MHz; (d) RMS error of f^k versus size of look-up-table with SNR = 20 dB, TLS L = 8, and fk = 50 MHz.
(a)
(b)
(c)
(d)Functional hardware simulations have been conducted to compare the performance of the frequency deduction algorithm and the FFT algorithm. The input frequency is an arbitrary value between 0 and 800 MHz. The input data are divided into 8 phases, and the operation clock isfclock = fs/8 = 200 MHz. The latency is simulated by Quartus II and Modelsim. The RMSE is simulated by Matlab using 12,000 trials. The results are shown in Table 2. We can conclude that the frequency deduction algorithm performs better performance in accuracy and response time.Table 2
Implementation results for FFT and frequency deduction algorithms.
SNR (dB)
L
Comparison between FFT and frequency deduction
FFT
Frequency deduction
N_FFT
RMSE (MHz)
Latency (ns)
X
RMSE (MHz)
Latency (ns)
10
8
1024
15.773
10790
4
10.239
105
2048
15.735
21080
8
8.477
125
16
1024
5.269
10795
4
6.831
110
2048
5.181
21085
8
6.237
130
20
8
1024
9.500
10790
4
2.262
105
2048
9.522
21080
8
1.670
125
16
1024
3.161
10795
4
1.219
110
2048
3.039
21085
8
0.967
130
30
8
1024
8.665
10790
4
0.676
105
2048
8.689
21080
8
0.487
125
16
1024
2.847
10795
4
0.339
110
2048
2.739
21085
8
0.251
130Table3 compares the computation complexity for frequency deduction and FFT with the parameters N_FFT = 1024 and L = 16. For the FFT method, its complexity can be found easily from [6]. The N_FFT-point FFT requires to calculate 2 × N_FFT × log2N_FFT real multiplications, 2 × N_FFT × log2N_FFT real additions, and N_FFT-point peak search. The complexity of frequency deduction is derived from (8) to (11) in this paper, which needs 3L + 4 multiplications, 3(L − 1) + 2 additions, 1 division, and some other operations. It is obvious that the proposed method requires far fewer resources than FFT.Table 3
Computation complexity comparison of the FFT and frequency deduction.
Method
MULT.
ADD.
DIV.
Search numbers
Square root
Look-up-table
FFT
20480
20480
0
1024
0
0
Frequency deduction
52
47
1
0
1
1
MULT. represents multiplication operation, ADD. represents addition operation, and DIV. represents division operation.


## 3.1.2. Performance Comparison between Frequency-Band Division and CRT

CRT described in [16, 21] applies only to the complex exponential signal, i.e., bk = 1. In the real signal case, the value of bk must be estimated in advance. In [29], a method to estimate bk by employing a triggering circuit to detect the waveform’s zero crossing point ‘O’ is proposed. The method needs extra hardware. For simplicity, we assume that the value of bk is estimated correctly for CRT. The configuration of simulation is K = 3, fs1=MΓ1=1500MHz, fs2=MΓ2=1600MHz, fs3=MΓ3=1700MHz, and fc is a random number and is uniformly distributed between 0 to 5,000 MHz. The folding frequency f^k is obtained by f^k=fk+Δfk, and fk is the accurate folding frequency. 12,000 trials are implemented for different SNR. The SNR is defined as SNR=20logM/σ which is the same as [21]. So, we can derive that(30)τ=3M⋅10−SNR/20,where M = 100 MHz is the greatest common divisor of fs1, fs2, and fs3.From the analysis of Section2.2, we can see that accurate judgment of mk is quite important. In Figure 4(a), we simulated the probability for judging mk correctly in different SNR. The result shows that the proposed algorithm can judge the folding integers mk more accurately than CRT. Since both robust CRT and closed-form CRT need to calculate mk one by one, the calculation process is complicated especially in low SNR, and it is difficult to guarantee the calculation accuracy of each mk. However, the proposed algorithm judges mk by comparing the errors of each sub-band, so it is less sensitive to noise.Figure 4
Performance comparison between robust CRT, closed-form CRT, and the proposed algorithm for the real signal: (a) probability of correction versus SNR; (b) RMSE versus SNR; (c) mean error versus SNR.
(a)
(b)
(c)Figures4(b) and 4(c) show the curves of RMSE value and mean error versus SNR. The theoretical RMSE value of CRT is given by Wang and Xia [21] as(31)fc,RMSE=σK.Since our algorithm is less sensitive to noise, it has better RMSE and mean error performance at low SNR. Additionally, we can see that when the SNR is larger than the bound which from the error boundτ=M/4, RMSE and mean error decrease linearly with SNR due to the robustness of CRT. In this case, the bound is about 16.8 dB, indicated by the black line in figures. However, the bound of our proposed algorithm is about 11 dB, which is lower than the CRT. So, it shows better robustness.Table4 compares the computation complexity between the frequency-band division algorithm and robust CRT and closed-form CRT. Through careful analysis (17), (18), and (28) in [16], we find that the robust CRT method needs to calculate 4Γ3(K − 1) +3 (Γ2 + Γ1) multiplications, 4Γ3 (K − 1) +3 (Γ2+Γ1) divisions, 4Γ3(K − 1) +4 (Γ2 + Γ1) additions, and (K − 1) Γ1+Γ2+⋯+ΓK searches. Similarly, the closed-form CRT in [21] needs to calculate 4K multiplications, 3K divisions, 2K + 1 additions, K rounding of integers, and 2K modular operations. Since the values of mkfsk can be registered in advance, there is no need for multiplication in our algorithm. The frequency-band division deblurring method needs only 6Q additions, 3Q operations for finding absolute values, and Q searches. Therefore, we can conclude that the deblurring algorithm proposed is faster and more effective, so it is more suitable for the application of DIFM.Table 4
Computation complexity comparison of the CRT and frequency-band division.
Method
MULT.
DIV.
Mod.
Round.
ADD.
ABS.
Search numbers
Robust CRT [16]
205
205
0
0
236
0
63
Closed-form CRT [21]
8
6
4
2
5
0
0
Frequency-band division
0
0
0
0
108
54
18
Mod. represents a modular operation, Round. represents rounding integer operation, and ABS. represents absolute value operation.


## 3.2. Experiment

Experiments have also been conducted on a digital radio frequency memory (DRFM) module (Figure5(a)) to verify the practicability of the DIFM technique. This module is originally designed for ECM and can be reused in our experiments with minor modifications. Signal source 1 is used to generate sampling clocks of different frequencies for the ADC. Signal source 2 is used to generate a monochromatic pulsed signal with a frequency from 0 to 5,000 MHz and a step size of 100 MHz. After being converted to the differential signal through Anaren B0430J50100A00 Balun, the input signal is injected into National Semiconductor ADC12D1800 ADC. The ADC can digitize the input signal with a sampling clock of up to 1,800 MHz and a data bit width of 12 bits. In the experiment, we set the sampling rates to be the same as in the simulation. The sampling data are transmitted to the Altera Stratix4SGX360 FPGA. Due to the limitation of experimental conditions, we only sample the input signal of one channel at a time and calculate folding frequency. Then, a deblurring operation is completed by combining the three channel results. All the algorithms were completed in FPGA. Quartus II 13.1 firmware fitting and timing analysis software shows a maximum clock speed of 200 MHz and FPGA logic resource usage of less than 5% for the 4SGX360 device. We have also used 12,000 trials for implementation. Figure 5(b) shows the experimental results of the mean deviation, maximum absolute error, and RMS error. It demonstrates that the DIFM technique proposed in this paper can measure the frequency from 0 to 5,000 MHz with high accuracy. Also, it can be observed that RMS errors of more than 90% frequency points are less than 0.5 MHz. The latency of this implementation is simulated as 38 FPGA clocks, which at the given clock rate yields(32)TDIFM=38∗fclk=190ns.Figure 5
Experiment platform and results. (a) DIFM hardware for the implementation of DIFM. (b) Experiment results of DIFM.
(a)
(b)


## 4. Conclusions

In this paper, the feasibility of the DIFM technique based on multichannel sub-Nyquist sampling is studied by combining the frequency deduction algorithm with the frequency-band division deblurring algorithm. Compared with the FFT algorithm, frequency deduction requires much fewer resources and can calculate folding frequency faster. Compared with CRT, the frequency-band division algorithm has much less computational complexity and is less sensitive to noise. Numerical simulations demonstrate that the performance of the frequency deduction algorithm can be significantly improved by increasing the sampling points, look-up-table size, and the moving average points.This DIFM receiver can directly digitalize the RF signal and measure the carrier frequency in a wide range with several low-rate ADCs. Owing to the limitation of the FPBW of the ADC, the IBW could not extend to the least common multiple offs1,fs1,…,fsk. If the frequency extends to more than the full power bandwidth, the amplitude decreases rapidly. Since the DIFM technique is amplitude-insensitive, we can measure the frequency from 0 to 2,800 MHz within 3 dB amplitude loss and 2,800 MHz to 5,000 MHz within acceptable amplitude loss in the environment of the ADC12D1800 device. Advancements in ADC will lead to the realization of greater IBW. As the excellent performance in response time, this DIFM technique can also be applied in linear frequency-modulated waveforms when the chirp rate is not too large.



---

*Source: 5089761-2020-03-24.xml*