Enhanced Wideband Frequency Estimation via FFT: Leveraging Polynomial Interpolation and Array Indexing

Kiran Jayarama^*, Chien-In Henry Chen

Digital Receiver Laboratory, Department of Electrical Engineering, Wright State University, Dayton, OH, USA.

DOI: 10.4236/jcc.2024.121003                

Abstract

Accurate frequency estimation in a wideband digital receiver using the FFT algorithm encounters challenges, such as spectral leakage resulting from the FFT’s assumption of signal periodicity. High-resolution FFTs pose computational demands, and estimating non-integer multiples of frequency resolution proves exceptionally challenging. This paper introduces two novel methods for enhanced frequency precision: polynomial interpolation and array indexing, comparing their results with super-resolution and scalloping loss. Simulation results demonstrate the effectiveness of the proposed methods in contemporary radar systems, with array indexing providing the best frequency estimation despite utilizing maximum hardware resources. The paper demonstrates a trade-off between accurate frequency estimation and hardware resources when comparing polynomial interpolation and array indexing.

Keywords

Scalloping Loss, Goertzel’s Algorithm, Super-Resolution, Fast Fourier Transform (FFT), Decimation in Frequency (DIF), Decimation in Time (DIT), Spectral Leakage, Frequency Estimation

1. Introduction

The Fast Fourier Transform (FFT) is a fundamental signal processing tool for efficiently computing the frequency content of signals. As technology advances, the demand for more accurate frequency estimates becomes crucial to extract valuable information from FFT output data. Traditional techniques such as peak detection and interpolation have proven effective in many applications [1] [2] [3] [4] . However, in noisy environments and with closely spaced signals, implementing these methods for frequency determination becomes extremely challenging. While methods like Wavelet Transform [5] and Short-Time Fourier Transform (STFT) [6] show promise in identifying non-stationary signals, they are less effective in accurately estimating frequencies in real-world scenarios where signals fluctuate over time. Adaptive frequency estimation [7] dynamically adjusts parameters based on input signal characteristics. Several adaptive approaches have been developed over time to precisely estimate frequency [8] [9] [10] . While these techniques offer adaptability to address non-stationary signals with varying frequencies, they do have limitations. Many adaptive frequency estimation methods feature complex algorithms, leading to increased computational complexity, especially in hardware implementations. Iterative methods refine frequency estimates through successive iterations. After determining the main frequency, interpolation techniques are applied to fit a curve to the peak region in the FFT spectrum, enhancing the precision of the frequency estimate.

This paper introduces two innovative methods for enhancing frequency estimation precision, employing polynomial interpolation and array indexing. The organization of the paper is as follows: Section 2 discusses the wideband frequency estimation flow via FFT, followed by spectrum analysis of scalloping loss. Section 3 provides a brief description of the super-resolution method in section. A detailed explanation of the proposed methods is provided in Section 4. Section 5 presents the results, and the paper concludes in Section 6.

2. Proposed Wideband Frequency Estimation Flow via FFT

Digital receivers play a critical role in various applications such as warfare, radar, communication, and biomedical signal processing. Receiver-on-chip (ROC) designs are becoming exceedingly robust due to recent hardware advancements [11] [12] . In Figure 1, we illustrate the design flow using the Xilinx RFSoC [13] . The Tektronix 5200 serves as the signal generator, producing the necessary signals for the experiments. A 12-bit analog-to-digital (ADC) within the RFSoC samples the signal at 2.56 GHz. The digitized data is collected and aligned for pipelined parallel processing. To introduce pseudo-periodicity, a window function is selected [14] . The ADC returns 12-bit digitized samples, and the Multiple Input Selections (MIS) scheme is employed to maintain 12-bit precision for the input to the 512-point decimation in frequency (DIF) FFT built using a dynamic kernel [15] . In the dynamic kernel, the unit circle is scaled by the powers of two, and the real and imaginary components of the twiddle factors are mapped to closed integer points to minimize approximation errors. Adders and shifters replace the multipliers in the butterfly blocks to achieve a higher throughput rate and minimize hardware computational complexity. Decimation in frequency (DIF) is utilized to realize the FFT, as the approximated kernel function has an advantage over decimation in time (DIT) during hardware implementation. The peak magnitude $X\left[m\right]$ and the bin number m are recorded, along with $X\left[m\pm 1\right]$ and its respective bin numbers for frequency detection. This paper

