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

Multichannel sub-Nyquist sampling is an eﬃcient technique to break through the limitation of the Nyquist sampling theorem for the wideband digital instantaneous frequency measurement (DIFM) receiver. The signiﬁcant 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 eﬃcient 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 conﬁrm the eﬀectiveness of our algorithms. The simulation results show that our algorithms have better performance than FFTand CRT in accuracy and response time. Board-level measurement results verify that the proposed DIFM technique can capture the frequency from 0 to 5,000MHz with a maximum of 0.8MHz root-mean-squared error in less than 200ns.


Introduction
Instantaneous frequency measurement (IFM) receiver plays an important role in military and civilian applications, such as electronic countermeasures (ECM), electronic countercountermeasures, 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]. e 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 general, researchers focus on developing a variety of interpolation algorithms [4][5][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.
is will increase hardware complexity and cost. Monobit ADC sampling [9,10] and photonic technology [11][12][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. ese technologies cannot deliver superior performance in terms of IBW, response time, precision, and hardware complexity at the same time.
e 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 f k (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][19][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][19][20] ignores that the remainders may be folding frequency f k or image frequency f skf k (f sk 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. e 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 frequency f c and f k 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.
is deblurring algorithm does not require that all the sampling rates are coprime as CRT. So, it is convenient for engineering. e 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. e simulation results show that compared with FFT, the frequency deduction algorithm can measure the frequency faster with the same accuracy. Also, the proposed frequencyband 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.
e remainder of this paper is organized as follows. In Section 2, 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.

Sub-Nyquist Sampling.
When the input signal is sampled by an ADC with sample frequency f sk , 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.
is can be mathematically expressed as where d is a nonnegative integer. For each channel, if f c ∈ (df sk , (d + (1/2)f sk )], f k is the remainder frequency as described in CRT; otherwise, f k is the image frequency. Obviously, f k corresponds to multiple f c , 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.

Frequency Deduction Based on the TLS Estimation
Algorithm. Assume a monochromatic sinusoidal signal of the kth channel which is defined as where A is the signal amplitude, φ is the initial phase, ω(t) ∼ N(0, σ 2 ), and N(0, σ 2 ) denotes a Gaussian distribution with a zero mean and variance σ 2 . e signal-to-noise ratio (SNR) of the input signal is defined as SNR � (A 2 /2σ 2 ).
If the sample period is T sk (T sk � 1/f sk ), then 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 where P � cos(2πf k T sk ) is a constant corresponding to f k one by one, |P| ≤ 1. erefore, the calculation of f k 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: 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 e error ΔP can be expressed as Dividing the numerator and denominator of the right side of (7) by 2σ, we obtain where Η � (ω(n + 1)/2σ) + (ω(n − 1)/2σ) − 2Pω(n)/2σ and Τ � A/σ cos(2πf k nT sk ) + ω(n)/σ denote the numerator and denominator of the fraction, respectively.
Since ω(n) ∼ N(0, σ 2 ), we get to know 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 | ����� 2SNR √ cos(2πf k nT sk )|. e estimation of P using (6) is based on an assumption that the denominator s k (n) is nonzero. However, in practice, due to the noise and quantization level error, we cannot ensure s k (n) is nonzero. Also, if s k (n) is close to zero, the estimation result will be significantly affected by noise.
To solve the problem and estimate P correctly, we construct a linear mathematical model with L sampling points as where e problem of estimating the value of P is converted to fit an optimal straight line that passes through the origin and minimizes the total errors of L points. Usually, the leastsquares (LS) method is adopted to estimate P. e LS method involves finding P LS such that

Mathematical Problems in Engineering 3
where ΔY is an L vector representing the errors in Y and ‖ · ‖ 2 is the l 2 norm given by en, the optimal fitting for P using the LS method is where U and V are defined as However, since both Y 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. e error of the TLS method takes into account both the Y errors and the X errors. So, the TLS method involves finding P where ‖ · ‖ F denotes the Frobenius norm given by e 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: e sum of the squared errors of L points is where W is defined as Setting f ′ (P TLS L ) � 0 (′ represents differential operation), we obtain We get the value of P TLS L that makes f(P TLS L ) the minimal value as .
Obviously, U and V can hardly be zero especially for large L. rough 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. e 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.
In addition, we note that estimating P 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: Coordinated rotation digital computer (CORDIC) algorithm or look-up-table can be used to calculate cos − 1 (P). 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.
e frequency error is defined as Δf k � f k − f k , which is influenced by |P|, | ����� 2SNR √ cos(2πf k nT sk )|, L, SNR, S, and the quantization level of look-up-table. For purposes of analysis, we consider Δf k being uniformly distributed over a certain range. e 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.

Frequency-Band Division
Deblurring. Based on the Nyquist sampling theorem, only the frequency under f s /2 can be exactly recovered without aliasing. Ambiguity occurs when the frequency is extended to more than f s /2. To reconstruct the unambiguous frequency, we employ a multichannel sampling technique. (1) can also be expressed as where m k is an unknown integer, b k � 1 when f k is the remainder frequency, b k � -1 when f k is the image frequency, and M k is the integer part of f up /f sk , where f up is the upper bound of measurable frequency. Usually, M k 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. e maximum unambiguous frequency determined by K channels is derived and can be expressed as where LCM denotes least common multiple. If K � 2, it is easy to find that f max 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, f max � {LCM (f s1 , f s2 ) + f s3 }/2 (f s1 < f s2 < f s3 ). Proper setting of sampling rates of three channels can guarantee that the value of f max is larger than the FPBW of ADC. erefore, all the frequencies that fall into the FPBW can be uniquely reconstructed.
From (23), the reconstruction of f c is transformed into the calculation of f k , m k , and b k . e calculation of f k has been introduced in Section 2.2. In this section, we study the fast algorithm for determining m k and b k 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 m 1 , m 2 , m 3 , b 1 , b 2 , b 3 . e values of m k and b k can be calculated in advance as where ⌊ · ⌋ denotes the flooring operation and r k is the remainder of f c modulo f sk . Table 1 presents an example for f s1 � 1500, f s2 � 1600, and f s3 � 1700. Usually, Q will not be too large since M k is small. So, we do not need to search for too many numbers. For example, if f c ∈ [0, 2f s3 ), Q � 12. So, the frequency range can be divided into 12 sub-bands. We just have to judge the subband which the carrier frequency falls into. e optimal sub-band is determined by calculating the error of each sub-band. We list all the values of m 1 , m 2 , m 3 , b 1 , b 2 , b 3 and carry out subtraction of m 1 f s1 + b 1 f 1 , m 2 f s2 + b 2 f 2 , and m 3 f s3 + b 3 f 3 . e sum of error between arbitrary two channels can be described as where f i error is the total error of three channels and i is the sub-band index.
All Q data are compared by the pipeline technique to get the minimum value.
is process can be expressed as follows: (1) Divide Q data into Q/2 groups and do a pairwise comparison among all Q data. (2) Choose smaller Q/2 from the results of (1), and continue the next comparison. is processing Mathematical Problems in Engineering procedure will continue until minimum f l error is achieved.
(3) Register the values of m l 1 , m l 2 , m l 3 , b l 1 , b l 2 , b l 3 that give the least value of f l error . When the right sub-band is determined, f c is calculated using the average of three channels as (27) For all carrier frequencies, we can obtain erefore, the error of carrier frequency measurement is (29) Assume that Δf k is uniformly distributed between [− τ, τ], τ is the error bound. e variance of Δf k is σ 2 � τ 2 /3. Only if the sub-band is determined correctly, we can ensure that |Δf c | ≤ τ.
It should be noted that the proposed deblurring method is still valid in the case of multiple sinusoidal signals. However, the computation complexity is T 3 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. erefore, 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.

Simulation for Frequency Deduction.
e performance of the frequency deduction algorithm is verified by a single channel frequency simulation. e configuration of simulation is f s � 1600 MHz and pulse width � 20 us. Figures 3(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 |cos(2πf k nT s )| and SNR, which validates our theory. A few sampling points for both LS and TLS methods can significantly weaken the influence caused by |P|, |cos(2πf k nT s )|, 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.
Functional hardware simulations have been conducted to compare the performance of the frequency deduction algorithm and the FFT algorithm. e input frequency is an arbitrary value between 0 and 800 MHz. e input data are divided into 8 phases, and the operation clock is f clock � f s / 8 � 200 MHz. e latency is simulated by Quartus II and Modelsim. e RMSE is simulated by Matlab using 12,000 trials. e results are shown in Table 2. We can conclude that the frequency deduction algorithm performs better performance in accuracy and response time. Table 3 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]. e N_FFT-point FFT requires to calculate 2 × N_FFT × log2N_FFT real multiplications, 2 × N_FFT × log2N_FFT real additions, and N_FFTpoint peak search. e 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.

Performance Comparison between Frequency-Band
Division and CRT. CRT described in [16,21] applies only to the complex exponential signal, i.e., b k � 1. In the real signal case, the value of b k must be estimated in advance. In [29], a method to estimate b k by employing a triggering circuit to detect the waveform's zero crossing point 'O' is proposed.  The size of look-up-table (bits)  e method needs extra hardware. For simplicity, we assume that the value of b k is estimated correctly for CRT. e configuration of simulation is K � 3, f s1 � MΓ 1 � 1500 MHz, f s2 � MΓ 2 � 1600 MHz, f s3 � MΓ 3 � 1700 MHz, and f c is a random number and is uniformly distributed between 0 to 5,000 MHz. e folding frequency f k is obtained by f k � f k + Δf k , and f k is the accurate folding frequency. 12,000 trials are implemented for different SNR.
e SNR is defined as SNR � 20 log(M/σ) which is the same as [21]. So, we can derive that where M � 100 MHz is the greatest common divisor of f s1 , f s2, and f s3 .  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.  From the analysis of Section 2.2, we can see that accurate judgment of m k is quite important. In Figure 4(a), we simulated the probability for judging m k correctly in different SNR. e result shows that the proposed algorithm can judge the folding integers m k more accurately than CRT. Since both robust CRT and closed-form CRT need to calculate m k one by one, the calculation process is complicated especially in low SNR, and it is difficult to guarantee the calculation accuracy of each m k . However, the proposed algorithm judges m k by comparing the errors of each subband, so it is less sensitive to noise. Figures 4(b) and 4(c) show the curves of RMSE value and mean error versus SNR. e theoretical RMSE value of CRT is given by Wang and Xia [21] as 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. Table 4 compares the computation complexity between the frequency-band division algorithm and robust CRT and closed-form CRT. rough 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 closedform 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 m k f sk can be registered in advance, there is no need for multiplication in our algorithm. e frequency-band division deblurring method needs only 6Q additions, 3Q operations for finding absolute values, and Q searches. erefore, we can conclude that the deblurring algorithm proposed is faster and more effective, so it is more suitable for the application of DIFM.

Experiment.
Experiments have also been conducted on a digital radio frequency memory (DRFM) module ( Figure 5(a)) to verify the practicability of the DIFM technique. is 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. e 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. e sampling data are transmitted to the Altera Strat-ix4SGX360 FPGA. Due to the limitation of experimental conditions, we only sample the input signal of one channel at a time and calculate folding frequency. en, 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.
e latency of this implementation is simulated as 38 FPGA clocks, which at the given clock rate yields T DIFM � 38 * f clk � 190 ns. (32)

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. is 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 of f s1 , f s1 , . . . , f sk . 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.

Data Availability
All the data used for the numerical simulations and comparison purpose have been reported in the tables included and visualized in the graphical illustrations.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.