A fast DFT method for generally k sparse signals recovery

In this paper, a new method to fast compute DFT of generally sparse signals is presented. Firstly, the original signal is downsampled with different time shifts, and the discrete Fourier Transforms (DFTs) of downsampled signals are calculated by FFT. Then the DFTs are used to determine and measure the K non-zero (significant) freq. grids by combining the moment preserving problem with the BigBand method. The proposed method is hardware-friendly, and simulation results show that the proposed method has better recovery performance compared with other methods.


Introduction
Nowadays, Fast Fourier transform (FFT) is the significant approach in fast computing Discrete Fourier Transform (DFT) of the signals depending on its time complexity O(N log N), where N denotes the signal length. Recently, FFT is widely used all over the world for communications and signal processing [1][2][3]. However, it's still a big challenge to how FFT can be outperformed and receives attention. Sparsity is the main feature of the signals to speed up FFT in the literature. Signal with N length is known as exactly K-sparse, where K non-zero frequencies with K < N. Secondly, the signal is generally K-sparse if all frequencies are non-zero [4,5]. But we were only concerned about keeping the first K-largest (essential) frequencies in term of magnitudes and ignored the remainder, instead of estimating all frequencies [6].
Previously, H. Hassanieh et al. proposed the sparse Fast Fourier Transform (sFFT) [7][8][9]. The main idea of sFFT was to sample fewer frequencies (proportional to K) since most of the frequencies are zero or insignificant. In [9], BigBand was designed for the typical case of spectrum usage, where the occupied frequencies were randomly distributed with sparsity K = O( √ N) whereas sFFT proved to be a worst-case distribution of occupied frequencies (sparsity), K = O(N) . BigBand was designed under fewer constraints than sFFT so that it can be utilized more effectively compared to sFFT. Importantly, BigBand works with off-theshelf low-speed analog to digital converters [9]. The impression behind sFFT is to sample less frequency (shortened as freq. hereafter) grids (proportional to K) instead of Keeping all freq. grids since most freq. grids are zero or insignificant. FFT based on such a downsampling strategy which will only cost O(K log K). Nevertheless, as the locations and values of the K non-zero freq. grids are unknown, downsampled freq. grids frequently lead to data loss and cannot achieve perfect reconstruction. Deal with this difficulty, sFFT utilizes the strategies of filtering and permutation, which can increase the probability of capturing useful information from downsampled freq. grids. sFFT [10] is state of the art, but there are some limitations, which is followed.1) Filtering and permutation operated on x and its related to N. The implementation of sFFT for generally k-sparse signals is very complicated, and it involves so many parameters that are difficult to set.
Recently, Ghazi et al. in MIT proposed another sFFT version [6] that costs O(KlogK) operations for exactly K-sparse signals. The basic idea is similar, first downsample original signals before recovering K non-zero freq. grids from the downsampled signals via error correction techniques, where [6] uses Reed-Solomon code, which is equivalent to the Moment-preserving problem considered in [7]. The key difference is that Ghazi et al. 's method recovers all K non-zero freq. grids once.
Hsieh et al. [10] proposed a new concept about sFFT by downsampling in the time domain (sFFT-DT)for the exactly K-sparse signals, assuming that distribution of the non-zero frequency grid is uniform. Their focus was to downsample the original input signal at the beginning; and then, directly conducts operations on downsampled signals, where the length of downsampled signals was kept proportional to O(K ) . Downsampling probably leads to "aliasing," which means that different frequency grids of the original signal map into the same grid of the downsampled signal. In [10][11][12], the aliasing problem is reformulated as a momentpreserving problem (MPP) [13][14][15] and solved via existing approaches.
We proposed a new fast DFT method for generally sparse signals. The proposed method consists of three steps. First is downsampling of the original signal . Second is the calculation of the discrete Fourier Transform (DFT) of downsampled signals by using FFT, and third is the use of DFTs to determine and measure the K non-zero (significant) freq. grids . Downsampling will lead to "aliasing, " where different frequencies become indistinguishable in terms of their positions and values. According to, possible positions of nonzero freq. grid can be obtained from the solution of MPP. Based on the possible positions, several over-determined systems can be constructed, and the accurate positions and values of the aliasing terms can be determined by solving these systems. For generally sparse signals, the proposed method has a better recovery performance than the sFFT-DT method [10], and its complexity is reduced compared with the SFFT-DT method [10].

Methodology of proposed method
The proposed method consists of three steps. From Eq. 1 it can be observed that each frequency grid d is the sum of d terms . When more than two terms on the right side of Eq. 1 are non-zero, "aliasing "occurs ( Fig. 1).
For the solution of the aliasing problem, the shift property of DFT plays an important role [10]. Let where l represents the shift factor. Each freq. grid of d,l , which is the DFT of d,l , can be described as When aliasing occurs, we shall consider Eq. 2 in another aspect. Let the product of known values d and X d,l [k] be denoted by m l , i.e. m l = dX d,l [k] . Let p j and z l j describe the unknown X[s j ] and e i2 s j l∕ N , respectively, with 1 ≤ j ≤ a and Then (2) can be written as (2) (3) m l = p 1 z l 1 + p 2 z l 2 + ⋯ + p a z l a , 0 ≤ l ≤ 2a − 1

