FFT-Based Simultaneous Calculations of Very Long Signal Multi-Resolution Spectra for Ultra-Wideband Digital Radio Frequency Receiver and Other Digital Sensor Applications

The discrete Fourier transform (DFT) is the most commonly used signal processing method in modern digital sensor design for signal study and analysis. It is often implemented in hardware, such as a field programmable gate array (FPGA), using the fast Fourier transform (FFT) algorithm. The frequency resolution (i.e., frequency bin size) is determined by the number of time samples used in the DFT, when the digital sensor’s bandwidth is fixed. One can vary the sensitivity of a radio frequency receiver by changing the number of time samples used in the DFT. As the number of samples increases, the frequency bin width decreases, and the digital receiver sensitivity increases. In some applications, it is useful to compute an ensemble of FFT lengths; e.g., 2P−j for j=0, 1, 2, …, J, where j is defined as the spectrum level with frequency resolution 2j·Δf. Here Δf is the frequency resolution at j=0. However, calculating all of these spectra one by one using the conventional FFT method would be prohibitively time-consuming, even on a modern FPGA. This is especially true for large values of P; e.g., P≥20. The goal of this communication is to introduce a new method that can produce multi-resolution spectrum lines corresponding to sample lengths 2P−j for all J+1 levels, concurrently, while one long 2P-length FFT is being calculated. That is, the lower resolution spectra are generated naturally as by-products during the computation of the 2P-length FFT, so there is no need to perform additional calculations in order to obtain them.