Figure 1. Proposed wideband frequency estimation flow via FFT.

focuses on the accurate frequency estimation.

3. Frequency Estimation via FFT

The frequency spectrum of a signal represents the distribution of signal energy across different frequencies. This fundamental concept in signal processing is typically visualized through a graph or plot. The frequency spectrum offers valuable insights into the frequency content of a signal, facilitating the analysis and characterization of signals in diverse applications. The fast Fourier transform (FFT) is employed to compute the discrete frequency spectrum, yielding a set of discrete frequency components along with their corresponding amplitudes.

Frequency estimation is the process of determining the frequency content of the signal. The primary function of the FFT is to transform the time-domain signal to frequency-domain, presenting the signal in terms of its frequency components. Detecting the frequency becomes straightforward when the signal frequency aligns with an integer multiple of the frequency resolution. Challenges arise when the input signal falls between frequency bins. The subsequent sections outline four methods for determining frequency when the input signal is not precisely on a frequency bin.

3.1. Scalloping Loss

The FFT method effectively analyzes a finite-duration segment of the input signal. To mitigate spectral leakage, caused by the finite duration, the input signal is often multiplied by a window function. This process tapers the edges of the finite-duration segment to zero, reducing the effects of discontinuity. Scalloping loss, characterized by a reduction in the peak signal’s amplitude during the windowing process, is a notable outcome to identify signals not integer multiples of the frequency resolution.

In Figure 2, the plot depicts frequency on the x-axis and magnitude on the y-axis. A 512-point FFT analysis was conducted for frequencies ranging from 50 MHz to 60 MHz in increments of 500 KHz. The signals were sampled at 2.56 GHz, resulting in a frequency resolution of 5 MHz. Signals falling precisely on a frequency bin exhibit maximum magnitude, while those between bins have lower magnitudes. In this representation, 50 MHz, 55 MHz, 60 MHz, and 65 MHz fall precisely on frequency bins, resulting in higher magnitudes compared to 52.5 MHz, 57 MHz, and 62 MHz. The FFT spectrum in Figure 2 exhibits a series of peaks and valleys, resembling a scalloping pattern. Maximum scalloping loss is observed when the frequency lies precisely halfway between two bins.

Figure 3 shows the magnitude response of the frequencies between two frequencies centered on the bins. The mathematical representation of the magnitude of the frequency component is given in Equation (1). The peak magnitude $X\left[m\right]$ of a sinusoidal wave can lie anywhere between:

$\begin{array}{l}\frac{0.637\ast A\ast N}{2}\le X\left[m\right]\le \frac{A\ast N}{2}\hfill \end{array}$ (1)

where A is the amplitude of the time-domain signal and N is the length of the FFT [16] . To determine the signal’s time-domain amplitude A from Equation (1) by measuring the FFT spectral peak magnitude $X\left[m\right]$ is given by Equation (2).

Figure 2. Scalloping loss.

Figure 3. Normalized bin-to-bin magnitude response.

$\begin{array}{l}A=\frac{2\ast X\left[m\right]}{N}\hfill \end{array}$ (2)

The amplitude calculated using Equation (2) results in an error of 36.3%, which is equivalent to 3.9 dB. This level of error, while significant, may be acceptable for certain applications.

3.2. Reducing Scalloping Loss

One possible method to mitigate the frequency-dependent observed amplitude uncertainty inherent in FFT is to employ an N-sample flat-top window function on the original N time-domain samples. followed by executing a new FFT on the processed data. Flat-top window functions are designed to address the scalloping loss introduced during windowing. By using a flat-top window, the reduction in accuracy associated with scalloping can be minimized, given the importance of accurate amplitude representation. This involves multiplying the N-sample flat-top window function with the N-sampled time-domain data. It’s worth noting that multiplication in the time domain is equivalent to convolution in the frequency domain. To implement frequency domain convolution, we need to compute a single flat top windowed $X_{ft}\left(m\right)$ spectral sample from $X\left(m\right)$ spectral sample, using

