{'fraction_non_alphanumeric': 0.04707049684179957, 'fraction_numerical': 0.009752910142542533, 'mean_word_length': 4.14147069802165, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse Fast Fourier Transform for Exactly and Generally K-Sparse Signals by Downsampling and Sparse Recovery

Fast Fourier Transform (FFT) is one of the most important tools in digital signal processing. FFT costs O(N log N) for transforming a signal of length N. Recently, researchers at MIT have proposed Sparse Fast Fourier Transform (sFFT) [1][2] as a breakthrough with computational complexity O(K log N) and O(K log N log N/K) for exactly K-sparse signal (with only K non-zero frequency grids) and generally K-sparse signal (with K significant frequency grids), respectively, to outperform classic FFT. In this paper, a new sparse Fast Fourier Transform by downsampling in the time domain (sFFT-DT) is proposed, based on the assumption that the distribution of the non-zero frequency grids is uniform. The idea behind sFFT-DT is to downsample the original input signal at the beginning; then, subsequent processing operates under downsampled signals, where signal lengths are proportional to O(K). Downsampling, however, possibly leads to"aliasing". By the shift property of DFT, we recast the aliasing problem as a"moment-preserving problem (MPP),"which is solvable. We prove two theorems related to initializing the downsampling factors under different conditions to have computational complexity, O(K log K) and O(K^(5/4) log K). Moreover, for generally K-sparse signals, solutions to the MPP are inaccurate due to interference from non-significant frequency grids. We observe that significant frequency grids in aliasing are"sparse". This property is exploited, and a new sparse signal recovery algorithm in the context of compressive sensing is presented to refine the solution to MPP. The computational complexity still costs O(K log K) but requires a larger Big-O constant compared to the case of exactly K-sparse signals. We conduct theoretical complexity analyses and simulations to demonstrate that our method (sFFT-DT) outperforms classic FFT and MIT's sFFT.


I. INTRODUCTION
A. Background [2], which was proved to outperform FFT. sFFT relies on the assumption that input signals exhibit "sparsity" in the frequency domain. Let x ∈ C N be the input signal in the time domain, and let X ∈ C N be the Discrete Fourier Transform (DFT) of x. x is exactly K-sparse if there are K non-zero entries in X and K < N . On the other hand, x is a generally K-sparse signal if there are more than K non-zero entries in X, but it is desirable to keep only the first K-largest (significant) entries in terms of magnitude and ignore the remainder. For both types of signals, sFFT costs O(K log N ) and The idea behind sFFT is to sample fewer frequency (abbreviated 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 subsampling strategy will only cost O(K log K) calculations ideally. Nevertheless, as the locations and values of the K non-zero entries are unknown, subsampled freq. grids often lead to data loss and cannot achieve perfect (100%) reconstruction. To cope with this difficulty, sFFT utilizes the strategies of filtering and permutation, which can increase the probability of capturing useful information from subsampled freq. grids. These operations cost O(K log N ). According to [1] [2], sFFT is faster than FFTW [3] (a very fast C subroutine library for computing FFT) when X is exactly K-sparse with K ≤ N 2 6 . sFFT also outperforms previous works, such as [4] [5].
Nevertheless, for a generally K-sparse signal, subsampled freq. grids certainly face interference from neighboring insignificant freq. grids. To alleviate this interference, which is assumed to be random noise in [1], sFFT permutes X randomly in order to make neighboring entries randomly distributed. In this way, sFFT produces a number of subsampled signals under different permutations, where each subsampled signal is considered a candidate and the median is identified from these candidates.
Even though sFFT [1] [2] is outstanding, there are some limitations, summarized as follows. 1) Filtering and permutation are operated on x. Since x ∈ C N , these operations are related to N . Thus, the complexity of sFFT still involves N and cannot achieve the theoretical ideal complexity O(K log K). 2) sFFT only guarantees that it succeeds with a constant probability (e.g., 3) The implementation of sFFT for generally K-sparse signals is very complicated as it involves too many parameters that are difficult to set. † Recently, Ghazi et al. in MIT proposed another sFFT version [6] that costs O(K log K) operations for exactly K-sparse signals. The basic idea is similar to our previous work [7]. Both methods 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, while we propose a top-down strategy to solve K non-zero freq. grids iteratively.
The comparison between these two methods in terms of computational complexity and successful probability will be discussed later in Sec. II-D.
In addition to a series of results regarding sFFT from MIT, Pawar and Ramchandran [8] proposed an algorithm, called FFAST (Fast Fourier Aliasing-based Sparse Transform), which focuses only on exactly K-sparse signals. Their approach is based on filterless downsampling of the input signal x using a small set of co-prime downsampling factors guided by the Chinese Remainder Theorem (CRT). These aliasing patterns of different downsampled signals are formulated as parity-check constraints of good erasure-correcting sparse-graph codes. They also show that computing X is equivalent to decoding these sparse-graph codes. As claimed in [8], FFAST costs O(K log K) and outperforms sFFT [1] [2]. Nevertheless, FFAST depends on co-prime downsampling factors. In other words, the signal length N must be a product of a number (typically 3 to 9) of distinct primes; otherwise FFAST cannot work. Moreover, the smallest downsampling factor bounds FFAST's computational cost. For example, if N = 2 20 3 2 and K = 2 16 , the smallest downsampling factor is 3 2 . In this case, the computational cost of calculating FFT of a downsampled signal with length N/3 2 is (far) higher than O(K log K). Actually, these limitations are possibly harsh.