Introduction
Multi-resolution frequency spectrum analysis is widely used in many areas to process signals from different sources, such as signals from mechanical structures, acoustic sensors, human body signal sensors, and digital RF/microwave receivers.A few examples are given in [1][2][3][4].Multi-resolution analysis (MRA) is a powerful tool used in image and signal processing.In wavelet analysis [5], a mother wavelet function [6] is scaled in frequency and shifted in time to analyze images or signals at different resolutions [6,7].Recently, Ref. [8] gives a comprehensive review of the developments, challenges, and opportunities of wavelet analysis, and also covers the evolution of wavelet neural networks.In parallel with the developments in wavelet-based MRA, Fourier transform (FT)-based approaches are also being developed.This framework includes the fractional FT [9][10][11], and the relatively recently developed quadratic-phase FT [12,13], which offers a unified approach to processing transient and non-transient signals and is very useful for analyzing signals with time-varying properties.Although many FT-based MRA methods have been developed, in modern digital receiver development the conventional FT is still widely used.In general, for a fixed digitizer sampling rate, calculating DFTs having different numbers of samples results in spectra having different frequency resolutions.In the FT-based MRA, the short-time Fourier transform (STFT) is often used to calculate signal time-frequency Sensors 2024, 24, 1207 2 of 10 distributions.There are a number STFT-based methods, such as the multi-resolution STFT [14], adaptive STFT [15,16], and STFT frequency-domain (STFT-FD) [17] methods.The idea of the multi-resolution STFT is based on using different window sizes to achieve different frequency resolutions.It uses bandpass filters to divide the signal into frequency bands.Therefore, multi-resolution STFT relies on the use of filters to obtain the multiresolution spectra.The adaptive STFT adjusts the window size for each time instant based on the local signal characteristics.The STFT-FD computes a transform using the basic idea of the STFT; however, it uses a fixed window size in the frequency domain.This approach is similar to that of the multi-resolution STFT and wavelet methods, but with a different methodology.The key of the STFT-FD is to define the window size as the number of cycles of the frequency component being considered rather than defining it in time.The STFT-FD needs neither filtering nor the evaluation of local signal characteristics as in the multi-resolution STFT and the adaptive STFT, respectively.
It is well known that any STFT-based time-to-frequency transform has a limitation.This is the inherent trade-off between time and frequency resolution, i.e., the short-time window has good time resolution, but poor frequency resolution, and vice versa.In [18], the piecewise-constant window blocking fast Fourier transform (PCW-BFFT) method was introduced for ultra-wideband digital receivers in order to achieve a fine frequency resolution and a high receiver sensitivity.By applying a number of window functions with different widths in the time domain, the PCW-BFFT method can simultaneously produce multi-time-resolution spectrograms with a very fine frequency resolution for a very long-time sequence (e.g., N = 2 20 ); however, this method cannot produce multifrequency-resolution spectra concurrently.
The DFT-based time-to-frequency transform is the most widely used signal processing method in ultra-wideband digital sensors and is commonly implemented using the fast Fourier transform (FFT) on field programmable gate arrays (FPGAs) in real-time applications.Ref. [19] presents several methods using the FFT to detect and analyze intercepted signals with multi-resolution frequency spectral lines (MR-FSL).However, up to now, to obtain the MR-FSL of a signal, multiple FFT calculations with different signal lengths (e.g., N = 2 P−j where j = 0, 1, 2, . . ., J) must be performed one by one.The disadvantages of the traditional FFT-based methods are: (1) for large signal length (e.g., P = 20 and j = 0, 1, 2, 3), four different long-length FFTs have to be calculated; (2) one has to wait for a large number of time samples to be generated by the analog-to-digital converter (ADC) in the receiver before the FFTs can be performed; (3) many FPGA resources are needed to store those time samples before the FFTs can be performed; and (4) more importantly, the computational effort spent on the calculation of the longest FFT cannot be reused in the calculations of any of the shorter length FFTs.
To mitigate the problem of requiring a large amount of hardware resources to compute very long length FFTs and effectively calculate the frequency spectrum of long-sequence signals in an embedded device in real time, the blocking FFT (BFFT) method was first introduced in [20].It also demonstrated the implementation of the method on Xilinx's ZYNQ-7000 device.In order to perform an N-point frequency bin calculation, the BFFT conceptually divides N time samples into K consecutive timeslot blocks consisting of M samples each.As each block of M (= N/K) samples is generated by the ADC, the BFFT method performs K independent M-point FFTs, where suitable phase rotation factors are applied to the M input samples of each of the K FFTs.Each of these FFTs can be computed in parallel (e.g., on an FPGA or GPU).The resulting K × M = N frequency data are accumulated with the results from previous timeslots to form an N-bin spectrum line.Each spectrum line can be viewed as a group of K equal-length frequency blocks where the k th frequency block has frequency bins f kq , where k = 0, 1, 2, . . ., K − 1, and q = k, M + k, . . ., (K − 1)M + k.The final N-bin spectrum line is calculated as the accumulation of the frequency results from processing all K timeslots.The merit of the BFFT method is that instead of having to first store a large number of N time samples before being able to compute a single, long N-point FFT, which may be a challenge even on today's high-end FPGAs, the BFFT performs K many M-point FFTs in parallel, as soon as M samples are available from the current timeslot.This tremendously reduces the hardware resource requirements, produces useful intermediate results with low latency, and also enables the development of other useful methods for ultra-wideband digital receiver designs, such as the multi-time-resolution piecewise-constant window blocking Fourier transform (PCW-BFFT) method [18], the accumulatively increasing receiver sensitivity (AIRS) method [21], and the multi-frequency-resolution spectrum line calculation method introduced in this paper.
Taking advantage of the aforementioned BFFT process, this paper for the first time introduces an FFT-based MR-FSL calculation method, called the multi-frequency-resolution BFFT (MFR-BFFT).It can calculate the spectrum line having N 0 = 2 P frequency bins, while simultaneously producing all the spectrum lines having N j = 2 P−j frequency bins, where j = 0, 1, 2, . . ., J, without needing to perform extra FFT calculations on the M samples from the current timeslot.That is, the computational efforts to obtain the 2 P -bin spectrum line are reused to obtain the spectrum lines that have 2 P−1 , 2 P−2 , 2 P−3 , . . ., 2 P−J frequency bins at other frequency resolution levels.With a fixed ADC sampling rate (i.e., f s = 1 t s where t s is the time interval between samples), the MFR-BFFT uses 2 P−j time samples (where j = 0, 1, 2, . . ., J) to compute spectrum lines up to J = log 2 (K 0 ) levels with 2 j different spectrum lines at Level-j.Here, K 0 is the number of timeslots in Level-0.Since the ADC rate is fixed, the frequency resolutions at different spectrum levels are 2 j •∆ f , where ∆ f = f s /2 P is the Level-0 frequency resolution.The signal duration used to get different level spectrum lines is N j × t s .More details will be presented in the next section.
Compared with other aforementioned existing MRA methods, although the MFR-BFFT method does not have the flexibility to obtain different levels of resolutions as wavelet-based MRA methods do, it can tremendously reduce the computational effort on hardware to produce spectra with a fixed set of frequency resolutions.Compared with other FT-based MRA methods, the MFR-BFFT has the advantage of being able to process very long-length signal samples, and to produce a set of different frequency resolution spectra simultaneously.Note that the very long-length signal samples do not necessarily mean a long signal in the time domain.It can be a long signal sampled at a low ADC sampling rate, or a very short signal sampled at high sampling rate.For example, if an ADC samples at 2.5 GSample/s, within about 0.42 ms, it can produce 2 20 time samples.
The rest of the paper is organized as follows: The derivations and discussions of the MFR-BFFT methods are presented in the next section.Section 3 shows the MFR-BFFT results using data obtained from an ultra-wideband digital radio receiver and compares the results with those obtained from the "fft" function of Mathworks MATLAB (R2023b) to prove the correctness of the MFR-BFFT implementation in MATLAB.The conclusion can be found in the last section.The hardware implementation including an expansion of the comparative analysis between MFR-BFFT and conventional FFT, such as computational complexity, power consumption, and execution time, etc., will be reported in another full paper.