$\begin{array}{c}{X}_{ft}\left(m\right)=X\left(m\right)-\frac{1.88494}{2}\ast \left[X\left(m-1\right)+X\left(m+1\right)\right]\\ +\frac{0.88494}{2}\ast \left[X\left(m-2\right)+X\left(m+2\right)\right]\end{array}$ (3)

where $X\left(m\right)$ is the peak magnitude of the FFT samples and m is the FFT’s peak magnitude bin value. From Equation (1), $\left|X_{ft}\left(m\right)\right|$ is the peak magnitude of the FFT sample. Figure 4 presents a plot comparing scalloping loss (green) and without scalloping loss (red) for frequencies between 450 MHz and 455 MHz. In this representation, scalloping loss has been minimized to below 0.2 dB.

This aids in accurately determining the sine wave’s time-domain amplitude accurately. The coefficients $g_0=1.0$ , $g_1=\frac{-1.88494}{2}=-0.94247$ and

$g_2=\frac{0.88494}{2}=0.44247$ eliminate the amplitude gain loss without changing

their scalloping loss performance. Additionally, these coefficients have the potential to decrease computational complexity through the utilization of binary representation and the substitution of multipliers with arithmetic right-shifters. When applying Equation (3) to rectangular windowed $X\left(m\right)$ FFT samples and compute the flat-top windowed maximum FFT spectral peak magnitude $M_{peak}=\left|X_{ft}\left(m\right)\right|$ . In this scenario, the estimated value of A from Equation (2) will incur an error of no more than ±0.0003 dB. This minimal error is visually represented by the nearly ideal solid red curve in Figure 4.

3.3. Super-Resolution for Frequency Estimation

When the signal frequency is an integer multiple of frequency resolution, the

Figure 4. Reduction in scalloping loss.

signal frequency is estimated at the peak bin m with peak magnitude $X\left[m\right]=\left(A\ast N\right)/2$ , where A is amplitude of the input sinusoid, and N is FFT length. Often, the signal frequency does not fall on the FFT frequency bin, thereby reducing the magnitude to $X\left[m\right]<\left(A\ast N\right)/2$ . The super-resolution method enables computing the true signal frequency with the help of the adjacent bins of the peak bin. The ratio of the peak magnitude $X\left[m\right]$ and the higher magnitude among the two adjacent bins is utilized to determine R.

$R=\frac{X\left[m\right]}{\max \left(X\left[m-1\right]\mathrm{,}X\left[m+1\right]\right)}$ (4)

$\alpha =\frac{1}{1+R}$ (5)

α is computed from Equation (5). If α is close to zero then the true frequency is close to the peak bin m whereas, if α is close to 1 corresponds to the true frequency close to $m\pm 1$ . The value of α is quantized between the peak bin and the bin with a higher magnitude among the two adjacent bins to determine the actual bin value, which is used to estimate the actual frequency of the time domain signal [17] .

3.4. Super-Resolution after Reducing Scalloping Loss for Frequency Estimation

Super-resolution is embraced for its efficiency in utilizing minimal hardware and conducting fewer computations to estimate the actual frequency. However, potential challenges arise due to significant noise levels and the risk of scalloping loss associated with this method. The super-resolution technique described in section 4 is applied after compensating for scalloping loss. Figure 5 depicts the value of α and R with and without the scalloping loss technique implemented, respectively. As the super-resolution algorithm relies on linear approximation for frequency estimation, incorporating it after the scalloping loss technique yields a slight improvement. The heightened peak magnitude assists in more accurately isolating the signal frequency.

4. Polynomial Interpolation and Array Indexing

4.1. Polynomial Interpolation