B. Our Contributions
In our previous work [7], we proposed a new algorithm, called sFFT-DT, by downsampling in the time domain with time complexity of O(K log K) only for exactly K-sparse signals. The idea behind sFFT-DT is to downsample the input signal in the time domain before directly conducting all subsequent operations on the downsampled signals. By choosing an appropriate downsampling factor to make the length of a downsampled signal be O(K), no operations related to N are required in sFFT-DT. Downsampling, however, possibly leads to "aliasing," where different freq. grids become indistinguishable in terms of † In fact, according to our private communication with the authors of [1][2], they would not recommend implementing this code since it is not trivial. The authors also suggest that it is not easy to clearly illustrate which setting will work best because of the constants in the Big-O functions and because of the dependency on the implementation. The authors themselves did not implement it since they believed that the constants would be large and that it would not realize much improvement over FFTW.
their locations and values. To overcome this problem, the locations and values of these K non-zero entries are considered as unknown variables and the "aliasing problem" is reformulated as "moment-preserving problem (MPP)," which can be solved via existing approaches, such as that in [9]. Furthermore, sFFT-DT is conducted in a manner of a top-down iterative strategy under different downsampling factors, which can efficiently reduce the computational cost.
In this paper, the accurate computational cost and theoretical performance lower bound of sFFT-DT are proven further for exactly K-sparse signals. We also derive the Big-O constants of computational complexity of sFFT-DT and show that they are smaller than those of Ghazi et al.'s sFFT [6]. In addition, all operations of sFFT-DT are solved via analytical solutions but those of Ghazi et al.'s sFFT involve a numerical root finding algorithm, which is more complicated in terms of hardware implementation.
In previous works [1][2][6] [8], performance and computational complexity have been analyzed based on the assumption that sparsity K is known in advance. In practice, however, K is unknown and an input parameter decided by the user. If K is not guessed correctly, the performance is degraded and/or the computational overhead is higher than expected because the choice of some parameters depends on K. In this paper, we also propose a very simple solution to address this problem and relax this impractical assumption. We show that the cost for deciding K is the same as that required for sFFT-DT with known K.
Different from our previous study [7], we further extend sFFT-DT to generally K-sparse signals in this paper, in addition to conducting more advanced theoretical analyses. For generally K-sparse signals, since almost all freq. grids are non-zero, each freq. grid of a downsampled signal is composed of significant and insignificant entries due to aliasing from downsampling.
To extract significant entries from each freq. grid, the concept of sparse signal recovery from fewer samples, originating from compressive sensing (CS) [10], is employed since significant entries are "sparse". Under this circumstance, the insignificant freq. grids act like Gaussian random noise due to central limit theory and the signal-to-noise ratio (SNR) between the significant and insignificant freq. grids affects the successful probability of sparse signal recovery.
Finally, we conclude that our methods are easy to implement and are demonstrated to outperform the state-of-the-art in terms of theoretical analyses and simulation results.

C. Organization of This Paper
The remainder of this paper is organized as follows. In Sec. II, we describe the proposed method for exactly K-sparse signals. Our method for generally K-sparse signals will be expounded in Sec. III. Conclusions will be provided in Sec. IV.