FFT-Based Multi-Frequency-Resolution Spectrum Line Calculation
Based on the discrete Fourier transform (DFT), the frequency spectrum line at Level-0 can be expressed as: where v 0 = 0, 1, 2, . . ., N 0 − 1 are the frequency bin centers on the Level-0 spectrum line, , and x is a complex-valued signal.X 01 means the first spectrum line at Level-0, although at this level there is only one spectrum line.Using the BFFT, N 0 time samples are divided into K 0 equal timeslots and each timeslot has M samples; therefore, (1) can be written as: By substituting u = n − kM, then for the k th timeslot of M samples, the portion inside the curly braces can be rewritten as: We now define: which is the frequency information obtained from the k th timeslot of M samples for the Level-0 spectrum line.Using ( 4), ( 2) can be rewritten as: In order for the M samples in the k th timeslot to contribute to the overall N 0 bins, the BFFT also divides N 0 frequency bins into K 0 equal-length frequency blocks, where each frequency block also has M frequency bins.By substituting v 0 = m 0 + rK 0 , where m 0 = 0, 1, 2, . . ., K 0 − 1 to define the total K 0 frequency blocks at Level-0, and r = 0, 1, 2, . . ., M − 1 to specify the M frequency bins in each frequency block, then (4) becomes: where MK 0 = N 0 and W kMK 0 r N 0 = 1 are used.Hence, using the BFFT to calculate (1), we have: It can be seen that to contribute to all the frequency bins in X 01 using M samples in the k th timeslot, the samples need to be multiplied by phase rotation factors W m 0 (u+kM) N 0 , which are determined by the frequency block (m 0 ) and time index (n = u + kM) in the overall time sequence, before performing the M-point FFT (the operations specified in curly braces) indicated in (7).Since m 0 = 0, 1, 2, . . ., K 0 − 1, a total of K 0 many M-point FFTs are computed at Level-0.Now, let us consider the Level-1 spectrum line calculations with the same time sequence x(n) where (n = 0, 1, 2, . . ., N 0 − 1).Using the DFT, the spectrum lines can be computed as: where i = 1, 2 indicates that there are the two spectrum lines in Level-1.The first is calculated from samples 0 to N 1 − 1, and the second is calculated from samples N 1 to . Note that (1) these two spectrum lines have the same frequency bins v 1 = 0, 1, 2, . . ., N 1 − 1; (2) they only use half of the available samples in x(n) compared with what was used in Level-0; and (3) the frequency bin size at Level-1 is double that of Level-0, since the ADC sampling rate is fixed.Again, if we use the BFFT to calculate the Level-1 spectrum lines with the same number of samples (M) in a timeslot that is used at Level-0, we have: where the first and second K 1 (= K 0 /2) groups of timeslots are used to calculate X 11 (v 1 ) and X 12 (v 1 ), respectively.By substituting u = n − kM, then for the k th timeslot's M samples, we have: In order for the M samples in k th timeslot to contribute to the overall N 1 bins, the BFFT also divides the N 1 frequency bins into K 1 blocks, each of which contains M frequency bins.By substituting v 1 = m 1 + rK 1 , where m 1 = 0, 1, 2, . . ., K 1 − 1 to define K 1 frequency blocks in each spectrum line, and r = 0, 1, 2, . . ., M − 1 (i.e., there are M frequency bins in each frequency block), then from (10) we have: where MK 1 = N 1 and W kMK 1 r N 1 = 1 are used.Then, using the BFFT, the spectrum lines at Level-1 can be expressed as: In general, using the BFFT, the spectrum lines at Level-j can be expressed as: where j = 0, 1, 2, . . ., log 2 (K 0 ) defines the spectrum line levels; i = 1, 2, 3, . . ., 2 j denotes the spectrum lines in Level-j; and the i th spectrum line is obtained using timeslots from (i − 1)K j to iK j − 1, where K j specifies the number of timeslots in Level-j.
From the above discussions we can observe that: 1. with fixed N 0 samples used in Level-0, the maximum number of levels is limited by the number of samples (M) in a timeslot to J = log 2 (K 0 ) = log 2 (N 0 /M); 2.
although the M samples of a given timeslot can be used in all levels, the differences are: • the number of M-point FFTs is determined by K j = K 0 /2 j , and • the number of frequency bins is determined by N j = N 0 /2 j ; 3. more importantly, at different levels, the same M samples in the k th timeslot are multiplied with phase rotation factors W m j (u+kM) N j in (13), which can be rewritten as , given that N j = N 0 /2 j .Therefore, (13) can be expressed as: Sensors 2024, 24, 1207 6 of 10 Comparing the portion inside the curly braces in (14) with the portion inside the curly braces in (7), we find that after we calculate the DFT (the portion inside the curly braces of ( 7)) using M samples in a timeslot, the DFT results can not only be used for frequency bins at Level-0, but some of them can also be used for other-level spectrum lines, as long as: 1.
these other-level spectrum lines also use M samples from the current timeslot for the calculation; and 2.
the frequency block number m j satisfies the relation given in (15), which means these other-level spectrum lines not only share the timeslot, but also use the same phase rotations factors.
where again m 0 = 0, 1, 2, . . ., K 0 − 1, and m j = 0, 1, 2, . . ., K j − 1, where j = 0, 1, 2, . . ., log 2 (K 0 ) and K j = K 0 /2 j .The above discussion infers that the other-level spectrum lines can be obtained at the same time that the Level-0 spectrum line is being calculated using the BFFT.Since the ADC sampling rate is fixed, and the number of time samples is reduced by half at each progressively higher spectrum level, while the frequency resolution is also reduced by half.The number of spectrum lines is doubled when going from Level-j to Level-(j + 1) if all the time samples used in Level-0 are used in the computation of the spectrum lines at other levels.
To make the above discussion clear, let us see an example.Here, we use N 0 = 32 time samples with M = 4 samples in each timeslot.Therefore, there are K 0 = 8 timeslots, and there are 4 levels of spectrum lines.Table 1 shows the timeslots and time samples that are shared between spectrum lines (X ji ) at different levels.For example, 4 samples in TS = 0 are used in the calculations of X 01 , X 11 , X 21 , and X 31 , and 4 samples in TS = 4 are used in the calculations of X 01 , X 12 , X 23 , and X 35 .There are 8, 4, 2, and 1 frequency blocks at Level-0, -1, -2, and -3, respectively.Table 2 shows that the Level-0 M-point FFT results can also be used to generate other-level spectrum lines, as long as the samples in a timeslot 1. are shared between these spectrum lines (see Table 1), and 2.
are multiplied by the same phase rotation factors that meet the criterion described in (15).

