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.

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 [414]. DFT algorithm is a simple algorithm widely used in frequency estimation [48]. In general, researchers focus on developing a variety of interpolation algorithms [46] 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 [1113] 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, 1820]. 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, 1820] 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 frequency fc 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 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.

2. DIFM Algorithm

2.1. Sub-Nyquist Sampling

When the input signal is sampled by an ADC with sample frequency fsk, 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 aswhere d is a nonnegative integer.

For each channel, if , 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 the kth channel which is defined aswhere A is the signal amplitude, φ is the initial phase, , and denotes a Gaussian distribution with a zero mean and variance . The signal-to-noise ratio (SNR) of the input signal is defined as .

If the sample period is , thenwhere N is the length of the signal.

When considering the idea condition without noise, i.e., , it is easy to prove thatwhere is a constant corresponding to fk one by one, .

Therefore, the calculation of fk 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 , 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

The error can be expressed as

Dividing the numerator and denominator of the right side of (7) by , we obtainwhere denote the numerator and denominator of the fraction, respectively.

Since , we get to know

It is obvious that when the expectation of and the variance of are constant, the variance of is proportional to the variance of and inversely proportional to the expectation of . So, we conclude that the RMSE value of is proportional to and inversely proportional to .

The estimation of P using (6) is based on an assumption that the denominator is nonzero. However, in practice, due to the noise and quantization level error, we cannot ensure is nonzero. Also, if 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 aswhere , .

The 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 least-squares (LS) method is adopted to estimate P. The LS method involves finding such thatwhere is an L vector representing the errors in Y and is the l2 norm given by

Then, the optimal fitting for P using the LS method iswhere 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. The error of the TLS method takes into account both the Y errors and the X errors. So, the TLS method involves finding such thatwhere denotes the Frobenius norm given by

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:

The sum of the squared errors of L points iswhere is defined as

Setting (′ represents differential operation), we obtain

We get the value of that makes the minimal value as

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.

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 . 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 , we have done moving average of S points at the expense of increasing S clock response time.

The frequency error is defined as , which is influenced by , , L, SNR, S, and the quantization level of look-up-table. For purposes of analysis, we consider 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 under fs/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 aswhere 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 aswhere LCM denotes least common multiple.

If K = 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 . 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 . The values of mk and bk can be calculated in advance aswhere denotes the flooring operation and rk is the remainder of fc modulo fsk.

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

The optimal sub-band is determined by calculating the error of each sub-band. We list all the values of and carry out subtraction of , , and . The sum of error between arbitrary two channels can be described aswhere 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. This 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. This processing procedure will continue until minimum is achieved.(3)Register the values of that give the least value of .

When the right sub-band is determined, fc is calculated using the average of three channels as

For all carrier frequencies, we can obtain

Therefore, the error of carrier frequency measurement is

Assume that is uniformly distributed between , is the error bound. The variance of is . Only if the sub-band is determined correctly, we can ensure that .

It should be noted that the proposed deblurring method is still valid in the case of multiple sinusoidal signals. However, the computation complexity is T3 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 is fs = 1600 MHz and pulse width = 20 us.

Figures 3(a)3(c) show the change rule of RMSE value of to the input frequency, time index, and SNR in different values of L. It can be seen that the RMSE value of is proportional to and inversely proportional to and SNR, which validates our theory. A few sampling points for both LS and TLS methods can significantly weaken the influence caused by , , 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 . 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. The input frequency is an arbitrary value between 0 and 800 MHz. The input data are divided into 8 phases, and the operation clock is fclock = 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 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]. 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.

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, , , , and fc is a random number and is uniformly distributed between 0 to 5,000 MHz. The folding frequency is obtained by , and fk is the accurate folding frequency. 12,000 trials are implemented for different SNR. The SNR is defined as which is the same as [21]. So, we can derive thatwhere M = 100 MHz is the greatest common divisor of fs1, fs2, and fs3.

From the analysis of Section 2.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.

Figures 4(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

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 , 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. 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 (Γ21) divisions, 4Γ3(K − 1) +4 (Γ2 + Γ1) additions, and (K − 1) 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 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.

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

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

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 61971179).