A. Problem Formulation
The objective of this paper is to locate and estimate K non-zero freq. grids of X from the DFT X d of the downsampled signal x d .
Note that each freq. grid of X d is a sum of d terms of X. When more than two terms of X are non-zero, "aliasing" occurs, as illustrated in Fig. 1. Fig. 1(a) shows an original signal in the frequency domain, where only three freq. grids are non-zero (appearing at normalized frequencies = 0π, 0.5π, and π). Fig. 1(b) shows the downsampled signal in the frequency domain when d = 2, where the downsampled freq. grid at 0π incurs aliasing and is the sum of the original freq. grids at 0π and π, while no aliasing occurs at π. In Fig. 1(b), we solve all non-zero downsampled freq. grids once, no matter whether aliasing occurs or not. This procedure is called non-iterative sFFT-DT and will be discussed in detail later. Instead of solving all of the downsampled freq. grids once, Fig. 1(c) illustrates an example of iteratively solving freq. grids. At the first iteration, the downsampled freq. grid without aliasing at 1π is solved. This makes the remaining downsampled freq. grids more sparse.
Then, the signal is downsampled again with d = 4. At the second iteration, we solve the downsampled freq. grid with aliasing at 0π. This procedure, called iterative sFFT-DT, will be discussed further in Sec. II-E. grid of X d,l is denoted as: Thus, Eq. (2) degenerates to Eq. (1) when l = 0.
In practice, all we can obtain are X d,l [k]'s for different l's. For each downsampling factor d, there will be no more than d terms on the right side of Eq. (2), where each term contains two unknown variables. For example, X[k]e i2πkl/N is composed of two variables, X[k] and e i2πkl/N . Let a, 1 ≤ a ≤ d, denote the number of terms on the right side of Eq. (2). Therefore, we need 2a equations to solve these 2a variables, and l is within the range of [0, 2a − 1]. By taking the above into consideration, the problem of solving the 2a unknown variables on the right side of Eq. (2) can be formulated as: where X d,l [k] is known and denoted as m l while p j and z l j represent unknown X[s j ] and e i2πsj l/N , respectively, for s j ∈ To simplify the notation, we let S k = {k, k + N d , ... , k + (d − 1) N d } and It is trivial that no aliasing occurs if a = 1, irrespective of the downsampling factor. Under this circumstance, we have al.'s sFFT, which is a commonly used step in Reed-Solomon decoding [11]. The solution to MPP [9][12] based on orthogonal polynomials is useful and is discussed in the next subsection.

B. The Solution to Moment-Preserving Problem
Note that the moment-preserving problem (Eq. (3)) is nonlinear and cannot be solved by simple linear matrix operations.
On the contrary, we have to solve z j 's first, such that Eq. (3) becomes linear. Then, p i 's can be solved by matrix inversion.
Thus, the main difficulty is how to solve z j 's given known moments. According to [9], given the unique moments with m 0 , m 1 , ..., m 2a−1 , there must exist the corresponding orthogonal polynomial equation, P (z), with roots z j 's for 0 ≤ j ≤ a − 1.
That is, z j 's can be obtained as the roots of P (z). The steps for solving MPP are as follows.
Step (i): Let the orthogonal polynomial equation P (z) be: The relationship between P (z) and the moments is as follows:  Step (ii): Find the roots of P (z) in Eq. (4). These roots are the solutions of z 0 , z 1 ,...z a−1 , respectively.
Step (iii): Substitute all z j 's into Eq. (3), and solve the resulting equations to obtain p j 's.
Tsai [12] proposed a complete analytical solution composed of the aforementioned three steps for a ≤ 4, based on the constraint that p 0 + p 1 + ... + p a−1 = 1. Nevertheless, for the aliasing problem considered here, the constraint is , as indicated in Eq. (2). Thus, we derive the complete analytical solution for a = 2 ‡ as: The analytical solutions for an univariate polynomial with a ≤ 4 cost O(a 2 ) operations. Since there are N d freq. grids, the . For a > 4, Step (i) still costs O(a 2 ), according to the Berlekamp-Massey algorithm [13], which is well-known in Reed-Solomon decoding [11]. In addition, Step (iii) is designed to calculate the inverse matrix of a Vandermonde matrix and costs O(a 2 ) [14]. There is, however, no analytical solution of Step (ii) for a > 4. Thus, numerical methods of root finding algorithms with finite precision are required. A fast algorithm proposed by Pan [15] can approximate all of the roots with O(a(log log N ) O(1) ), where the detailed proof was shown in [6].
Step (ii) will dominate the cost of solving the MPP.
It is noted that the actual number of aliasing terms for each freq. grid, a, is unknown in advance. In practice, we let a m be user-defined and denoted as the maximum number of aliasing terms for all downsampled freq. grids during the whole ‡ The solutions for a = 3 and a = 4 are described in the Appendix.
process of solving the MPP. Therefore, 2a m moments are required to solve the MPP. If a's of all downsampled freq. grids are smaller than or equal to a m , the MPP perfectly recovers all of the freq. grids; i.e., it resolves all non-zero values of X and their corresponding locations. Otherwise, some non-zero entries of X cannot be recovered due to insufficient information.
Although a larger a m guarantees better recovery performance, it also means that more moments and higher computational cost are required.
In sum, the cost of the MPP consists of two parts. Since the size of a downsampled signal is N d , the cost of generating the where either of which is defined as the "P2 cost of MPP".
So far, our method of solving all downsampled freq. grids is based on fixing downsampling factor d (and a m ), as an example illustrated in Fig. 1 (b). In this case, we call this approach, non-iterative sFFT-DT. Its iterative counterpart, iterative sFFT-DT, which is the focus of this paper, will be described later in Secs. II-E and II-F.