Comparisons of MFR-BFFT Results with FFT Results
Figure 1 shows a segment of I/Q data recorded by an ultra-wideband digital receiver.There is a total of 2 20 samples at Level-0, and M = 256 samples in each timeslot.Four spectrum levels are calculated with 2 20 , 2 19 , 2 18 , and 2 17 frequency bins at each of the different levels.The ADC sampling rate is 533.333MHz, so that the frequency resolutions are 508.6,1017.3,2034.5, and 4069.0Hz for Level-0 through Level-3, respectively.

Comparisons of MFR-BFFT Results with FFT Results
Figure 1 shows a segment of I/Q data recorded by an ultra-wideband digita There is a total of 2 20 samples at Level-0, and  = 256 samples in each time spectrum levels are calculated with 2 20 , 2 19 , 2 18 , and 2 17 frequency bins at e different levels.The ADC sampling rate is 533.333MHz, so that the frequency r are 508.6,1017.3,2034.5, and 4069.0Hz for Level-0 through Level-3, respectivel Figure 2 shows the MFR-BFFT results at different spectrum levels.To verify the MFR-BFFT method, the results obtained using the MATLAB "fft" function are overlaid in the figure (with an arbitrary offset in order to improve visibility).The error distributions between the results in Figure 2 are shown in Figure 3 using histogram plots.The very small errors, which are due to differences in the order in which numerical rounding has occurred, show the correctness of our MFR-BFFT method.Figure 2 shows the MFR-BFFT results at different spectrum levels.To verify the MFR-BFFT method, the results obtained using the MATLAB "fft" function are overlaid in the figure (with an arbitrary offset in order to improve visibility).The error distributions between the results in Figure 2 are shown in Figure 3 using histogram plots.The very small errors, which are due to differences in the order in which numerical rounding has occurred, show the correctness of our MFR-BFFT method.