Fig. 1 Aliasing and its iterative solver
It is noticeable that m l it is the "l"th moment with m l = ∑ a j=1 p j z l j . The problem of solving p j 's and z j 's with given different moments ( m l 's) is known as the MPP. Additionally, MPP is also equivalent to the error locator polynomial problem [13] reported by Ghazi et al. 's sFFT, which is a commonly used procedure in Reed-Solomon decoding [12]. Based on orthogonal polynomials [14], the solution of MPP is provided [15], and complete analytical solution for a ≤ 4 was derived by Tsai et al. [15]. For example, the complete analytical solution a = 2 is Then the unknown positions s j 's and values X [s j ] of the aliasing non-zero terms on the right side of Eq. 2 can be determined from p j and z l j with j = 1, 2.

Proposed method
In this article, we propose a new method to fast compute DFT of generally K sparse signals. The schematic illustration of our proposed method a = 2 is demonstrated in Fig. 2, and this method is feasible for a ≤ 4 the analytic solutions of MPP. For a = 2 , we use three time shifts to get four downsampled signals d,l with 0 ≤ l ≤ 3.The DFT d,l of the downsampled signals can be expressed as where s 1 and s 2 denote positions of the two aliasing terms, and X [s 1 ], X[s 2 ] denote their values.
According to the above section, we can estimate the large, the solution of MPP can approach to the real value of (s 1 , s 2 ) . However, when the original signal is with low SNR, the estimated positions (ŝ 1 ,ŝ 2 ) based on the MPP method might no longer belong to the set S k . Let s ′ 1 and s ′ 2 denote the values in S k that nearest to ŝ 1 and ŝ 2 respectively. Then we can consider that the real positions (s 1 , s 2 ) satisfy.
For simplicity, define 1 Then we can get nine possible (s 1 , s 2 ) pairs, by choosing one value from 1 and choose the other one from 2 . In order to determine the accurate positions (s 1 , s 2 ) , we use the possible pairs to construct over-determined systems. Then the accurate positions and values of the aliasing terms can be recovered by solving these systems [8].
For each possible pair of (s 1 , s 2 ) , Eq. 5 becomes an overdetermined linear system of equations, in which there exist four equations and two unknowns X[ T . Then Eq. 5 can be written into matrix form as = , where As the matrix A is with full column rank, the unknown z can be solved by ̂ = † , where † = ( H ) −1 H is the Moore-Penrose pseudo-inverse of A.
The recovery error can be defined as e = ‖ −̂‖ 2 . After solving the over-determined system for each possible pair, the accurate positions (s 1 , s 2 ) , as well as the corresponding values X[s 1 ], X [s 2 ] , can be determined by finding the minimum recovery error.
In the BigBand method proposed by H. Hassanieh et al. [9], there are d(d − 1)∕ 2 possible (s 1 , s 2 ) pairs, which means that d(d − 1)∕ 2 over-determined systems need to be solved. Meanwhile, with the assistance of the solution to MPP, only nine possible pairs need to be considered in our method. As normally d >> 5 , our method requires less digital computation compared to the BigBand method and sFFT-DT.

MMP error analysis problem
When we are dealing with generally k sparse signal, the solved roots no longer belong to the set S k of real roots, as mentioned above. Which correspond to candidate correct  Figure 3 shows the real roots (in red) and solves roots (in blue) in complex coordinates under different SNR values. The real roots must locate with in the unit circle with radius = 1, sine they belong to S k , and each element of S k satisfies Nevertheless, the solved roots could deviate from the real roots and cannot be located within the unit circles. In order to close the gap between the candidate positions and estimated positions, the proposed method uses the minimum mean square method.

Proposed algorithm flow chart
See Fig. 4.

Result and discussion
During simulations in the following work, the signal , in time domain was generated as followed. Firstly, generate a frequency-domain signal 0 , which consists of K non-zero entries and N − K zeros where N = 1024, Then is obtained by adding the inverse FFT of 0 with white Gaussian noise.
Simulations are conducted to compare the proposed method with SFFT-DT and BigBand in terms of recovery performance and computational time. For all three methods, the downsampling factor d = 16 is fixed, and we suppose that the maximum number of non-zero aliasing terms is a m = 2 . Define successful rate as the fraction of occupied frequencies, which are successfully recovered. For different sparsity K, the successful rate of the proposed method versus SNR has been measured and reported in Fig. 5. Experimental results show that the recovery performance improves gradually with the decrease of K, as larger K causes more aliasing situations. Figure 6 shows the results of reconstruction accuracy versus SNR for the three methods mentioned before with K = 25. It is observed that our proposed method algorithm has a better recovery performance with the other two methods. In the BigBand method proposed by H.Hassanieh et al. [9], there are d(d − 1)/2 possible pairs, (s 1 , s 2 ) which mean that d(d-1)/2 over-determined systems need to be solved. Meanwhile, with the assistance of the solution to MPP, only nine possible pairs are needed to be considered in our method. As normally d ≫ 5, our method requires less digital computation compared to BigBand as shown in Fig. 7].
The false-alarm probability is defined as the fraction of empty frequencies, incorrectly reported as occupied. Figure 7 shows that the proposed method has lower false-alarm rate compared with the other two methods under different SNRs. The false-alarm probability of proposed method stays below 0.04 with SNR ≥ 25 dB (Fig. 8).

Conclusion
For generally sparse signals, this paper presents a new fast computing DFT method. This method is based on downsampling and combines MPP with the BigBand method. All operations of solving MPP are linear with analytical solutions involved. The BigBand method is utilized to modify the error of the solution to MPP, caused by the interference of insignificant freq. grids. Theoretical complexity analysis and simulation results demonstrate that the proposed method is hardware-friendly and has a better performance as compared to other reported algorithms.