C. MPP in sFFT-DT and Ghazi et al.'s sFFT
On the other hand, Ghazi et al.'s sFFT is to set d = N log K K and a m = C log K (C is a constant), and the total cost is ). In fact, these two kinds of parameter settings lead to different computational cost and recovery performance. We will discuss this issue further in the next subsection.
Moreover, in order to reduce computational cost in our framework further, a top-down iterative strategy for sFFT-DT will be presented in Sec. II-E.

D. Analysis
In this section, we first will study the relationship between a m and d and analyze the probability of a downsampled freq.
grid with an aliasing number larger than a m . Second, we will discuss computational complexity and recovery performance of our non-iterative sFFT-DT. Third, we will compare our method with Ghazi et al.'s sFFT [6]. In addition, the Big-O constant of complexity is induced in order to highlight the computational simplicity of non-iterative sFFT-DT. Finally, we will conclude by presenting an iterative sFFT-DT to reduce computational cost further.

1) Relationship between Maximum Aliasing Number and Downsampling
Factor: Now, we consider the relationship between a m and d. If a m is set to d, then we always can recover any X without errors but the computational cost will be larger than that of FFT. Thus, it is preferable to set 1 ≤ a m < d. In fact, small a m is feasible when X is uniformly distributed. For each freq. grid, the number of aliasing, a, will be small with high probability if N dK is larger than 1.
) am+1 , where e is Euler's. The probability that X can be perfectly reconstructed is Proof: For each downsampled freq. grid, the probability of a > a m is Σ d . Under this circumstance, the probability of at least a downsampled freq. grid with a > a m is bounded by . Thus, we can derive: 2) Computational Cost and Recovery Performance: From Lemma 1, non-iterative sFFT-DT obtains perfect recovery with probability δ = 1 − P r(d, a m ) at least. Now, we first discuss our approach with fixed a m = 4, under which most freq. grids can be solved, as an example illustrated in Fig. 2. The computational cost and recovery performance of non-iterative sFFT-DT are depicted further in Theorem 1 and Theorem 2, respectively, based on different settings of d.

Theorem 1.
If non-zero freq. grids of X distribute uniformly, given a m = 4 and δ, the upper bound of d is

sFFT-DT perfectly recovers X with probability at least δ and the computational cost is
Proof: According to Eq. (7), we can derive: In addition, given d = O( N K ), the computational cost can be derived to be O(K log K), as described in Sec. II-B.
In summary, Theorem 1 guarantees perfect recovery with probability being larger than δ and computational cost being bounded by O(K 5 4 log K), whatever K is, while Theorem 2 achieves the ideal computational cost O(K log K) with recovery performance being degraded along with the increase of K.
3) Comparison with [6]: We will compare our results in Theorem 1 and Theorem 2 with the one described in Theorem 3 below.
Theorem 3. (Rewritten from Theorem 4.6 [6]) If non-zero freq. grids of X distribute uniformly, given a m = C log K and d =  solution is easier to implement than a numerical procedure. We conclude that there are two main advantages in Theorem 1 and Theorem 2, compared to Theorem 3 [6]. First, the Big-O constants in Theorem 1 and Theorem 2 are smaller than those of Ghazi et al.'s sFFT. Second, our analytical solution is hardware-friendly in terms of implementation.
On the other hand, when the signal is not so sparse with K approaching N (e.g., K = N 8 and N + = 2 0 ), the cost of 8 FFTs in a downsampled signal is almost equivalent to that of one FFT in the original signal. To further reduce the cost, a top-down iterative strategy is proposed in the next subsection.
It also should be noted that the above discussions (and prior works) are based on the assumption that K is known. In practice, K is unknown in advance. Unfortunately, how to automatically determine K is ignored in the literature. Instead of skipping this problem, in this paper, we present a simple but effective strategy in Sec. II-H to address this issue.