The super-resolution technique employs eight quantization points to estimate a new bin value when the actual frequency does not align with an integer multiple of the frequency resolution. Enhancing the number of quantization points is anticipated to yield improved frequency estimation. To achieve this, we fit a polynomial to the frequency spectrum affected by scalloping loss, effectively increasing the quantization points. Given the symmetry of the frequency spectrum around the midpoint between two bins, MATLAB’s linear interpolation using the fit function is chosen as the fitting type. Figure 6 illustrates the magnitude response with 50 quantization points.

Figure 5. Variations in α and R.

Figure 6. Curve fitting.

The ratio of the peak magnitude of the FFT spectrum to the peak magnitude after reducing scalloping loss yields a point m_b close to the quantization points derived from the interpolation. Utilizing the magnitudes of the adjacent bins from the peak bin, the quantized point is added or subtracted from the peak bin to generate a more accurate bin number corresponding to the input frequency. Suppose $X\left[m+1\right]>X\left[m-1\right]$ , m_b is added to the peak bin m or else subtracted from m to obtain the true bin value. The actual bin value is utilized to get the input frequency.

4.2. Array Indexing

At low Signal-to-Noise Ratio (SNR) values, the methods discussed in section 3 continue to face challenges in accurately estimating the frequency. To address the limitations posed by negative SNRs, a new approach named array indexing employing Goertzel’s algorithm [18] is presented. After detecting the bin value for the peak magnitude using the general FFT formula, the array indexing specifically designed for non-integer multiples of the frequency resolution.

Goertzel’s algorithm is a computational technique designed for efficiently calculating individual terms of the Discrete Fourier Transform (DFT). It proves beneficial for identifying the presence of specific frequency components in a signal. When a particular frequency is chosen for signal analysis, Goertzel’s algorithm iterates efficiently through each sample in the signal, computing coefficients based on the selected frequency and the signal’s length. The algorithm provides real and imaginary parts of the frequency components, close to the regular DFT or FFT. Magnitude and phase components can then be calculated from the obtained real/imaginary pair. An optimized version of the Goertzel algorithm, while faster and more straightforward than the basic Goertzel algorithm, lacks the provision of real and imaginary frequency components. Instead, it provides the relative magnitude squared [18] .

$Bi{n}_{value}=\left[00.80.60.40.200.20.40.60.80\right]$ (6)

The optimized Goertzel’s algorithm is particularly well-suited for our problem, given our focus on the magnitude spectrum. The peak magnitude for the frequencies at bin $m-1$ , i.e., 450 MHz, to bin $m+1$ , i.e., 460 MHz, is computed with increments of 1 MHz. For example, when an input signal of 453 MHz is sampled at 2.56 GHz using 512-point FFT, the peak magnitude $X\left[m\right]$ and its adjacent bin magnitudes $X\left[m\pm 1\right]$ are plotted. In Figure 7, the bin number is plotted on the x-axis, and its corresponding magnitude is plotted on the y-axis. The peak magnitude falls on bin number 91, and the frequency is estimated as 91*5 MHz = 455 MHz, but the input frequency is 453 MHz, and the estimated frequency is 454.65 MHz obtained using a super-resolution algorithm. Figure 7 illustrates the effectiness of Goertzel’s algorithm when the target frequency is swept from 450 MHz to 460 MHz with increments of 1 MHz. The magnitude values are stored in a $G_{mag}$ array. The max index value ( $G_{mag}$ ) is computed to determine a new bin value, which can estimate the frequency with higher

Figure 7. Goertzel Algorithm vs Ideal FFT.

accuracy. The $G_{mag}$ array has eleven values, and the $Bi{n}_{value}$ array comprises precomputed quantized bin values, as shown in Equation (6).

The first 0 indicates the bin $m-1$ and the second 0 indicates the bin m and the third 0 corresponds to the bin $m+1$ . From Figure 7, the peak magnitude value will be stored in the fourth position in $G_{mag}$ array. The value stored at the fourth position in $Bi{n}_{value}$ array is 0.4 which is utilized to determine the actual frequency. If the index value is less than 6, then the values from $Bi{n}_{value}$ array are subtracted from m else it is added to m. Here, the value 0.4 obtained from $Bi{n}_{value}$ array corresponding to the 4^th index value, we subtract the quantized bin value with m to obtain the true bin value. The new bin value is used to determine the frequency. For example, $m=91$ and $Bi{n}_{value}=0.4$ since $G_{mag}=4$ , subtract $Bi{n}_{value}$ from m to get the actual bin value $B=91-0.4=90.6$ . Multiplying B with the frequency resolution $\delta f$ presents us with the actual frequency (90.6*5 MHz = 453 MHz).