Conclusions
This communication introduces a new multi-resolution spectrum line calculation method, called the multi-frequency-resolution blocking fast Fourier transform (MFR-BFFT).It can simultaneously calculate a set of very long multi-resolution signal spectra by reusing the BFFT results on all samples used to calculate the finest frequency resolution spectrum.The paper gives the detailed derivation of the MFR-BFFT method and uses actual recorded data to demonstrate the use of the method.The MFR-BFFT results are also compared with conventional FFT results.They are identical for all practical purposes.The success of the MFR-BFFT in simultaneously calculating multi-frequency-resolution spectrum lines using the same set of samples leads to a new and efficient approach for analyzing long-length digitized samples in real-time systems, such as FPGA-based hardware.Although this paper discusses the MFR-BFFT method for applications with long-length signals, this method can be considered for any applications needing multi-resolution spectra and using FPGAs for near-real-time calculations.In follow-on work, the MFR-BFFT method will be combined with the PCW-BFFT method to develop a new DFT-based timefrequency distribution calculation method that can concurrently produce multi-frequency-resolution and multi-time-resolution spectrograms, while consuming digitized data in  smaller blocks of  samples produced by high-speed ADCs.

Conclusions
This communication introduces a new multi-resolution spectrum line calculation method, called the multi-frequency-resolution blocking fast Fourier transform (MFR-BFFT).It can simultaneously calculate a set of very long multi-resolution signal spectra by reusing the BFFT results on all samples used to calculate the finest frequency resolution spectrum.The paper gives the detailed derivation of the MFR-BFFT method and uses actual recorded data to demonstrate the use of the method.The MFR-BFFT results are also compared with conventional FFT results.They are identical for all practical purposes.The success of the MFR-BFFT in simultaneously calculating multi-frequency-resolution spectrum lines using the same set of samples leads to a new and efficient approach for analyzing long-length digitized samples in real-time systems, such as FPGA-based hardware.Although this paper discusses the MFR-BFFT method for applications with long-length signals, this method can be considered for any applications needing multi-resolution spectra and using FPGAs for near-real-time calculations.In follow-on work, the MFR-BFFT method will be combined with the PCW-BFFT method to develop a new DFT-based time-frequency distribution calculation method that can concurrently produce multi-frequency-resolution and multitime-resolution spectrograms, while consuming digitized data in K smaller blocks of M samples produced by high-speed ADCs.

Figure 1 .
Figure 1.A segment of a digitized I/Q signal with a total of 2 20 samples.

Figure 1 .
Figure 1.A segment of a digitized I/Q signal with a total of 2 20 samples.

Figure 2 .
Figure 2. Spectrum lines calculated using different numbers of samples using the MFR-BFFT (blue) and the MATLAB "fft" function (red).Since results are very close, an arbitrary offset of 100 was added to the magnitude of the MATLAB results in order to improve visibility.

Figure 2 . 10 Figure 3 .
Figure 2. Spectrum lines calculated using different numbers of samples using the MFR-BFFT (blue) and the MATLAB "fft" function (red).Since results are very close, an arbitrary offset of 100 was added to the magnitude of the MATLAB results in order to improve visibility.

Author Contributions:
The conceptualization, algorithm, funding, and original draft preparation tasks were carried out by C.W. The concept discussions, reviewing, and editing were carried out by M.L.All authors have read and agreed to the published version of the manuscript.Funding: This research was supported by Defence Research and Development Canada-Ottawa Research Centre, Department of National Defence, Government of Canada.Institutional Review Board Statement: Not applicable.

Figure 3 .
Figure3.Histograms of the differences between the results in Figure2obtained using the MFR-BFFT method and the MATLAB "fft" function at different frequency resolution levels.

Table 1 .
Timeslots (TSs) and time samples used in different levels spectrum line (SL) calculations.

Table 2 .
The M-point FFT results can be reused at different levels