E. Top-Down Iterative Strategy for Iterative sFFT-DT
According to Fig. 2, the probability of aliasing decreases fast with the increase of a when N + ≥ 1. Thus, we can solve these downsampled freq. grids with a = 1 before subtracting these solved freq. grids from X d . This will make X d more sparse. Then, under this circumstance, d subsequently can be set to be larger values to reduce computational cost without sacrificing the recovery performance. Recall that, when d is increased, a is increased as well. In this section, an iterative strategy is proposed to solve the aliasing problem with an iterative increase of the downsampling factor d. Fig 1 illustrates such an example. In Fig. 1(b), if we try to solve all aliasing problems in the first iteration, 4 FFTs are required, since the maximum value of a is 2. On the other hand, if we first solve the downsampled freq. grids with a = 1 (at normalized frequency = π), it costs 2 FFTs, as shown in Fig. 1(c). Since 2 FFTs are insufficient for solving the aliasing problem completely under a = 2, extra 2 FFTs are required to solve a more sparse signal.
The key is how to calculate the 2 extra FFTs in the above example with lower cost. The idea motivated by sFFT is to discard the solved freq. grids in the frequency domain. Thus, if, for example, K ′ freq. grids are subtracted from the original X, the sparsity of the remaining signal is K − K ′ . Since a more sparse signal is generated in an iterative manner, d can be set to be larger under fixed N + . As shown in Fig. 1(d), 2 extra FFTs can be done quickly with a larger d (=4) to solve the downsampled freq. grid with a = 2 (at normalized frequency = 0π). Consequently, d is doubled iteratively in our method and the total cost is dominated by that required at the first iteration.
The proposed method with the top-down iterative strategy is called iterative sFFT-DT.

F. Iterative sFFT-DT: Algorithm for Exactly K-Sparse Signals
In this section, our method, iterative sFFT-DT, is developed based on Theorem 2 and is depicted in Algorithm 1, which is composed of three functions, main, SubFreq, and MPP. In fact, the only difference between Theorem 1 and Theorem 2 is how to initialize d. Thus, the algorithm iterative sFFT-DT is modified easily to satisfy Theorem 1. Basically, iterative sFFT-DT solves downsampled freq. grids from a m = 1 to 4 with an iterative increase of d. Note that, its variation, non-iterative sFFT-DT, solves all downsampled freq. grids with a m = 4 and fixed d.
At the initialization stage, the sets S and T , recording the positions of solved and unsolved freq. grids, respectively, are set to be empty. a m and d are initialized via Theorem 2. The algorithmic steps are explained in detail as follows.
Function main, which is executed in a top-down manner by doubling the downsampling factor iteratively, is depicted from Line 1 to Line 16. In Lines 3-4, the input signal x is represented by two shift factors 2l and 2l + 1. Then they are used to perform FFT to obtain X d,2l and X d,2l+1 in Lines 5-6. In Line 7, the function SubFreq, depicted between Line 17 and Line 22, is executed to remove freq. grids from X d,2l and X d,2l+1 that were solved in previous iterations. The goal of function SubFreq is to make the resulting signal more sparse. grids are solved under correct a's. Finally, the downsampling factor is doubled, as indicated in Line 14, to solve the unsolved freq. grids in an iterative manner. This means that the downsampled signal in the next iteration will become shorter and can be dealt faster than that in the previous iterations.

Algorithm 1 Iterative sFFT-DT for exactly K-sparse signals.
Input: x, K; Output: X; SubFreq(X d,2l , X d,2l+1 , X, d, l, S);  In sum, the proposed algorithm, iterative sFFT-DT, is dominated by "FFT" and costs O(K log K) operations. Now, we discuss Big-O constants for operations of addition and multiplication, respectively. Since d is doubled iteratively, the P1 cost of MPP gradually is reduced in the later iterations. The total cost is and those for multiplication range from 6 to 12.
As for recovery performance in iterative sFFT-DT, since the downsampling factor d is doubled along with the increase of iterations, a question, which naturally arises, is if a larger downsampling factor leads to more new aliasing artifacts. If yes, these newly generated aliasing artifacts degrade the performance of iterative sFFT-DT. If no, the iterative style is good since it reduces computational cost and maintains recovery performance.
In Lemma 2, we prove that the probability of producing new aliasing artifacts after a sufficient number of iterations will approach zero.