5. Experimental Results

The proposed algorithms have been implemented and tested using MATLAB. A 512-point Decimation in Frequency (DIF) FFT is employed, with a sampling frequency of 2.56 GHz, resulting in a frequency resolution of $\delta f=5\text{MHz}$ The input signal-to-noise ratio varies from −10 dB to 10 dB for each frequency. The input frequency is swept from 40 MHz to 1240 MHz in incremental steps of 3 MHz, with the addition of white Gaussian noise. The average frequency error and maximum frequency for all the methods are recorded. A total of (400 frequency values) × (4 methods) × (10 different SNRs) = 16,000 signal test cases are performed for the case study. The frequency error value for each signal (i.e., 400 frequency values) is calculated, and the average and maximum frequency errors are plotted in Figure 8 and Figure 9 respectively.

We observe a marginal improvement between super-resolution algorithms with scalloping loss and after reducing the scalloping loss. However, the polynomial interpolation method demonstrates substantial enhancement, particularly as we transition to positive SNR values. It’s noteworthy that this method exhibits insignificance when SNR values are more negative [19] . Among all the techniques, the array indexing method outperforms, showcasing significant improvement in both average and maximum frequency errors for both positive and negative SNRs.

The 512-point DIF FFT utilizing the Dynamic kernel incorporates adders and shifters to replace multipliers, enhancing the throughput and performance of the FFT. Table 1 outlines the hardware estimates for the four design methodologies. The polynomial interpolation method uses 50 registers to store pre-calculated quantized bin values and the hardware presented in Table 1. Array indexing utilizes the maximum hardware resources but yields the best results among the four methods. Goertzel’s algorithm can be implemented using an IIR filter but may be slower for designs requiring high-performance throughput. Nevertheless, this method proves highly efficient for design purposes and does not demand tight timing constraints.

Figure 8. Average frequency error.

Figure 9. Maximum frequency error.

Table 1. Hardware estimation.

6. Conclusion

Accurately estimating frequencies in a wideband digital receiver using the FFT algorithm poses challenges, notably spectral leakage resulting from the FFT’s assumption of signal periodicity. This paper introduces two innovative methods for precise frequency precision: polynomial interpolation and array indexing, comparing their performance against super-resolution and scalloping loss. The experimental configuration involves a signal-to-noise ratio ranging from −8 dB to 10 dB and a signal frequency sweeping from 40 MHz to 1240 MHz, incorporating white Gaussian noise. Maximum frequency precision errors for four methods’ super resolution with scalloping loss, super resolution after reducing scalloping loss, polynomial interpolation, and array indexing’ are reported. Super resolution with scalloping loss exhibits a maximum frequency error ranging from 5.5 to 3.7 MHz, while super resolution after reducing scalloping loss shows a maximum frequency error ranging from 5.5 to 3.5 MHz. Polynomial interpolation has a maximum frequency error ranging from 5.5 to 1.4 MHz, and array indexing ranges from 0.8 to −0.2 MHz. Among the four hardware implementations, array indexing, while utilizing maximum hardware resources, yields the most accurate frequency estimation. The paper goes deeper into exploring the balance between getting accurate frequency measurements and using computer resources wisely. It takes a closer look at finding the right mix between two methods’ polynomial interpolation and array indexing. The goal is to understand how to achieve precise frequency results while making the best use of the hardware resources. This analysis provides insights that can be easily applied in real-world scenarios, helping design systems that are both accurate and efficient.

Conflicts of Interest

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

References