Lemma 2.
Suppose K non-zero entries of X distribute uniformly (i.e., with probability K N ). Let P r ali l be the probability that new aliasing artifacts are produced at the l'th iteration in iterative sFFT-DT. Let K l be the number of freq. grids with a ≥ l + 1 at the l'th iteration (0 ≤ l ≤ 3, Proof: According to Algorithm 1, after the first iteration (l = 0), all downsampled freq. grids with only a = 1 aliasing term are solved. Thus, we focus on discussing the probability of producing new aliasing artifacts under l ≥ 1. By the same idea of Lemma 1, we can define N d to be the probability that there is a downsampled freq. grid with a ≥ 2 aliasing terms. Under the condition of l ≥ 1 and a ≥ 2, however, some non-zero freq. grids have been solved in previous iterations.
Thus, the number of remaining non-zero freq. grids are no longer K and K N and should be modified. In other words, the number of downsampled freq. grids with a ≥ 2 must be less than K l for l ≥ 1 and N d According to our iterative sFFT-DT algorithm, let d l = 2 l d for l ≥ 1. We can derive: Eq. (10) converges to 0 only when K l ≤ K 2 l . As shown in Fig. 2, by initializing d properly, almost K freq. grids can be solved in the first few iterations. This makes K l ≤ K 2 l easy to be satisfied. Under this circumstance, N + ∈ [2 1 2 3 ] would be a good choice. By replacing K l with K 2 l , we can derive P r ali l ≤ 1 2 l+1 ( d N )K 2 . When l increases to be large enough, the probability of P r ali l will be small because of Lemma 2 indicates the probability of producing new aliasing artifacts in an asymptotic manner. This provides us the information that the probability of producing new aliasing artifacts finally converges to zero. In our simulations, we will show that the exact probability with new aliasing is very low if d is selected based on either Theorem 1 or Theorem 2. In other words, the iterative approach can reduce the computational cost and maintain the recovery performance when there is no new aliasing generated at later iterations.

H. A Simple Strategy for Estimating Unknown Sparsity K
As previously described, the sparsity K of a signal is important in deciding the downsampling factor d. Nevertheless, K is, in general, unknown. In this section, we provide a simple bottom-up strategy to address this issue.
First, we set a large downsampling factor d = N , and then run sFFT-DT. If there is any downsampled freq. grid that cannot be solved, then d is halved and sFFT-DT is applied to solve X again. When d is halved iteratively until the condition in either Theorem 1 or Theorem 2 is satisfied, sFFT-DT guarantees one to stop with the probability indicated in either Theorem 1 or Theorem 2. This strategy needs the same computational complexity required in sFFT-DT with known K because the cost

I. Simulation Results for Exactly K-Sparse Signals
Our method, iterative sFFT-DT, was numerically verified and compared with FFTW (http://www.fftw.org/), sFFT-v3 [    III. SFFT-DT FOR GENERALLY K -SPARSE SIGNALS For sparse FFT of a generally K-sparse signal x, the goal is to compute an approximate transform X approx satisfying: where X approx is exactly K-sparse and X is generally K-sparse. Without loss of generality, we assume that all freq. grids in X are non-zero. Similar to exactly K-sparse signals, we assume that K significant freq. grids (with the first K largest magnitudes) of X distribute uniformly.
Due to generally K-sparsity of X, the right-hand side of Eq. (1) will contain d terms. If we try to solve the MPP with a = d, the total computational cost will exceed FFT, thereby creating nonsense. Therefore, our policy is to solve the MPP with only contains two significant terms, we formulate this aliasing problem with a = 2. Nevertheless, X d [k] is, in fact, composed of d terms instead of 2 terms. In this case, from Eq. (2), the sum of d − 2 insignificant terms can be considered as a normal distribution by central limit theory and act as noise, leading to inaccurate estimation in Eq. (11). Apparently, the degree of inaccuracy depends on the signal-to-noise ratio (SNR) SN R ori between the significant terms and insignificant terms. More precisely, let X opt be the optimal solution of Eq. (11). SN R ori is defined as: where M SE(·) is the function of calculating the mean squared error, X opt is equivalent to significant terms, and X − X opt represents insignificant terms.
Due to the impact of insignificant terms, the solved roots (estimated positions), z's, of P (z) = 0 (Eq. (4)), no longer belong to the set U k of real roots, as defined in the paragraph below Eq. (3), which correspond to candidate/correct positions.
Nevertheless, if SN R ori is sufficiently large (> 20db), our empirical observations show that z's can approach the real roots. The real roots must locate within the unit circle with radius = 1 since they belong to U k and each element of U k satisfies e i2π(k+ N d )l/N 2 = 1. Nevertheless, the solved roots could deviate from the real roots and not be located within the unit circles. In order to close the gap between the candidate positions and estimated positions, a refinement strategy is proposed in terms of sparse signal recovery in the context of compressive sensing; then, it is combined with non-iterative sFFT-DT because iterative sFFT-DT cannot work well in the aforementioned scenario. More specifically, the idea behind iterative sFFT-DT is to remove the solved freq. grids in an iterative manner to make the remaining unsolved signal more sparse. Nevertheless, we cannot guarantee an exact solution can be attained at each iteration to avoid error propagation into subsequent iterations. 21 The refinement strategy by sparse signal recovery is motivated by the following facts. 1) The magnitudes of significant terms must be larger than those of insignificant terms.
2) The number of significant terms is less than that of insignificant terms.
Thus, estimating the locations and values of significant terms is reformulated as a sparse signal recovery problem.

A. Refinement via Sparse Signal Recovery
Instead of considering the locations (z l j ) and values (p j ) as variables in Eq. (3), only values (p j ) are thought of as variables here. Then, we reformulate the aliasing problem in terms of a matrix form m = W b as:  . . . . . .
where W k (t, l) = e i2π(k+t N d )l/N and t ∈ [0, d − 1], b t is the value at (k + t N d )'th freq. grid, and W ∈ C r×d . By comparing when SN R ori is large enough. With the assumption of sparsity, Eq. (13) is equivalent to a sparse signal recovery problem.
Conventionally, two strategies [10], ℓ 1 -minimization and greedy approaches, are used for sparse signal recovery. Computational complexity and recovery performance, however, are two main issues that need further attention.
First, the complexities of these standard CS approaches require polynomial time. For example, Orthogonal Matching Pursuit (OMP) [17] is one of the greedy algorithms and requires O(ard) for solving Eq. (13). OMP runs at most N d times and leads to the total cost of sFFT-DT being O(arN ). This is unacceptably slower than conventional FFT. Fig. 5 illustrates the fact that the distance between the solved root and the real root is related to SN R ori . In other words, the real roots that are far from the corresponding solved root can be ignored when SN R ori is sufficiently large (> 20db from our empirical observations).
Thus, a large number of columns corresponding to the real roots that are far from the solved roots can be deleted. Moreover, the number of remaining columns is related to the number of solved roots, which is O(a), and OMP costs O(ra 2 ) now.
Nevertheless, how to decide a for each freq. grid is difficult. In addition, the details about how the solved roots are affected will be addressed later in Sec. III-B.
Second, how to set r is important in terms of computational complexity, O(ra 2 ), and performance. Candes and Wakin [18] pointed out that r required to recover sparse signal must satisfy r ≥ Ca log s a , where C is a constant and s is the column is not related to s, then r ≥ Ca log N aK . In other words, N is also a parameter that impacts the size of r. This will make the cost of solving Eq. (13)  Third, from the theory of sparse signal recovery, the recovery performance also depends on mutual coherence of W , which is defined as: where W i is the i'th column of W and µ is expected to be as small as possible. In this case, the phase difference between the W k (t, l) and W k (t + 1, l) is 2π l d , as defined in Eq. (13). Recall l ∈ {n 0 , n 1 , ..., n r−1 }. In the exactly K-sparse case, the shift operator ranges from 0 to 2a − 1. In other words, we set n 0 = 0, n 1 = 1, ..., n r−1 = 2a − 1. Even though the maximum shift operator l = 2a − 1 is encountered, the phase difference between the W k (t, l) and W k (t + 1, l) still approaches 0 with a). Under this circumstance, µ → 1 and perfect sparse recovery will become impossible. This is why the elements of {n 0 , n 1 , ..., n r−1 } should be chosen randomly within the set [0, d − 1]. In this case, W will form a random partial Fourier matrix and its mutual coherence will be small, as shown in [19] [20].
In sum, sFFT-DT for generally K-sparse signals must satisfy the three aforementioned requirements simultaneously and perform more FFTs of downsampled signals with different shift factors. Under this situation, W ∈ C r×d in Eq. (13) is simplified to W ∈ C O(a)×O(a) and the acceptable cost of OMP for recovering sparse signals is within O(a 3 ). In experiments, we suggest r = 3a.
Consequently, 3a moments with random shift operators are required to minimize the mutual coherence of W for refining the solved roots. Furthermore, we still need the moments with shift factors ranging from 0 to 2a − 1 to obtain solved roots by MPP. In sum, the total number of moments required to solve an aliasing problem containing a significant terms is 5a. As a result, for the generally K-sparse case, it requires 5a m moments, where a m denotes the maximum number of significant terms within all of the downsampled freq. grids.

B. How to Decide the Number of Solved Roots (a) for Each Freq. Grid
The previous section for refinement via sparse signal recovery exploits the assumption that a is known for each downsampled freq. grid. In reality, though, a is unknown. Thus, deciding a for each freq. grid is an important task.
where n j = j like the setting in exactly K-sparse signals. M ns acts like a "noise" matrix produced by insignificant components.
Thus, it is apparent that M s F M ns F , where · F denotes the Frobenius norm, is proportional to SN R ori , as defined in Eq. (12). In addition, letting n j = n 0 + vj for j ∈ [0, 2a m − 2], where v is a constant (discussed later), M can also be expressed as: where  (3) The formulation of SVD is from. Consequently, we are interested in finding the first K largest singular values to resolve the K non-zero frequency grids.

C. Non-iterative sFFT-DT: Algorithm for Generally K-Sparse Signals
For exactly K-sparse signals, sFFT-DT adopts the proposed top-down iterative strategy for reducing the computational cost.
Nevertheless, it is not appropriate for generally K-sparse signals since we cannot guarantee exact solutions that can be obtained at each iteration, leading to error propagation. Thus, for generally K-sparse signals, sFFT-DT solves the aliasing problem once, as shown in Algorithm 2. Algorithm 2 sFFT-DT for generally K-sparse signals.
Input: x, K; Output: X approx ; Consequently, the computational cost of sFFT-DT for generally K-sparse signals still is impacted by a m and d as in the exactly-K sparse case. If significant freq. grids distribute uniformly, according to Fig. 2, both a m and d can be set based on Theorem 1 or Theorem 2. In this case, since a m = 4, the computational cost is bounded by Part 1, which is the same as in the exactly K-sparse case. It should be noted that the Big-O constants of the generally K-sparse case are larger than those of the exactly K-sparse case because the former needs more moments.

E. Simulation Results for Generally K-Sparse Signals
The simulation environment is similar to the one described in Sec. II-I. We only compare sFFT-DT with FFTW because sFFT [1] [2] does not release the code and the code of sFFT for the generally K-sparse case is difficult to implement (as mentioned in the footnote on Page 3). Therefore, no experimental results for generally K-sparse signals are shown in their papers or websites.
Here, the test signals were generated from the mixture Gaussian model as: where p = K N is the active probability that decides which Gaussian model is used and σ on > σ of f . For each test signal, its real root is defined as X opt , as described in Sec. III, and X approx is the output signal obtained from sFFT-DT. Due to the fixed noise model, when N K becomes larger, SN R ori decreases accordingly.
The comparison of computational time is depicted in Fig. 6. Fig. 6(a) shows the results of computational time versus signal sparsity under fixed N . It is observed that sFFT-DT is remarkably faster than FFTW, except for the cases with K ≥ 2 15 . Fig.   6(b) shows the results of computational time versus signal dimension under fixed K. It is apparent that the computational time of sFFT-DT is not related to N . and a m = 3. Recovery accuracy was measured in terms of SNR, defined in Eq. (12). The experimental results show that sFFT-DT can approximate the first K-largest freq. grids accurately when SN R ori is large enough. For low SN R ori with large N K , the performance degrades.
In sum, compared with [1] [2], the proposed sFFT-DT for generally K-sparse signals is the first algorithm with the reasonable Big-O constants and is verified to be faster than FFTW.

IV. CONCLUSIONS AND FUTURE WORK
We have presented new sparse Fast Fourier Transform methods based on downsampling in the time domain (sFFT-DT) for both exactly k-sparse and generally K-sparse signals in this paper. The accurate computational cost and theoretical performance lower bound of sFFT-DT are proven for exactly K-sparse signals. We also derive the Big-O constants of computational complexity of sFFT-DT and show that they are smaller than those of MIT's methods [1][2] [6]. In addition, sFFT-DT is more hardware-friendly, compared with other algorithms, since all operations of sFFT-DT are linear and involved in an analytical solution. On the other hand, previous works, such as [1][2] [6], are based on the assumption that sparsity K is known in advance. To address this issue, we proposed a simple solution to estimate K and relax this impractical assumption. We showed that the extra cost for deciding K is the same as that required for sFFT-DT with known K. Moreover, we extended sFFT-DT to generally K-sparse signals in this paper. To solve the interference from insignificant freq. grids in aliasing, we presented a new sparse signal recovery algorithm to refine the solution to the MPP.
Overall, theoretical complexity analyses and simulation results demonstrate that our sFFT-DT outperforms the state-of-the-art.
In future work, we want to extend sFFT-DT to a real application, such as GPS [22]. In this case, the sparsity constraint is 1 but X is interfered with massive noises. The most important thing is to identify the location with the maximum magnitude.
A future direction needing further attention is to combine the mechanism of location estimation in our algorithm and sFFT, as shown in Fig. 8. Fig. 8(a) shows an extremely sparse signal X with only one non-zero entry (K = 1). The red region of Fig.   8(b) is the estimated locations using sFFT (detailed descriptions are ignored due to limited space). The red spikes in Fig. 8(c) are the possible d locations derived from our method. By combining both algorithms, the possible location is reduced to the red circle, as shown in Fig. 8(a). Such a new mechanism is more efficient than the one proposed in [22].