[1]	Rabiner, L.R. and Gold, B. (1975) Theory and Application of Digital Signal Processing. Prentice-Hall, Englewood Cliffs.
[2]	Oppenheim, A.V. (1978) Applications of Digital Signal Processing. Prentice Hall, Hoboken, New Jersey.
[3]	Roberts, R.A. and Mullis, C.T. (1987) Digital Signal Processing. Addison-Wesley Longman Publishing Co., Inc., Boston.
[4]	Smith, S.W., et al. (1997) The Scientist and Engineer’s Guide to Digital Signal Processing. http://www.dspguide.com/pdfbook.htm
[5]	Todorovska, M.I. (2001) Estimation of Instantaneous Frequency of Signals Using the Continuous Wavelet Transform. University of Southern California, Department of Civil Engineering, Los Angeles.
[6]	Griffin, D. and Lim, J. (1984) Signal Estimation from Modified Short-Time Fourier Transform. IEEE Transactions on Acoustics, Speech, and Signal Processing, 32, 236-243. https://doi.org/10.1109/TASSP.1984.1164317
[7]	Etter, D. and Hush, D. (1987) A New Technique for Adaptive Frequency Estimation and Tracking. IEEE Transactions on Acoustics, Speech, and Signal Processing, 35, 561-564. https://doi.org/10.1109/TASSP.1987.1165146
[8]	Abdoush, Y., Pojani, G. and Corazza, G.E. (2019) Adaptive Instantaneous Frequency Estimation of Multicomponent Signals Based on Linear Time-Frequency Transforms. IEEE Transactions on Signal Processing, 67, 3100-3112. https://doi.org/10.1109/TSP.2019.2912132
[9]	Shen, T.-A., Li, H.-N., Zhang, Q.-X. and Li, M. (2017) A Novel Adaptive Frequency Estimation Algorithm Based on Interpolation FFT and Improved Adaptive Notch Filter. Measurement Science Review, 17, 48-52. https://doi.org/10.1515/msr-2017-0006
[10]	Pin, G., Wang, Y., Chen, B. and Parisini, T. (2019) Identification of Multi-Sinusoidal signals with Direct Frequency Estimation: An Adaptive Observer Approach. Automatica, 99, 338-345. https://doi.org/10.1016/j.automatica.2018.10.026
[11]	Jha, S.S. and Darak, S.J. (2022) Design & Performance Analysis of Wireless Radar Receiver for Joint Radar Communication Systems on RFSoC. https://repository.iiitd.edu.in/jspui/bitstream/handle/123456789/1120/Shragvi%20Sidharth%20Jha.pdf?sequence=1&isAllowed=y
[12]	Rudersdorfer, R., Graf, U. and Zahnd, H. (2013) Radio Receiver Technology: Principles, Architectures and Applications. John Wiley & Sons, Hoboken, New Jersey. https://doi.org/10.1002/9781118647882
[13]	XILINX (2023) Zynq Ultrascale+ Device Technical Reference Manual UG1085. https://www.origin.xilinx.com/support/documents/user_guides/ug1085-zynq-ultrascale-trm.pdf
[14]	Tsui, J.B. (2010) Special Design Topics in Digital Wideband Receivers. Artech House, Washington, DC.
[15]	Lee, Y.-H.G. (2009) Dynamic Kernel Function for High-Speed Real-Time Fast Fourier Transform Processors. Ph.D. Dissertation, Wright State University, Dayton, OH.
[16]	Lyons, R.G. (1997) Understanding Digital Signal Processing, 3/E. Pearson Education India.
[17]	Chen, C.-I., George, K., McCormick, W., Tsui, J.B., Hary, S.L. and Graves, K.M. (2005) Design and Performance Evaluation of a 2.5-Gsps Digital Receiver. IEEE Transactions on Instrumentation and Measurement, 54, 1089-1099. https://doi.org/10.1109/TIM.2005.847206
[18]	Gazi, O. (2018) Understanding Digital Signal Processing. Springer, Singapore. https://doi.org/10.1007/978-981-10-4962-0
[19]	Lyons, R. (2011) Reducing FFT Scalloping Loss Errors without Multiplication. IEEE Signal Processing Magazine, 28, 112-116. https://doi.org/10.1109/MSP.2010.939845