Next Article in Journal
Position Control for Soft Actuators, Next Steps toward Inherently Safe Interaction
Previous Article in Journal
Temperature-Sensitivity of Two Microwave HEMT Devices: AlGaAs/GaAs vs. AlGaN/GaN Heterostructures
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On Performance of Sparse Fast Fourier Transform Algorithms Using the Aliasing Filter

School of Mechanical and Electrical Engineering and Automation, Shanghai University, Shanghai 200072, China
*
Author to whom correspondence should be addressed.
Submission received: 27 March 2021 / Revised: 1 May 2021 / Accepted: 4 May 2021 / Published: 9 May 2021
(This article belongs to the Section Circuit and Signal Processing)

Abstract

:
Computing the sparse fast Fourier transform (sFFT) has emerged as a critical topic for a long time because of its high efficiency and wide practicability. More than twenty different sFFT algorithms compute discrete Fourier transform (DFT) by their unique methods so far. In order to use them properly, the urgent topic of great concern is how to analyze and evaluate the performance of these algorithms in theory and practice. This paper mainly discusses the technology and performance of sFFT algorithms using the aliasing filter. In the first part, the paper introduces the three frameworks: the one-shot framework based on the compressed sensing (CS) solver, the peeling framework based on the bipartite graph and the iterative framework based on the binary tree search. Then, we obtain the conclusion of the performance of six corresponding algorithms: the sFFT-DT1.0, sFFT-DT2.0, sFFT-DT3.0, FFAST, R-FFAST, and DSFFT algorithms in theory. In the second part, we make two categories of experiments for computing the signals of different SNRs, different lengths, and different sparsities by a standard testing platform and record the run time, the percentage of the signal sampled, and the L 0 , L 1 , and L 2 errors both in the exactly sparse case and the general sparse case. The results of these performance analyses are our guide to optimize these algorithms and use them selectively.

1. Introduction

The widely popular algorithm to compute DFTs is the famous fast Fourier transform (FFT) [1] invented by Cooley and Tukey, which can reduce the computational complexity of discrete Fourier transform significantly from O( N 2 ) to O( N log N ). The demand for big data computing and low sampling ratios motivates the requirement for the new algorithms to replace previous FFT algorithms that can compute DFT in sublinear time and only use a subset of the input data. The new sFFT algorithms take advantage of the signal’s inherent characteristic that a large number of signals are sparse in the frequency domain, and only K ( K < < N ) frequencies are non-zeros or are significantly large. Under this assumption, people can reconstruct the spectrum with high accuracy by using only K most significant frequencies. Because of its excellent performance and generally satisfied assumptions, the technology of sFFT was named one of the ten breakthrough technologies in MIT Technology Review in 2012.
Through the first stage in the sFFT frequency bucketization, N frequency coefficients are hashed into B buckets, and the length of one bucket is denoted by L. The effect of the Dirichlet kernel filter is to make the convoluted signal into a rectangular window in the time domain; it can be equivalent to the signal multiplied by a Dirichlet kernel window of size L ( L < < N ) in the frequency domain. The effect of the flat window filter is to make the signal multiply a mixed window in the time domain; it can be equivalent to the signal convoluted by a flat window of size L ( L < < N ) in the frequency domain. The effect of the aliasing filter is to make the signal multiply a comb window in the time domain; it can be equivalent to the signal convoluted by a comb window of size B ( B K ) in the frequency domain. After bucketization, the algorithm only needs to focus on the non-empty buckets and locate the positions and estimated values of the large frequency coefficients in those buckets in what we call the identifying frequencies or the spectrum reconstruction. Concerns regarding the performance of these sFFT algorithms include runtime complexity, sampling complexity, and robustness performance. This paper will provide complete answers in theory and practice.
The first sFFT algorithm, called the Ann Arbor fast Fourier transform (AAFFT) algorithm, using the Dirichlet kernel filter, is a randomized algorithm with runtime and sampling complexity O ( K 2 poly ( log N ) ) . The performance of the AAFFT0.5 [2] algorithm was later improved to O ( K poly(log N ) ) in the AAFFT0.9 [3,4] algorithm through the use of unequally-spaced FFTs and the use of a binary search technique for spectrum reconstruction.
The sFFT algorithms using the flat window filter can compute the exact K-sparse signal with runtime complexity O ( K log N ) and general K-sparse signal with runtime complexity O ( K log N log ( N / K ) ) . In the one-shot framework, the sFFT1.0 [5] and sFFT2.0 [5] algorithms can locate and estimate the K largest coefficients in one shot. In the iterative framework, the sFFT3.0 [6] and sFFT4.0 [6] algorithm can locate the position by using only 2B or more samples of the filtered signal inspired by the frequency offset estimation both in the exactly sparse case and general sparse case. Later, a new robust algorithm, so-called the Matrix Pencil FFT (MPFFT) algorithm, was proposed in [7] based on the sFFT3.0 algorithm. The paper [8] summarizes the two frameworks and five reconstruction methods of these five corresponding algorithms both in theory and practice.
There are three frameworks for sFFT algorithms using the aliasing filter. The algorithms of the one-shot framework are the so-called sFFT by downsampling in the time domain (sFFT-DT)1.0 [9], sFFT-DT2.0 [10], and sFFT-DT3.0 algorithms. The algorithms of the peeling framework are the so-called fast Fourier aliasing-based sparse transform (FFAST) [11,12] and R-FFAST (Robust FFAST) [13,14] algorithms. The algorithm of the iterative framework is the so-called deterministic sparse FFT (DSFFT) [15] algorithm. This paper mainly discusses the technology and performance of these six algorithms, and all the details will be mentioned later in the paper.
Under the assumption of arbitrary sampling (while utilizing only a fraction of the FFT’s required samples), the Gopher fast Fourier transform (GFFT) [16,17] algorithm and the Christlieb Lawlor Wang sparse Fourier transform (CLW-SFT) [18,19] algorithm can compute the exact K-sparse signal with runtime complexity O ( K log K ) . They are aliasing-based search deterministic algorithms guided by the Chinese remainder theorem (CRT). The DMSFT [20] (generated from GFFT) algorithm and CLW-DSFT [20] (generated from CLW-SFT) algorithm can compute the general K-sparse signal with runtime complexity O ( K 2 log K ) . They use the multiscale error-correcting method to cope with high-level noise.
The paper [21] summarizes a three-step approach in the stage of spectrum reconstruction and provides a standard testing platform to evaluate different sFFT algorithms. There have also been some researches that tried to conquer the sFFT problem from other aspects: complexity [22,23], performance [24,25], software [26,27], hardware [28], higher dimensions [29,30], implementation [31,32], and special setting [33,34] perspectives.
The identification of different sFFT algorithms can be seen through a brief analysis as above. The algorithms using the Dirichlet kernel filter are not efficient because it only bins some frequency coefficients into one bucket one time. The algorithms using the flat filter are probabilistic algorithms with spectrum leakage [21]. In comparison to them, the algorithms using the aliasing filter are very convenient and have no spectrum leakage, whether N is a product of some co-prime numbers or N is a power of two. This type of algorithm is suitable for the exact sparse case and the general sparse case, but it is not easy to solve the worst case because there may be many frequency coefficients in the same bucket accidentally for the reason that the scaling operation is of no use to the filtered signal. As for the application of the algorithms, there are many practical examples, such as the R-FFAST algorithm using in the magnetic resonance imaging (MRI) application. The paper [14] provides empirical experiments on a real MRI dataset, to demonstrate the feasibility of using R-FFAST sampling with LASSO reconstruction. In the experiments, the peak signal to noise ratio (PSNR) with the FFAST sampling was slightly lower than that with Poisson-disk sampling, and visually the images appeared very similar. s The paper is structured as follows. Section 2 provides a brief overview of the basic sFFT technique using the aliasing filter. Section 3 introduces and analyzes three frameworks and six corresponding spectrum reconstruction methods. In Section 4, we analyze and evaluate the performance of six corresponding algorithms in theory. In the one-shot framework, the sFFT-DT1.0, sFFT-DT2.0 and sFFT-DT3.0 algorithms use the CS solver with the help of the moment preserving problem (MPP) method. In the peeling framework, the FFAST and R-FFAST algorithms use the bipartite graph with the help of the packet erasure channel method. In the iterative framework, the DSFFT algorithm uses the binary tree search with the help of the statistical characteristics. In Section 5, we perform two categories of comparison experiments. The first kind of experiment is to compare the algorithms with each other. The second is to compare them with other algorithms. The analysis of the experiment results satisfies the inferences obtained in theory.

2. Preliminaries

In this section, we initially present some notations and basic definitions of sFFTs.

2.1. Notation

The N-th root of unify is denoted by ω N = e 2 π i / N . The DFT matrix of size N is denoted by F N C N × N as follows:
F N [ j , k ] = 1 N ω N j k
The DFT of a vector x C N (consider a signal of size N) is a vector x ^ C N defined as follows:
x ^ = F N x x ^ i = 1 N j = 0 N 1 x j ω N i j
In the exactly sparse case, spectrum x ^ is exactly K-sparse if it has exactly K non-zero frequency coefficients while the remaining N K coefficients are zero. In the general sparse case, spectrum x ^ is general K-sparse if it has K significant frequency coefficients while the remaining N K coefficients are approximately equal to zero. The goal of sFFT is to recover a K-sparse approximation x ^ by locating K frequency positions f 0 , , f K 1 and estimating K largest frequency coefficients x ^ f 0 , , x ^ f K 1 .

2.2. Frequency Bucketization

The first stage of sFFT is encoding by frequency bucketization. The process of frequency bucketization using the aliasing filter is achieved through shift operation and subsampling operation.
The first technique is the use of shift operation. The offset parameter is denoted by τ R . The shift operation representing the original signal multiplied by matrix S τ . S τ R N × N is defined as follows:
S τ [ j , k ] = 1 , j τ k ( mod N ) 0 , o . w .
If a vector x C N , x = S τ x , x ^ = F N S τ x , such that:
x i = x ( i τ ) x i + τ = x i x ^ i = x ^ i ω τ i
The second technique is the use of the aliasing filter. The signal in the time domain is subsampled such that the corresponding spectrum in the frequency domain is aliased. It also means frequency bucketization. The subsampling factor is denoted by L Z + , representing how many frequencies are aliasing in one bucket. The subsampling number is denoted by B Z + , representing the number of buckets ( B = N / L ) . The subsampling operation representing the original signal multiplied by matrix D L . D L R B × N is defined as follows:
D L [ j , k ] = 1 , k = j L 0 , o . w .
Let vector y ^ B , τ C B be the filtered spectrum obtained by shift operation and subsampling operation. If y ^ B , τ = F B D L S τ x , we obtain Equation (6) in bucket i.
y ^ B , τ [ i ] = j = 0 L 1 x ^ j B + i ω τ ( j B + i )
As we see above, frequency bucketization includes three steps: shift operation ( x = S τ x , it costs 0 runtime), subsampling operation ( y B , τ = D L S τ x , it costs B samples), Fourier transform ( y ^ B , τ = F B D L S τ x , it costs B log B runtime). Hence one round of total frequency bucketization costs B log B runtime and B samples.

3. Techniques

In this section, we introduce an overview of the techniques and frameworks that we will use in the sFFT algorithms based on the aliasing filter.
As mentioned above, frequency bucketization can decrease runtime and sampling complexity because all operations are calculated in B dimensions ( B = O ( K ) , B < < N ) . After frequency bucketization, it needs spectrum reconstruction performed by identifying frequencies that are isolated in their buckets. The aliasing filter may lead to frequency aliasing where more than one significant frequency are aliasing in one bucket. This increases the difficulty in recovery because finding frequency position and estimating frequency values becomes indistinguishable in terms of their aliasing characteristics. There are three frameworks to overcome this problem, namely the one-shot framework, the peeling framework, the iterative framework.

3.1. One-Shot Framework Based on the CS Solver

Firstly we introduce the one-shot framework, which can recover all K significant frequencies in one shot. The block diagram of the one-shot framework is shown in Figure 1. The concepts, technology, and framework involved in this section were proposed by the Taiwan university in the paper [9,10].
The first stage of sFFT is encoding by frequency bucketization. Suppose there are at most an a m number of significant frequencies aliasing in every bucket, after running 3 a m times for the set τ = { τ 0 = a m , τ 1 = a m + 1 , , τ a m = 0 , , τ 3 a m 1 = 2 a m 1 } , calculate y ^ B , τ = F B D L S τ x , representing the filtered spectrum by encoding. Suppose that in bucket i, the number of significant frequencies is denoted by a; there is a high probability that a a m , we obtain simplified Equation (9) from Equations (7) and (8). In Equations (8) and (9), p j = x ^ f j respecting effective frequency values for | p 0 | | p 1 | | p a 1 | other values of p j , z j = ω f j respecting effective frequency position for f j i , i + L , , i + ( L 1 ) B ) , m k = y ^ B , k [ i ] respecting filtered spectrum in bucket i for k a m , , 2 a m 1 . In most cases, a = 0 respecting sparsity. In a small number of cases, a = 1 respecting only one significant frequency in the bucket. Only in very few cases, a > = 2 respecting frequencies aliasing in the bucket. It is very unlikely that a > a m . In the exactly sparse case, the approximately equal sign becomes the equal sign in Equations (7)–(9), (11), (14), and (18).
The spectrum reconstruction problem in bucket i is equivalent to obtaining unknown variables including a , z 0 , , z a 1 , p 0 , , p a 1 , as we have known variables including m a , , m 2 a 1 . The aliasing problem is reformulated as a moment preserving problem (MPP). The MPP problem formulated by Bose–Chaudhuri–Hocquenghem (BCH) codes [35] can be divided into three subproblems: how to obtain a, how to obtain z j ’s, and how to obtain p j ’s in every bucket. Below, we will solve these three subproblems step by step.
y ^ B , τ [ i ] = j = 0 L 1 x ^ j B + i ω τ ( j B + i ) j = 0 a m 1 x ^ f j ω τ f j
m k j = 0 a m 1 p j z j k j = 0 a 1 p j z j k
z 0 0 z 1 0 z a 1 0 z 0 1 z 1 1 z a 1 1 z 0 2 a 1 z 1 2 a 1 z a 1 2 a 1 p 0 p 1 p a 1 m 0 m 1 m 2 a 1
Step 1: Obtain the number of significant frequencies.
Solution: Suppose a a m ; this means m k is composed of at most an a m number of Prony component p j ’s. Let vector z ^ j C a m × 1 defined as z ^ j = [ z j 0 , z j 1 , , z j a m 1 ] T and Matrix M a m C a m × a m defined as Equation (10), then the relationship between M a m , z ^ j and z ^ j T satisfies Theorem 1.
M a m = m 0 m 1 m a m 1 m 1 m 2 m a m m a m 1 m a m m 2 a m 1 a m × a m
Theorem 1. 
M a m j = 0 a m 1 p j z ^ j z ^ j T
Proof of Theorem 1. 
p 0 z ^ 0 z 0 0 + p 1 z ^ 1 z 1 0 + + p a m 1 z ^ a m 1 z a m 1 0 [ m 0 , m 1 , , m a m 1 ] T p 0 z ^ 0 z 0 1 + p 1 z ^ 1 z 1 1 + + p a m 1 z ^ a m 1 z a m 1 1 [ m 1 , m 2 , , m a m ] T p 0 z ^ 0 z 0 a m 1 + p 1 z ^ 1 z 1 a m 1 + + p a m 1 z ^ a m 1 z a m 1 a m 1 [ m a m 1 , m a m , , m 2 a m ] T
Based on the properties as mentioned above, we obtain Equation (11). □
Equation (11) is similar to the symmetric singular value decomposition (SSVD). Nevertheless, there are some differences. (1) p j ’s are complex but not real. (2) The columns of [ z ^ 0 , , z ^ a m 1 ] are not mutually orthogonal normalized vectors. It is easy to perform a transformation that M a m j = 0 a m 1 p j z ^ j z ^ j T , where p j ’s are real and the absolute value of p j is directly proportional to the p j , and the columns of [ z ^ 0 , , z ^ a m 1 ] are mutually orthogonal normalized vectors. The paper [36] proved that for the symmetric matrix, the Σ obtained from the SVD is equal to the Σ obtained from the SSVD. For example, the SVD of 0 1 1 0 is 0 1 1 0 = 0 1 1 0 1 0 0 1 1 0 0 1 and the SSVD of 0 1 1 0 is 0 1 1 0 = 1 2 1 i 1 i 1 0 0 1 1 2 1 1 i i ; the Σ values gained via these two methods are the same. After knowing this, we can compute the SVD of M a m and obtain a m singular values, then perform the principal component analysis (PCA), σ 1 σ 2 σ a m . In these a m number of singular values σ j ’s, the amount of large singular values indicates the amount of efficient Prony components p j ’s; it also indicates how many significant frequencies are in bucket i.
Step 2: Obtain effective frequency position f j ’s in bucket i.
Solution: Let the orthogonal polynomial formula P ( z ) be defined as Equation (12) and P ( z ) 0 . Let Matrix M a C a × a be defined as Equation (13). Let vector C be defined as C = [ c 0 , c 1 , , c a 1 ] T . Let vector M s be defined as M s = [ m a , m a + 1 , , m 2 a 1 ] T . The moments’ formula satisfies Theorem 2.
P ( z ) = z a + c a 1 z a 1 + + c 1 z + c 0
M a = m 0 m 1 m a 1 m 1 m 2 m a m a 1 m a m m 2 a 2 a × a
Theorem 2. 
M a C M s
m 0 m 1 m a 1 m 1 m 2 m a m a 1 m a m m 2 a 2 a × a c 0 c 1 c a 1 a × 1 m a m a + 1 m 2 a 1 a × 1
Proof of Theorem 2. 
c 0 m 0 ( p 0 z 0 0 + p a 1 z a 1 0 ) c 0 c a 1 m a 1 ( p 0 z 0 a 1 + p a 1 z a 1 a 1 ) c a 1 c 0 m 0 + + c a 1 m a 1 p 0 ( c 0 z 0 0 + c a 1 z 0 a 1 ) + + p a 1 ( c 0 z a 1 0 + c a 1 z a 1 a 1 ) ( p 0 z 0 a ) + + ( p a 1 z a 1 a ) m a
The first element of M s has been proven and other elements of M s can also be proven. □
For the convenience of matrix calculation, for a matrix, the superscript “T” denotes the transpose, the superscript “+” denotes the Moore–Penrose inverse or pseudoinverse, the superscript “*” denotes the adjoint matrix, and the superscript “−1” denotes the inverse. Through Theorem 2, we can obtain C ( M a 1 ) M s . After gaining C, there are three ways to obtain z j ’s through Equation (12). The first approach is the polynomial method, the second approach is the enumeration method, and the last approach is the matrix pencil method. After knowing z j ’s, we can obtain the approximate positions f j ’s through z j = ω f j .
(Method 1) Polynomial method: In the exactly sparse case, the a number of roots of P ( z ) = 0 is the solution of z j ’s. For example, if a = 1 , through P ( z ) = z + c 0 = 0 , then z 0 = c 0 . If a = 2 , through P ( z ) = z 2 + c 1 z + c 0 = 0 , then z 0 = ( c 1 ( c 1 2 4 c 0 ) 0.5 ) / 2 , z 1 = ( c 1 + ( c 1 2 4 c 0 ) 0.5 ) / 2 .
(Method 2) Enumeration method: For f j i , i + L , , i + ( L 1 ) B ) , try every possible position z j = ω f j in Equation (12), the first a number of the smallest | P ( z ) | is the solution of z j ’s. L attempts are needed in one solution.
(Method 3) Matrix pencil method: The method was proposed in paper [37]. The matrix pencil method, like the Prony method, is a standard technique for mode frequency identification for computing the maximum likelihood signal under Gaussian noise and evenly spaced samples. For our problem, let the Toeplitz matrix Y be defined as Equation (15), let Y 1 be Y with the rightmost column removed and be defined as Equation (16), and let Y 2 be Y with the leftmost column removed and be defined as Equation (17). The set of generalized eigenvalues of Y 2 λ Y 1 satisfies Theorem 3.
Y = m 0 m 1 m a m 1 m 0 m a + 1 m a m a 1 m 0 ( a + 1 ) × ( a + 1 )
Y 1 = m 0 m 1 m a + 1 m 1 m 0 m a + 2 m a m a 1 m 1 ( a + 1 ) × ( a )
Y 2 = m 1 m 2 m a m 0 m 1 m a + 1 m a 1 m a 2 m 0 ( a + 1 ) × ( a )
Theorem 3. 
The set of generalized eigenvalues of Y 2 λ Y 1 are the z j ’s we seek.
Proof of Theorem 3. 
Let the diagonal matrix C C a × a be defined as C = diag ( p j ) . Let the Vandermonde martix U a + 1 C ( a + 1 ) × a be defined as follows: U a + 1 = 1 1 1 z 0 z 1 z a 1 z 0 a z 1 a z a 1 a ( a + 1 ) × a . Y , Y 1 , Y 2 has a Vandermonde decomposition, we can obtain Y = 1 a + 1 U a + 1 C U a + 1 * , Y 1 = 1 a + 1 U a + 1 C U a * , Y 2 = 1 a + 1 U a + 1 C ( diag ( z j ) * ) U a * . For example, if a = 1, Y = m 0 m 1 m 1 m 0 = p 0 p 0 z 0 1 p 0 z 0 1 p 0 = 1 z 0 p 0 1 z 0 1 , and if a = 2, Y = m 0 m 1 m 2 m 1 m 0 m 1 m 2 m 1 m 0 = p 0 + p 1 p 0 z 0 1 + p 1 z 1 1 p 0 z 0 2 + p 1 z 1 2 p 0 z 0 1 + p 1 z 1 1 p 0 + p 1 p 0 z 0 1 + p 1 z 1 1 p 0 z 0 2 + p 1 z 1 2 p 0 z 0 1 + p 1 z 1 1 p 0 + p 1 = 1 1 z 0 z 1 z 0 2 z 1 2 p 0 p 1 1 z 0 1 z 0 2 1 z 1 1 z 1 2 , Y 1 = m 0 m 1 m 1 m 0 m 2 m 1 = 1 1 z 0 z 1 z 0 2 z 1 2 p 0 p 1 1 z 0 1 1 z 1 1 , Y 2 = m 1 m 2 m 0 m 1 m 1 m 0 = 1 1 z 0 z 1 z 0 2 z 1 2 p 0 p 1 z 0 1 z 1 1 1 z 0 1 1 z 1 1 . Using the Vandermonde decomposition, and we can obtain Y 2 λ Y 1 = 1 a + 1 U a + 1 C ( diag ( z j ) * λ I ) U a * , so the Theorem 3 can be proven. □
If the rank ( Y ) = a , the set of generalized eigenvalues of Y 2 λ Y 1 is equal to the set of nonzero 2 ordinary eigenvalues of ( Y 1 + ) Y 2 . It is most likely that for the rank ( Y ) < a , it is necessary to compute the SVD of the Y , Y = V ˜ Σ V * , and then we can use the matrix pencil method to deal with the matrix V afterword. For details, please refer to paper [7,37].
Step 3: Obtain effective frequency values p j ’s in bucket i.
Solution: In order to use the CS method, we need several random samplings. Thus, for a P number of random numbers r 1 , , r P , we calculate y ^ B , τ = F B D L S τ x for the set τ = τ 0 = r 1 , τ 1 = r 2 , , τ P 1 = r P in a P times’ round. Suppose that in bucket i, the number of significant frequencies a and approximate effective frequency position z j ’s have been known by step 1 and step 2; we can obtain Equation (18). (There may be errors for z j ’s obtained by step 2 because of the interference of another L a number of Prony components).
m 1 m P P × 1 = z 0 r 1 z L 1 r 1 z 0 r P z L 1 r P P × L p 0 p L 1 L × 1 z 0 r 1 z a 1 r 1 z 0 r P z a 1 r P P × a p 0 p a 1 a × 1
Equation (18) is very likely similar to the CS formula. The model of CS is formulated as y Φ S , where S is a sparse signal, Φ is a sensing matrix, and y is the measurements. In Equation (18), y is a vector of P × 1 by P measurements, Φ is a matrix of P × L , and S is a vector of L × 1 that is a-sparse. It should be noted that Φ must satisfy either the restricted isometry property (RIP) or mutual incoherence property (MIP) for a successful recovery with high probability. It has been known that the Gaussian random matrix and partial Fourier matrix are good candidates to be Φ , so the Φ of Equation (18) meets the criteria. Furthermore, the number of measurements P one collects should be more than a log 10 L , so that these measurements will be sufficient to recover signal x.
In order to obtain p j ’s using the CS solver, we use the subspace pursuit method. The process is as follows: (1) Through the positions f j ’s gained by step 2, obtain the possible value of z j ’s as follows: z j = ω f j , z j = ω f j + B , z j = ω f j B , then obtain 3a vectors as follows: { z 0 r 1 , , z 0 r P } T , { z 0 r 1 , , z 0 r P } T , { z 0 r 1 , , z 0 r P } T , { z a 1 r 1 , , z a 1 r P } T . An (over-complete) dictionary can be characterized by a matrix D , and it contains the 3 a vectors listed above. (One wishes one-third vectors of them form a basis). Each (column) vector in a dictionary is called an atom. (2) From 3 a atoms of the dictionary matrix D , find a number of atoms that best match the measurements’ residual error. Select these a number of atoms to construct a new sensing matrix Φ ^ . (3) Obtain S ^ by the support of S ^ = argmin y Φ ^ S ^ 2 through the least square method (LSM). (4) If the residual error y Φ ^ S ^ 2 meets the requirement, or the number of iterations reaches the assumption, or the residual error becomes larger, the iteration will be quit; otherwise continue to step 2. After computing, we obtain the final sensing matrix Φ ^ of size P × a and sparse signal S ^ of size a just in the right part of Equation (18). Thus we obtain effective frequency positions z j ’s and effective frequency values p j ’s in bucket i.

3.2. Peeling Framework Based on the Bipartite Graph

Secondly, we introduce the peeling framework, which can recover all K significant frequencies layer by layer. The block diagram of the peeling framework is shown in Figure 2. The concepts, technology, and framework involved in this section were proposed by Berkeley university in the paper [11,12,13,14].
In the peeling framework, we require the size N of signal x to be a product of a few (typically three or more) co-prime numbers p i ’s. For example N = p 1 p 2 p 3 , where p 1 , p 2 , p 3 are co-prime numbers. With this precondition, we determine L i ’s and B i ’s gained from p i ’s in each bucketization cycle. In every cycle, we use the same set τ = τ 0 = r 1 , τ 1 = r 2 , , τ R 1 = r R inspired by different spectrum reconstruction methods introduced later to calculate y ^ B , τ = F B D L S τ x in R times’ rounds (this also means there are R delay chains for one bucket). Suppose there are d cycles, after d cycles, stage 1 encoding by bucketization is completed. In order to solve the aliasing problems, the filtered spectrum in bucket i of the no. j’ cycle is denoted by vector y B j , i C R × 1 as follows:
y B j , i = y ^ B j , τ 0 [ i ] y ^ B j , τ R 1 [ i ] = k = 0 L j 1 x ^ ( k B j + i ) ω τ 0 ( k B j + i ) k = 0 L j 1 x ^ ( k B j + i ) ω τ R 1 ( k B j + i )
We use a simple example to illustrate the process of encoding by bucketization. Consider a signal x of size (N = 20) that has only five (K = 5) significant coefficients, x ^ 1 , x ^ 3 , x ^ 5 , x ^ 10 , x ^ 13 > > 0 , while the rest of the coefficients are approximately equal to zero. With this precondition, there are two bucketization cycles. In the first cycle, for B 1 = 4 and L 1 = 5 , we obtain four vectors, { y 4 , 0 , y 4 , 1 , y 4 , 2 , y 4 , 3 } , respecting the filtered spectrum in four buckets for the set τ in R rounds. In the second cycle, for B 2 = 5 and L 2 = 4 , we obtain five vectors, { y 5 , 0 , y 5 , 1 , y 5 , 2 , y 5 , 3 , y 5 , 4 } . After the bucketization, we can construct a bipartite graph shown in Figure 3 through Equation (19). In Figure 3, there are N = 20 variable nodes on the left (referring to the 20 coefficients of x ^ ) and N b = B 1 + B 2 = 4 + 5 = 9 parity check nodes on the right (referring to nine buckets in two cycles). The values of the parity check nodes on the right are approximately equal to the complex sum of the values of variable nodes (its left neighbors) through Equation (19). In these check nodes, some have no significant variable node with no left neighbor, which is called a “zero-ton” bucket (three blue-colored check nodes). Some have exactly one significant variable node, as one left neighbor, which is called a “single-ton” bucket (three green-colored check nodes). Others have more than one significant variable nodes, as more than one left neighbors, which is called a “multi-ton” bucket (three red-colored check nodes).
After bucketization, the subsequent problem is how to recover the spectrum from all buckets gained from several cycles. Through the identification of vector y B j , i , we can determine the characteristics of bucket B j [ i ] . If the bucket is a “zero-ton” bucket, then x ^ ( k B j + i ) = 0 for k [ 0 , L j 1 ] , so the problem of frequency recovery in this bucket can be solved. If the bucket is a “single-ton” bucket, suppose the one effective frequency position f 0 and the one effective frequency value x ^ f 0 can be obtained by the afterward methods, then x ^ ( k B j + i ) = 0 for ( k [ 0 , L j 1 ] ) ( k B j + i f 0 ) x ^ f 0 for k B j + i = f 0 , and thus the frequency recovery in this bucket can be solved. If the bucket is a “multi-ton” bucket, it is necessary to separate the “multi-ton” bucket into the “single-ton” bucket to realize the decoding of the “multi-ton” bucket. For example, after the first peeling, in the bucket i of the no. j cycle, suppose x ^ f 0 , x ^ f n 1 1 have been obtained via another “single-ton” bucket, then the vector y B j , i respecting the original bucket changes to y B j , i , where y B j , i = y B j , i poly ( x ^ f 0 ) poly ( x ^ f n 1 1 ) respecting the remaining frequencies in the bucket. Through the identification of vector y B j , i , we can analyze the remaining elements in this bucket; if it is a “zero-ton” bucket or a “single-ton” bucket, we can stop peeling and obtain all frequencies in this bucket. If not, continue the second peeling, and suppose x ^ f n 1 , x ^ f n 2 1 can be obtained by another “single-ton” bucket through new peeling. We can identify y B j , i to continue to analyze the bucket, where y B j , i = y B j , i poly ( x ^ f n 1 ) poly ( x ^ f n 2 1 ) . After q times’ peeling, the problem of frequency recovery in the “multi-ton” bucket can be solved as follows:
x ^ ( k B j + i ) = 0 for ( k [ 0 , L j 1 ] ) ( k B j + i f 0 ) ( k B j + i f n q 1 ) x ^ f 0 for k B j + i = f 0 x ^ f n q 1 for k B j + i = f n q 1
If the frequency recovery in all buckets can be solved, we can finish the spectrum reconstruction. The peeling-decoder successfully recovers all of the frequencies with high probability under the three following assumptions: (1) “Zero-ton”, “single-ton” and “multi-ton” buckets can be identified correctly. (2) If the bucket is a “single-ton” bucket, the decoder can locate the effective frequency position f 0 and estimate the effective frequency value x ^ f 0 correctly. (3) It is sufficient to cope with “multi-ton” buckets via the peeling platform.
Subproblem 1: How to identify vector y B j , i to distinguish bucket B j [ i ] ?
Solution: In the exactly sparse case, if y B j , i 2 = 0 , the bucket is a “zero-ton” bucket. If the bucket is not a “zero-ton” bucket, make the second judgment. The way to make the second judgment about whether the bucket is a ”single-ton” bucket or not is to judge y B j , i poly ( x ^ f 0 ) 2 = 0 or 0 , where x ^ f 0 is gained from the solution of subproblem 2. If the bucket B j [ i ] is not a “zero-ton” bucket nor a “single-ton” bucket, it is a “multi-ton” bucket. In the general sparse case, the two equations are translated to y B j , i 2 T 0 and y B j , i poly ( x ^ f 0 ) 2 T 1 , where T 0 and T 1 are the identification thresholds. It costs O ( R ) runtime to identify vector y B j , i by knowing x ^ f 0 in advance. After d cycles, O ( R ( B 1 + + B d ) ) runtime is needed for the identification in the first peeling.
Subproblem 2: Suppose the target is a “single-ton” bucket, how to recover the one frequency in this bucket?
Solution: In the exactly sparse case, we use R = 3 and the set τ = τ 0 = 0 , τ 1 = 1 , τ 2 = 2 to calculate y ^ B , τ = F B D L S τ x in three rounds. Suppose the bucket B j [ i ] is a ”single-ton” bucket, then three elements of the vector y B j , i are y ^ B j , 0 [ i ] = x ^ f 0 ω 0 · f 0 , y ^ B j , 1 [ i ] = x ^ f 0 ω 1 · f 0 and y ^ B j , 2 [ i ] = x ^ f 0 ω 2 · f 0 . Thus, only if y ^ B j , 0 [ i ] y ^ B j , 1 [ i ] = y ^ B j , 1 [ i ] y ^ B j , 2 [ i ] , it is a ”single-ton” bucket. If it is a ”single-ton” bucket, the position f 0 can be obtained by f 0 = y ^ B j , 1 [ i ] y ^ B j , 0 [ i ] · ( N 2 π ) , and the value x ^ f 0 can be obtained by x ^ f 0 = y ^ B j , 0 [ i ] . In all, it costs three samples and four runtimes to recover the one significant frequency in one bucket.
In the general sparse case, the frequency recovery method is the optimal whitening filter coefficient of the minimum mean squared error (MMSE) estimator and sinusoidal structured bin-measurement matrix for the speedy recovery. At first, the set of the one candidate f 0 is f 0 i , i + L , , i + ( L 1 ) B ) , and there are L possible choices. In the first iteration of the binary search, we choose a random number r 0 , then calculate y ^ B j , r 0 , y ^ B j , r 0 + 1 , , y ^ B j , r 0 + m 1 in m rounds. In fact, we obtain y ^ B j , r 0 [ i ] x ^ f 0 ω ( r 0 ) · f 0 , y ^ B j , r 0 + 1 [ i ] x ^ f 0 ω ( r 0 + 1 ) · f 0 , , y ^ B j , r 0 + m 1 [ i ] x ^ f 0 ω ( r 0 + m 1 ) · f 0 . Δ t respecting the phase difference is defined as Δ t = y ^ B j , r 0 + t + 1 [ i ] y ^ B j , r 0 + t [ i ] = ω f 0 + U t + 1 U t , where U t relates to the real error in no. t’ round. In paper [38], we can see the maximum likelihood estimate (MLE) ω ^ f 0 of ω f 0 is calculated by Equation (21), where w t is defined as Equation (22). After obtaining ω ^ f 0 , make a judgment of binary search; if ω ^ f 0 [ 0 , π ] , there is a new restriction that f 0 [ 0 , N / 2 1 ] , otherwise the restriction is f 0 [ N / 2 , N 1 ] the complement set of set [ 0 , N / 2 1 ] . The next iteration is very similar; we choose a random number r 1 , then calculate y ^ B j , r 1 [ i ] , y ^ B j , r 1 + 2 [ i ] , , y ^ B j , r 1 + 2 ( m 1 ) [ i ] in bucket i. After obtaining ω ^ 2 f 0 through Equations (21) and (22), we make a judgment of binary-search; if ω ^ 2 f 0 [ 0 , π ] , there is a new restriction that f 0 [ 0 , N / 4 1 ] [ N / 2 , 3 N / 4 1 ] , otherwise the restriction is the complement set of the previous set. After C iterations, in the advantage of restrictions, we can locate the only position f 0 from the original set i , i + B , , i + ( L 1 ) B ) . From the paper [13], if the number of iterations C = O ( log N ) and the number of rounds per iteration m = O ( log 1 / 3 N ) , the singleton-estimator algorithm can correctly identify the unknown frequency f 0 with a high probability. If the signal-to-noise ratio (SNR) is low, the m and C must be increased. After knowing f 0 , we can obtain x ^ f 0 by applying the LSM. In all, to recover the one approximately significant frequency in one bucket, it needs C m = O ( log 4 / 3 N ) samples and 3 C m = O ( log 4 / 3 N ) runtime. The runtime includes C m runtime to calculate all Δ t ’s and 2 C m runtime to calculate all ω ^ 2 j f 0 ’s. x ^ f 0 can be used to judge whether the bucket is a “single-ton” bucket or not through the discriminant y B j , i poly ( x ^ f 0 ) T 1 .
ω ^ f 0 = t = 0 m 2 w t Δ t
w t = 3 m 2 ( m 2 1 ) ( 1 4 ( 2 t m + 2 2 m ) 2 )
Subproblem 3: How to solve the “multi-ton” buckets by the peeling platform?
Solution with method 1: The genie-assisted peeling decoder by the bipartite graph is useful. As shown in Figure 3, the bipartite graph represents the characteristics of the bucketization. The variable nodes on the left represent the characteristics of the frequencies. The parity check nodes on the right represent the characteristics of the buckets. Every efficient variable node connects d different parity check nodes as its neighbors in d cycles; it respects the d edges connected to each efficient variable node. The example of the process of the genie-assisted peeling decoder is shown in Figure 4, and the steps are as follows:
Step 1: We can identify all the indistinguishable buckets. If the bucket is a “zero-ton” bucket, the frequency recovery in this bucket is finished. If the bucket is a “single-ton” bucket, we can obtain the frequency in the bucket through the solution of subproblem 2. In the graph, these operations represent to select and remove all right nodes with degree 0 or degree 1, and moreover, to remove the edges connected to these right nodes and corresponding left nodes.
Step 2: We should remove the contributions of these frequencies gained by step 1 in other buckets. For example, the first peeling in bucket i means we calculate a new vector y B j , i just as y B j , i = y B j , i poly ( x ^ f 0 ) poly ( x ^ f n 1 1 ) , where x ^ f 0 , x ^ f n 1 1 have been gained by step 1. In the graph, these operations represent removing all the other edges connected to the left nodes removed in step 1. When the edges are removed, their contributions are subtracted from their right nodes. In the new graph, the degree of some right nodes decreases as well.
Step 3: If all buckets have been identified, we successfully reconstruct the spectrum x ^ . Otherwise, we turn to step 1 and step 2 to continue identifying and peeling operations. In the graph, it means that if all right nodes are removed, the decoding is finished. Otherwise, turn to step 1 and step 2.
Solution with method 2: The sparse-graph decoder by the packet erasure channel method is more efficient. From Figure 4, we see all processes need 13 (9 + 3 + 1) identifications in three peelings, so it is not very efficient. As we can see, spectrum recovery can transform into a problem of decoding over sparse bipartite graphs using the peeling platform. The problem of decoding over sparse bipartite graphs has been well studied in the coding theory literature. From the coding theory literature, we know that several sparse-graph code constructions are of low-complexity and capacity-achieving by using the erasure channels method. Thus we use the erasure channels method to improve efficiency.
We construct a bipartite graph with K variable nodes on the left and N b check nodes on the right. Each efficient left node x ^ f connects exactly d right nodes y B j , i ’s in d cycles and the set of N b check nodes are assigned to d subsets with the no. i’s subset having B i check nodes. The example of the “balls-and-bins” model defined above is shown in Figure 2. In Figure 3, there are K ( K = 5 ) variable nodes on the left and N b ( N b = 9 ) check nodes on the right and each left node x ^ f connects exactly d ( d = 2 ) right nodes; the set of check nodes is assigned to two subsets ( { y 4 , 0 , y 4 , 1 , y 4 , 2 , y 4 , 3 } , { y 5 , 0 , y 5 , 1 , y 5 , 2 , y 5 , 3 , y 5 , 4 } ). If one left node is selected, a d number of its neighboring right nodes will be determined. If one right node is selected, an L number of its neighboring left nodes will be determined as well. These two corresponding nodes are connected through an edge (in the graph, if the variable nodes on the left are efficient nodes, the edge is a solid line. Otherwise, the edge is a dotted line or is omitted). A directed edge e in the graph is represented as an ordered pair of nodes such as e = { V , C } or e = { C , V } , where V is a variable node and C is a check node. A path in the graph is a directed sequence of directed edges such as e 1 , , e l , where the end node of the previous edge of the path is the start node of the next edge of the path. Depth l is defined as the length of the path, which also indicates the number of directed edges in the path. The induced subgraph N V l is the directed neighborhood of depth l of left node V. It contains all of the edges and nodes on paths e 1 , , e l starting at node V. It is a tree-like graph that can be seen in Figure 5. In Figure 5, it can be seen that subgraph N V l starting at x ^ 2 with depth l is equal to two, three, four, five, and six, and subgraph N V l starting at x ^ 7 with depth l is equal to six. Under these definitions, we obtain the steps of thepacket erasure channel method as follows: Step 1: Take O ( K ) random left nodes as starting points, and draw O ( K ) trees N V 1 from these starting points. As shown in Figure 5, we choose two left nodes x ^ 2 and x ^ 7 as starting points. The endpoints of N V 1 are check nodes; we then identify the characteristics of these check nodes. If the check node is a “zero-ton” bucket, such as y 5 , 2 , this path stops extending. If the check node is a “multi-ton” bucket, continue waiting until its left neighboring node is identified by other paths. If the check node is a “single-ton” bucket, such as y 4 , 2 and y 4 , 3 , continue to connect its only efficient left neighboring node. Then obtain the tree N V 2 through expending these left nodes. Their ( d 1 ) number of new right neighboring nodes will be determined through these left nodes; then, obtain the tree N V 3 by expending these right nodes.
Step p: For each start-point V, we have obtained tree N V 2 p 1 from the last step. The endpoints of N V 2 p 1 are check nodes. We should remove the contributions of some frequencies gained by the previous paths at first, then identify the characteristics of these modified check nodes. For example, before identifying y 5 , 0 , the endpoint of N x ^ 2 3 , we should remove the contributions of x ^ 10 . Furthermore, before identifying y 4 , 1 , the endpoint of N x ^ 2 5 , we should remove the contributions of x ^ 13 and x ^ 5 . If the modified check node is a “zero-ton” bucket, this path stops extending. If the modified check node is a “multi-ton” bucket, continue waiting until its left neighboring node is identified by other paths. If the modified check node is a “single-ton” bucket, continue to connect its only efficient left neighboring node. Then obtain the tree N V 2 p through expending these left nodes. Their ( d 1 ) number of new right neighboring nodes will be determined through these left nodes; then obtain the tree N V 2 p + 1 by expending these right nodes.
If the number of left nodes identified is equal to K, it means the spectrum recovery has been successful. The example is shown in Figure 5. In Figure 5, from the beginning of starting points x ^ 2 and x ^ 7 , we can obtain two trees, N x ^ 2 6 and N x ^ 7 6 (depth = 6), through three expanding by three steps. From the graph, we can obtain all five variable nodes. All processes need six (4 + 4 − 2) identifications, which is far less than the thirteen identifications of method 1.

3.3. Iterative Framework Based on the Binary Tree Search Method

Thirdly, we introduce the iterative framework based on the binary tree search method. The example is shown in Figure 6. The concepts, technology, and framework involved in this section were proposed by Georg-August university in the paper [15].
Consider a signal x ^ sized N = 64 that has only five ( K = 5 ) significant coefficients x ^ 1 , x ^ 6 , x ^ 13 , x ^ 20 , x ^ 59 > > 0 , while the rest of the coefficients are approximately equal to zero. The node of the first layer of the binary tree is y ^ 1 , 0 [ 0 ] = j = 0 N 1 x ^ j , with N aliasing frequencies. The nodes of the second layer of the binary tree are y ^ 2 , 0 [ 0 ] = j = 0 N / 2 1 x ^ 2 j and y ^ 2 , 0 [ 1 ] = j = 0 N / 2 1 x ^ 2 j + 1 , with N / 2 aliasing frequencies. The nodes of the third layer of the binary tree are y ^ 4 , 0 [ 0 ] , y ^ 4 , 0 [ 1 ] , y ^ 4 , 0 [ 2 ] , y ^ 4 , 0 [ 3 ] , with N / 4 aliasing frequencies. The nodes of the no. ( d + 1 ) ’s layer of the binary tree are y ^ 2 d , 0 [ 0 ] , y ^ 2 d , 0 [ 1 ] y ^ 2 d , 0 [ 2 d 1 ] , with N / ( 2 d ) aliasing frequencies.
The insight of the binary tree search is as follows: (1) The frequency of aliasing is gradually dispersed with the expansion of binary tree layers. The number of efficient nodes of the no. ( d + 1 ) ’s layer is defined as m d and the final number will be approximately equal to the sparse number K with the expansion. In Figure 6, m 0 = 1 , m 1 = 2 , m 2 = 4 , m 3 = 5 , and the final number of efficient nodes of no. 4’s layer m 3 is equal to the sparse number K. (2) If the parent node exists, there may be one or two child nodes. If the parent node does not exist, there is no need to continue the binary search for this node. For the example of Figure 6, parent node y ^ 1 , 0 [ 0 ] exists, so at least one of the two child nodes y ^ 2 , 0 [ 0 ] and y ^ 2 , 0 [ 1 ] exists. By the same reason, parent node y ^ 2 , 0 [ 0 ] and y ^ 2 , 0 [ 1 ] exist, so at least one of the two child nodes y ^ 4 , 0 [ 0 ] and y ^ 4 , 0 [ 2 ] exists, and at least one of the two child nodes y ^ 4 , 0 [ 1 ] and y ^ 4 , 0 [ 3 ] exists. On the contrary, parent node y ^ 8 , 0 [ 0 ] does not exist, so its two child nodes y ^ 16 , 0 [ 0 ] and y ^ 16 , 0 [ 8 ] do not exist either. Inspired by these two ideas, the steps of of binary tree search are as follows:
Step 1: Calculate the node y ^ 1 , 0 [ 0 ] of the first layer, obtain m 0 and efficient nodes’ distribution in this layer.
Step 2: According to the last result, calculate the node y ^ 2 , 0 [ 0 ] and y ^ 2 , 0 [ 1 ] of the second layer selectively, and then obtain m 1 and the efficient nodes’ distribution in this layer.
Step d + 1 : According to the last result, calculate the node y ^ 2 d , 0 [ 0 ] , y ^ 2 d , 0 [ 1 ] , …, y ^ 2 d , 0 [ 2 d 1 ] of the no. ( d + 1 ) ’s layer selectively, and then obtain m d and the efficient nodes’ distribution in this layer.
We do not need to start from step 1, we can start from step p ( p = O ( log K ) ) . With the binary tree search, if the m d gained by step d + 1 is approximately equal to the sparse number K, the binary tree search is finished. In Figure 6, m 3 = 5 = K , the binary tree search is finished after the fourth layer by step 4. According to the efficient nodes’ distribution of no. ( d + 1 ) ’s layer, we can solve the frequency recovery problem of each “single-ton” bucket. Finally, the K frequency recovery problem is solved by combining the efficient frequency in each bucket. In Figure 6, y ^ 8 , 0 [ 4 ] , y ^ 8 , 0 [ 6 ] , y ^ 8 , 0 [ 1 ] , y ^ 8 , 0 [ 5 ] , y ^ 8 , 0 [ 3 ] exists; it represents that each of the five buckets has exactly one effective frequency. Through the frequency recovery of each “single-ton” bucket, we obtain x ^ 1 , x ^ 6 , x ^ 13 , x ^ 20 , x ^ 59 individually in every bucket. Finally, we obtain x ^ of five elements. This algorithm involves three subproblems as follows:
Subproblem 1: How to calculate y ^ 2 d , 0 [ 0 ] , y ^ 2 d , 0 [ 1 ] , …, y ^ 2 d , 0 [ 2 d 1 ] selectively according to the last result?
Solution: The runtime of calculating all 2 d of y ^ 2 d , 0 [ i ] is 2 d log 2 d = d 2 d by using the FFT algorithm. The number of effective nodes in the upper layer is m d 1 , so the maximum number of effective nodes in this layer is 2 m d 1 . The formula for calculating each effective node is Equation (23), so the runtime of calculating these 2 m d 1 nodes is 2 m d 1 2 d . Therefore, when 2 m d 1 < = d , use Equation (23) to calculate 2 m d 1 nodes. Otherwise, use the FFT algorithm to calculate all nodes.
y ^ 2 d , 0 [ i ] = 1 2 d j = 0 2 d 1 ( D L S 0 x ) j ω 2 d i j
Subproblem 2: The condition of stopping the binary tree search.
Solution: If the stop condition is defined as m d = K , in the worst case, the condition cannot be satisfied until d is very large. For example, if x ^ 0 and x ^ N / 2 are significant frequencies, after the decomposition of the first layer’s node y ^ 1 , 0 [ 0 ] , the second layer’s node y ^ 2 , 0 [ 0 ] , the third layer’s node y ^ 4 , 0 [ 0 ] ,…, and their child node y ^ 2 , 0 [ 0 ] , y ^ 4 , 0 [ 0 ] , y ^ 8 , 0 [ 0 ] , are still aliasing until d = log N . Therefore, threshold = 0.75 K can be considered, which means that the K frequencies are put into 0.75 K buckets. The aliasing problem can be solved by using the SVD decomposition described in the last paragraph. Its extra sampling and calculating cost may be less than that of continuous expansion search if no more frequencies are separated by the next layer.
Subproblem 3: The frequency recovery problem of an approximately “single-ton” bucket.
Solution: In the exactly sparse case, the problem can be solved by using the phase encoding method described in the last paragraph. In the general sparse case, the problem can be solved by the CS method or the MMSE estimator described in the last paragraph.
The binary tree search method is especially suitable for small K or small support. In the case of small support, the frequencies must be scattered into every different bucket. For example, only eight consecutive frequencies of the 1028 frequencies x ^ 0 , , x ^ 1027 are significant frequencies, and their eight locations are equal to 0mod8, 1mod8, 2mod8, …, 7mod8. In this way, m 3 = 8 can be obtained in the fourth layer, and the binary tree search can be stopped by step 4. This method also has the advantage that it is a deterministic algorithm. Finally, the distribution of at most one frequency in one bucket must happen with the layer’s expansion, unlike some probabilistic algorithms such as sFFT-DT algorithms.

4. Algorithms Analysis in Theory

As mentioned above, the goal of frequency bucketization is to decrease runtime and sampling complexity in the advantage of low dimensions. After bucketization, the filtered spectrum y ^ L , τ , σ can be obtained by the original signal x. Then through the three frameworks introduced above, it is sufficient to recover the spectrum x ^ of the filtered spectrum y ^ L , τ , σ in their own way. In this section, we introduce and analyze the performance of corresponding algorithms, including the sFFT-DT1.0 algorithm, the sFFT-DT2.0 algorithm, the sFFT-DT3.0 algorithm, the FFAST algorithm, the R-FFAST algorithm, and the DSFFT algorithm. The theoretical analysis of the performance of the six algorithms can also be found in the relevant documents. The theoretical results of the sFFT-DT3.0, and DSFFT algorithms are given by the author for the first time, and the results of other algorithms are consistent with the original paper.

4.1. The sFFT-DT Algorithms by the One-Shot Framework

The block diagram of the sFFT algorithms by the one-shot framework is shown in Figure 1. We explain the details and analyze the performance in theory as below:
Stage 1 Bucketization: Run ( 3 a m + 3 P ) times’ round for set τ = { τ 0 = a m , , τ 3 a m 1 = 2 a m 1 } and set τ = τ 0 = r 1 , τ 1 = r 2 , , τ P 1 = r P , and calculate y ^ L , τ = F B D L S τ x . It costs ( 3 a m + 3 P ) B log B runtime and ( 3 a m + 3 P ) B samples.
Stage 2, Step 1: Obtain the number of significant frequencies in every bucket. By computing the SVD for the matrix M a m in every bucket, there are a m singular values for each bucket. Collect all B · a m singular values from all buckets and index each singular value. The first K largest singular values will indicate the number of significant frequencies in every bucket. It costs a m 3 B runtime.
Stage 2, Step 2: Locate the position by method 1: In the exactly sparse case, in a bucket waiting to be solved, first obtain C ( M a 1 ) M s , then obtain a number of roots of P ( z ) = 0 . The sFFT-DT1.0 algorithm uses this method, and it costs a m 2 K runtime in this step under the assumption of a m 4 .
Stage 2, Step 2: Locate the position by method 2: In a bucket waiting to be solved, first obtain C ( M a 1 ) M s , try every possible z j = ω f j in Equation (12), the first a number of smallest | P ( z ) | is the solution. The sFFT-DT2.0 algorithm uses this method, and it costs a m K L runtime in this step.
Stage 2, Step 2: Locate the position by method 3: In a bucket waiting to be solved, first compute the SVD of the matrix Y : Y = V ˜ Σ V * , then obtain the set of nonzero 2 ordinary eigenvalues of ( V 1 + ) V 2 . The sFFT-DT3.0 algorithm uses this method, and it costs a m 3 K runtime in this step.
Stage 2, Step 3: Estimate the value by CS solver: In a bucket waiting to be solved, (1) through the positions f j ’s gained by step 2, obtain the possible value of z j ’s as follows, z j , z j , z j , then obtain 3a vectors. The dictionary D contains 3 a vectors listed above. (2) From 3 a atoms of the dictionary matrix D , find a number of atoms that best matches the measurements’ residual error. Select these a number of atoms to construct a new sensing matrix Φ ^ . (3) Obtain S ^ with the support of S ^ = argmin y Φ ^ S ^ 2 through the LSM. (4) If the residual error meets the requirement, or the number of iterations reaches the assumption, or the residual error becomes larger, the iteration will be quit; otherwise continue to step 2. It costs 3 a m 2 P K runtime in this step.
Theorem 4. 
Suppose a m = 4 , B = 32 K , L = N / 32 K , P = 3 a m = 12 and N / K 32,000, it costs O ( K l o g K ) runtime and O ( K ) samples in the sFFT-DT1.0 algorithm. It costs O ( K l o g K + N ) runtime and O ( K ) samples in the sFFT-DT2.0 algorithm. It costs O ( K l o g K ) runtime and O ( K ) samples in the sFFT-DT3.0 algorithm.
Proof of Theorem 4. 
In the sFFT-DT1.0 algorithm, it costs in total ( 3 a m + 3 P ) B log B + a m 2 K + 3 a m 2 P K = O ( K log K ) runtime and ( 3 a m + 3 P ) B = O ( K ) samples. In the sFFT-DT2.0 algorithm, it costs in total ( 3 a m + 3 P ) B log B + a m K L + 3 a m 2 P K = O ( K log K + N ) runtime and ( 3 a m + 3 P ) B = O ( K ) samples. In the sFFT-DT3.0 algorithm, it costs in total ( 3 a m + 3 P ) B log B + a m 3 K + 3 a m 2 P K = O ( K log K ) runtime and ( 3 a m + 3 P ) B = O ( K ) samples. Notes: Theorem 4 has an assumption that P a m log 10 L because of the necessary number of measurements for CS recovery. This means that if N / 32 K 10 3 , B cannot continue to maintain itself equal to 32 K , and it has to increase to ensure that L = N / B 10 3

4.2. The FFAST Algorithms by the Peeling Framework

The block diagram of the sFFT algorithms by the peeling framework is shown in Figure 2. We explain the details and analyze the performance in theory as below:
Stage 1 Encoding by bucketization: With the assumption of K < N 1 / 3 , we run three bucketization cycles, and each cycle needs R rounds. The number of buckets is equal to B 1 , B 2 , B 3 individually in three cycles, and correspondingly the size of one bucket is equal to L 1 , L 2 , L 3 individually in three cycles. In the exactly sparse case, in the first peeling, we use the FFAST algorithm with R = 3 and the set τ = τ 0 = 0 , τ 1 = 1 , τ 2 = 2 in three rounds. In the general sparse case, in the first peeling, we use the R-FFAST algorithm with R = C m = O ( log 4 / 3 N ) and the set τ = r 0 , r 0 + 1 , , r 1 , r 1 + 2 , , r C 1 , r C 1 + 2 C 1 , , r C 1 + 2 C 1 ( m 1 ) in R rounds. After bucketization, we obtain the vectors y B 1 , 0 , , y B 1 , B 1 1 , y B 2 , 0 , , y B 2 , B 2 1 , y B 3 , 0 , , y B 3 , B 3 1 . It costs 3 ( B 1 log B 1 + B 2 log B 2 + B 3 log B 3 ) runtime and 3 ( B 1 + B 2 + B 3 ) samples in stage 1 for the FFAST algorithm. It costs ( log 4 / 3 N ) ( B 1 log B 1 + B 2 log B 2 + B 3 log B 3 ) runtime and ( log 4 / 3 N ) ( B 1 + B 2 + B 3 ) samples in stage 1 for the R-FFAST algorithm.
Stage 2—Decoding by several peelings: The subsequent calculation complexity is mainly composed of two parts: the judgment of the old bucket and the generation of the new bucket. The runtime of the generation of the new bucket can be ignored, so the computational complexity is mainly determined by the runtime of bucket judgment. Since the main buckets are identified as a “zero-ton” bucket or a “single-ton” bucket during the first peeling, the number of buckets that need to be calculated in the subsequent peeling iteration is not large, so the runtime is mainly determined at the first peeling. For the FFAST algorithm, there are ( B 1 + B 2 + B 3 ) buckets to be judged in the first peeling, and the runtime of each bucket is O ( 4 ) , as proven in Section 3. For the R-FFAST algorithm, the runtime of each bucket is O ( 4 log 4 / 3 N ) , as proven in Section 3.
Theorem 5. 
With the assumption of K < N 1 / 3 , suppose B 1 = O ( K ) , B 2 = O ( K ) , B 3 = O ( K ) ; it costs O ( K l o g K ) runtime and O ( K ) samples in the FFAST algorithm. It costs O ( K l o g 7 / 3 N ) runtime and O ( K l o g 4 / 3 N ) samples in the R-FFAST algorithm.
Proof of Theorem 5. 
For the FFAST algorithm, it costs in total 3 ( B 1 log B 1 + B 2 log B 2 + B 3 log B 3 ) + 4 ( B 1 + B 2 + B 3 ) = O ( K log K ) runtime and 3 ( B 1 + B 2 + B 3 ) = O ( K ) samples. For the R-FFAST algorithm, it costs in total ( log 4 / 3 N ) ( B 1 log B 1 + B 2 log B 2 + B 3 log B 3 ) + 4 ( log 4 / 3 N ) ( B 1 + B 2 + B 3 ) = O ( K log 7 / 3 N ) runtime and ( log 4 / 3 N ) ( B 1 + B 2 + B 3 ) = O ( K log 4 / 3 N ) samples. Notes: From the paper [12], with the assumption of K < N 1 / 3 , if the number of buckets B 1 , B 2 , B 3 is a co-prime number and approximately equal to O ( K ) , the chain number R = O ( log 4 / 3 N ) , and in every ”single-ton” estimator, the number of iterations C = O ( log N ) , and the number of rounds per iteration m = O ( log 1 / 3 N ) , then the algorithm can correctly recover the spectrum with a high probability (probability of successful ( 1 1 / K ) ) . □

4.3. The DSFFT Algorithm by the Binary Tree Search Framework

The first stage of the DSFFT algorithm is finding exactly O ( K ) efficient buckets through the binary tree search method. When the m d gained by step ( d + 1 ) is equal to the sparse number K, the binary tree search is finished. Now we start to calculate the probability in this case. The total number of buckets is 2 d = B , and the number of aliasing frequencies in each bucket is L = N / B = N / 2 d . In other words, the right case is K large frequencies of totally N frequencies allocated to B buckets, and there is no aliasing in the B distributed buckets. The probability P is calculated as follows: First, analyze the first bucket, the probability P 1 is the case of the non-aliasing case divided by the total case, P 1 = ( N K ) · ( N K 1 ) ( N K L + 2 ) ( N 1 ) · ( N 2 ) ( N L + 1 ) . Second, analyze the second bucket, the probability P 2 = ( N K L + 1 ) · ( N K L ) ( N K 2 L + 3 ) ( N L 1 ) · ( N L 2 ) ( N 2 L + 1 ) , . Finally, analyze the no. K’s bucket, P K = ( N ( K 1 ) L 1 ) ( N K L + 1 ) ( N ( K 1 ) L 1 ) ( N K L + 1 ) . Thus the complete probability P = P 1 P 2 P K = ( N L ) ( N 2 L ) ( N ( K 1 ) L ) ( N 1 ) ( N 2 ) ( N K + 1 ) = L ( K 1 ) ( ( B 1 ) ! ) ( ( N K ) ! ) ( ( N 1 ) ! ) ( ( B K ) ! ) . For example, consider N = 1000 , K = 5 , B = 10 , L = 100 and the complete probability P = ( 995 × 897 999 × 901 ) ( 599 × 501 599 × 501 ) = 995 · 990 · 985 · 980 999 · 998 · 997 · 996 . Table 1 shows the probability P B 0 , P B 1 , P B 2 , P B 3 of different K in the case of B 0 = K , B 1 = 2 K , B 2 = 4 K , B 3 = 8 K . When B 0 = K , P B 0 = ( N K ) ( K 1 ) ( ( K 1 ) ! ) ( ( N K ) ! ) ( N 1 ) ! . When B 1 = 2 K , P B 1 = ( N 2 K ) ( K 1 ) ( ( 2 K 1 ) ! ) ( ( N K ) ! ) ( ( N 1 ) ! ) ( ( K ) ! ) . When B 2 = 4 K , P B 2 = ( N 4 K ) ( K 1 ) ( ( 4 K 1 ) ! ) ( ( N K ) ! ) ( ( N 1 ) ! ) ( ( 3 K ) ! ) .
As we can see from Table 1, the DSFFT algorithm has high efficiency only when K is very small. The total complexity can be calculated according to the sum of probability times conditional complexity. Thus the total runtime complexity is P B 0 K log K + ( P B 1 P B 0 ) 2 K log 2 K + ( P B 2 P B 1 ) 4 K log 4 K + = j = 0 j = log ( N / K ) ( P B j P B j 1 ) 2 j K log ( 2 j K ) for P B 1 = 0 , and the total sampling complexity is P B 0 K + ( P B 1 P B 0 ) 2 K + ( P B 2 P B 1 ) 4 K + = j = 0 j = log ( N / K ) ( P B j P B j 1 ) 2 j K for P B 1 = 0 .

4.4. Summary of Six Algorithms in Theory

After analyzing three types of frameworks, six corresponding algorithms, Table 2 (the performance of algorithms using the aliasing filter is got as above. The performance of other algorithms is got from paper [7]. The analysis of robustness will be explained in the next section.) can be concluded with the additionl information of other sFFT algorithms and fftw algorithm.
From Table 2, we can see that the sFFT-DT1.0 algorithm and the FFAST algorithm have the lowest runtime and sampling complexity, but they are non-robust. Other algorithms using the aliasing filter have good sampling complexity but compared to the other sFFT algorithms they posses no advantage in the runtime complexity except for the sFFT-DT3.0 algorithm. In addition, the use of the aliasing algorithm is limited in some aspects. For the algorithms of the one-shot platform, the performance is not very efficient if N / K 32,000 because the necessary number of measurements is needed for CS recovery. For the algorithms of the peeling platform, the assumption of K < N 1 / 3 is required because if there are four or more cycles, it will be very complicated. For the algorithm of the binary tree search framework, K is required to be very small, because only in this case, there is no need to expand the scale of a large binary tree.
As to different filters, the algorithms using the aliasing filter have good sampling because of their very short support set, and they also have no spectrum leakage because of properties of the filter. The algorithms using the flat filter have good time complexity based on the simple calculation with a short support set. The algorithms using the Dirichlet kernel filter bank are the most complicated but can support random sampling on account of their estimation method.

5. Algorithms Analysis in Practice

In this section, we evaluate the performance of six sFFT algorithms using the aliasing filter: the sFFT-DT1.0, sFFT-DT2.0, sFFT-DT3.0, FFAST, R-FFAST, and DSFFT algorithms. We first compare these algorithms’ runtime, the percentage of the signal sampled and robustness characteristics. Then we compare some of these algorithms’ characteristics with other algorithms: the fftw, sFFT1.0, sFFT2.0, sFFT-3.0, sFFT-4.0, and AAFFT algorithms. The codes of the sFFT1.0, sFFT2.0 (the code is available at http://groups.csail.mit.edu/netmit/sFFT/), sFFT3.0, MPFFT (the code is available at https://github.com/urrfinjuss/mpfft), sFFT-DT2.0 (the code is available at https://www.iis.sinica.edu.tw/pages/lcs), R-FFAST (the code is available at https://github.com/UCBASiCS/FFAST), AAFFT (the code is available at https://sourceforge.net/projects/aafftannarborfa/), and fftw (the code is available at http://www.fftw.org/) algorithms are already open source. All experiments were run on a Linux CentOS computer with a 4 Intel(R) Core(TM) i5-4590 3.3 GGHz CPU and 8 GB of RAM.

5.1. Experimental Setup

In the experiment, the test signals were gained in a way that K frequencies were randomly selected from N frequencies and assigned a magnitude of 1 and a uniformly random phase. The remaining frequencies were set to zero in the exactly sparse case or combined with additive white Gaussian noise in the general sparse case, whose variance varies depending on the SNR required. The parameters of these algorithms were chosen so that they can make a balance between time efficiency and robustness. In each experiment, the platform could generate a signal with SNR, K, or N as required. The prepared signal and the value of N and K were transmitted to different algorithm libraries through a standard interface. The results of the algorithm library, the time and sampling ratio used in operation were also be returned to the platform. These algorithm libraries were generated by modifying the open source code of the implementation of the algorithms mentioned above. Each test record contains the run time, the sampling proportion, and the L 0 , L 1 , L 2 error between the calculation result and the best result through the algorithm library. The details of code, data, and reports are all open source (the R-FFAST algorithm does not pass through this platform because its signal length is required to be the product of several prime numbers). The new testing platform was developed from an old platform (https://github.com/ludwigschmidt/sft-experiments. The detail of codes, data, and reports are all open sources (https://github.com/zkjiang/-/tree/master/docs/sfftproject.

5.2. Comparison Experiment of Different Algorithms Using the Aliasing Filter

We plot Figure 7a,b representing run time vs. signal size and vs. signal sparsity for the sFFT-DT2.0, sFFT-DT3.0, R-FFAST, and DSFFT algorithms in the general sparse case (the general sparse case means SNR = 20 db. The exact case was very similar to the general sparse case except in the sFFT-DT1.0 algorithm and the FFAST algorithm. As mentioned above, the runtime is determined by two factors. One is how many buckets are there to cope with, and another is how much time does it cost to identify frequencies in efficient buckets. Thus, from Figure 7, we can see that: (1) The runtimes of these four algorithms were approximately linear in the log scale as a function of N and non-linear as a function of K. (2) The result of ranking the runtime of four algorithms was sFFT-DT3.0 > sFFT-DT2.0 > DSFFT > R-FFAST. The reason is the method of the R-FFAST algorithm is the most time-consuming, the method of DSFFT algorithm is also not efficient, the method of sFFT-DT2.0 is much better, and the method of sFFT-DT3.0 had the highest time performance. (3) K has a certain limitation. The DSFFT algorithm is very inefficient when K is greater than 50, and the R-FFAST algorithm does not work when K is greater than N 1 / 3 .
We plot Figure 8a,b representing the percentage of the signal sampled vs. signal size and vs. signal sparsity for the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithms in the general sparse case. As mentioned above, the percentage of the signal sampled is also determined by two factors: how many buckets there are and how many samples are sampled in one bucket. Thus from Figure 8, we can see the following: (1) The percentages of the signal sampled of these three algorithms were approximately linear in the log scale as a function of N and non-linear as a function of K. (2) The result of ranking the sampling complexity of the three algorithms is R-FFAST > sFFT-DT2.0 = sFFT-DT3.0 because the R-FFAST algorithm uses the principle of CRT, the number of buckets is independent of K, and the proportion decreases with the increase of N. (3) There is a limit for the sFFT-DT2.0 and sFFT-DT3.0 algorithms. When N / K is too large, B cannot maintain the linearity of the sparsity K. For the R-FFAST algorithm, there is a limit in which K cannot be greater than N 1 / 3 .
We plot Figure 9a,b representing run time and L 1 -error vs. SNR for the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithms. From Figure 9 we can see that: (1) The runtime of sFFT-DT2.0 and sFFT-DT3.0 was approximately equal vs. SNR, but the runtime of the R-FFAST algorithm increased with the noise. (2) To a certain extent, these three algorithms are all robust. When SNR was less than 0 db, only the R-FFAST algorithm satisfied the robustness. When SNR was medium, the sFFT-DT3.0 algorithm also met the robustness. Only when SNR was bigger than 20 db, could the sFFT-DT2.0 algorithm deal with the noise interference. The reason is that the method of binary search is better than other ways in terms of robustness.

5.3. Comparison Experiment of Algorithms Using the Aliasing Filter and Other Algorithms

We plot Figure 10a,b representing run time vs. signal size and vs. signal sparsity for the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, AAFF,T and fftw algorithms in the general sparse case (the sFFT1.0 algorithm is ignored because it is similar to the sFFT2.0 algorithm but not as efficient as it. The sFFT3.0 algorithm is ignored because it is not suitable for general sparsity. The MPSFT algorithm is ignored because it is similar to the sFFT4.0 algorithm but not as efficient as it). From Figure 10, we can see that: (1) These algorithms are approximately linear in the log scale as a function of N, except for the fftw algorithm. These algorithms are approximately linear in the standard scale as a function of K, except for the fftw and sFFT-DT3.0 algorithms. (2) The result of ranking the runtime complexity of these six algorithms was sFFT2.0 > sFFT4.0 > AAFFT > sFFT-DT3.0 > fftw > R-FFAST when N was large. The reason is the least absolute shrinkage and selection operator (LASSO) method used in R-FFAST costs a large amount of time. Additionally, the SVD method and the CS method used in the sFFT-DT also cost a large amount of time. (2) The result of ranking the runtime complexity of these six algorithms was fftw > sFFT-DT3.0 > sFFT4.0 > sFFT2.0 > AAFFT > R-FFAST when K was large. The reason is algorithms using the aliasing filter need much fewer buckets than algorithms using the flat filter when K is large.
We plot Figure 11a,b representing the percentage of the signal sampled vs. signal size and vs. signal sparsity for the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, AAFFT, and fftw algorithms in the general sparse case. From Figure 11, we can see that: (1) These algorithms are approximately linear in the log scale as a function of N, except for the fftw and sFFT-DT algorithms. The reason is that sampling in low dimensions in the sFFT algorithms can decrease sampling complexity and it forms a limit to the size of the bucket in the sFFT-DT algorithm by using the CS method. These algorithms are approximately linear in the standard scale as a function of K except for the R-FFAST, sFFT-DT, and fftw algorithms. The reason is that algorithms using the aliasing filter save samples by using a smaller number of buckets. (2) The result of ranking the sampling complexity of these six algorithms was R-FFAST > sFFT4.0 > AAFFT > sFFT2.0 > sFFT-DT3.0 > fftw when N was large. (3) The result of ranking the sampling complexity was sFFT-DT3.0 > sFFT4.0 > AAFFT > sFFT2.0 > fftw when K was large.
We plot Figure 12a,b representing run time and L 1 -error vs. SNR for the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, AAFFT, and fftw algorithms. From Figure 12, we can see that: (1) The runtime is approximately equal vs. SNR except the R-FFAST algorithm. (2) To a certain extent, these six algorithms are all robust, but when SNR was less than 0 db, only the fftw and R-FFAST algorithms satisfied the robustness. When SNR is medium, the sFFT2.0, AAFFT, and sFFT-DT3.0 algorithms can also meet the robustness. Only when SNR is bigger than 20 dB could the sFFT4.0 algorithm can deal with the noise interference.

6. Conclusions

The development of new technology of signal processing and the demand for big data computing and low sampling ratio of DFT accelerate the research of sparse Fourier algorithms. There are three types of filters used in various algorithms. The typical algorithms (sFFT1.0, SFFT2.0, sFFT3.0, sFFT4.0, MPSFT) based on the flat filter were proposed by MIT (Massachusetts Institute of Technology). The typical algorithms (sFFT-DT1.0, SFFT-DT2.0, FFAST, R-FFAST, DSFFT) based on the aliasing filter were proposed by Taiwan university, Berkeley university, and Georg-August University respectively. The typical algorithms (AAFFT0.5, AAFFT0.9) based on the Dirichlet kernel filter bank window were proposed by Michigan State University. This paper mainly discussed all of the theoretical techniques of algorithms using the aliasing filter and makes a complete comparative analysis between them and other algorithms.
In the first part, the paper introduced the techniques used in the sFFT algorithms including time-shift operation and subsampling operation. In the second part, we analyzed in detail six typical algorithms using the aliasing filter. By the one-shot framework based on the CS solver, the sFFT-DT1.0 algorithm uses the polynomial method, the sFFT-DT2.0 algorithm uses the enumeration method and the sFFT-DT3.0 algorithm uses the matrix pencil method. By the peeling framework based on the bipartite graph, the FFAST algorithm uses the phase encoding method and the R-FFAST algorithm uses the MMSE estimator method. By the iterative framework, the DSFFT algorithm uses the binary tree search method. We obtained the conclusion of the performance of these algorithms including runtime complexity, sampling complexity, and robustness in theory, as detailed in Table 2. In the third part, we made two categories of experiments for computing the signals of different SNRs, different N and different K using a standard testing platform and recorded the run time, percentage of the signal sampled, and L 0 , L 1 , L 2 error by using nine different sFFT algorithms both in the exactly sparse case and the general sparse case. The analysis of the experimental results satisfies theoretical inference. In the comparison of algorithms using different frameworks, the algorithms of one-shot framework had good sampling complexity because their support set is based on O ( K ) . The algorithms of peeling frameworks had good robustness because they use the binary search method. The algorithms of iterative frameworks were very inefficient when K was greater than 50. In the comparison of algorithms using different filters, the algorithms using the Dirichlet kernel filter were the inefficient deterministic algorithms. The algorithms using the flat filter had good running complexity because it is easy to calculate. The algorithms using the aliasing filter had good sampling complexity because of their very short support set.
The main contribution of this paper are as follows: (1) Three platforms and all related technologies of sFFT algorithms using the aliasing filter were analyzed in detail. Based on the above introduction, the theoretical performance of different algorithms using corresponding techniques was derived and the comparative studies could also be carried out as well. (2) A total of nine typical sFFT algorithms were implemented in the form of an fftw library, and the interface and codes are all open sources. We developed a standard testing platform, which can test more than ten typical sFFT algorithms in different situations on the basis of these algorithm libraries. (3) We obtained conclusions regarding the character and performance of the six typical sFFT algorithms using the aliasing filter. The comparison results of sFFT algorithms using the aliasing filter based on different platforms were obtained, and the comparison results of this type of algorithm and other algorithms were also obtained. (4) In addition to the above theoretical, software and practical contributions, other novel points of this paper are as follows: The sFFT-DT3.0 algorithm was proposed and the sFFT4.0 algorithm was implemented by the authors for the first time. The complete comparative analysis of six different typical sFFT algorithms using the aliasing filter and other algorithms was obtained step by step in theory and in practice for the first time. In all, the systematic analysis of different algorithms in this paper will be helpful for readers outside the signal processing community who wish to obtain a solid idea regardin SFFT algorithms using the aliasing filter and use them correctly and efficiently.

Author Contributions

Conceptualization, Z.J.; methodology, Z.J.; software, Z.J.; validation, Z.J.; formal analysis, Z.J.; investigation, Z.J.; resources, Z.J.; data curation, Z.J.; writing—original draft preparation, Z.J.; writing—review and editing, B.L.; visualization, Z.J.; supervision, Z.J.; project administration, B.L.; funding acquisition, J.C. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Youth Program of National Natural Science Foundation of China under Grant 61703263.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cooley, J.W.; Tukey, J.W. An Algorithm for the Machine Calculation of Complex Fourier Series. Math. Comput. 1965, 19, 297. [Google Scholar] [CrossRef]
  2. Gilbert, A.C.; Guha, S.; Indyk, P.; Muthukrishnan, S.; Strauss, M. Near-optimal sparse Fourier representations via sampling. Conf. Proc. Annu. ACM Symp. Theory Comput. 2002, 2, 152–161. [Google Scholar] [CrossRef] [Green Version]
  3. Iwen, M.A.; Gilbert, A.; Strauss, M. Empirical evaluation of a sub-linear time sparse DFT algorithm. Commun. Math. Sci. 2007, 5, 981–998. [Google Scholar] [CrossRef]
  4. Gilbert, A.C.; Strauss, M.J.; Tropp, J.A. A tutorial on fast Fourier sampling: How to apply it to problems. IEEE Signal Process. Mag. 2008, 25, 57–66. [Google Scholar] [CrossRef]
  5. Hassanieh, H.; Indyk, P.; Katabi, D.; Price, E. Simple and practical algorithm for sparse Fourier transform. In Proceedings of the Annual ACM-SIAM Symposium on Discrete Algorithms, Kyoto, Japan, 17–19 January 2012; pp. 1183–1194. [Google Scholar] [CrossRef] [Green Version]
  6. Hassanieh, H.; Indyk, P.; Katabi, D.; Price, E. Nearly optimal sparse Fourier transform. In Proceedings of the Annual ACM Symposium on Theory of Computing, New York, NY, USA, 20–22 May 2012; pp. 563–577. [Google Scholar] [CrossRef] [Green Version]
  7. Chiu, J. Matrix Probing, Skeleton Decompositions, and Sparse Fourier Transform. Ph.D. Thesis, Massachusetts Institute of Technology, Department of Mathematics, Cambridge, MA, USA, 2013. [Google Scholar] [CrossRef]
  8. Li, B.; Jiang, Z.; Chen, J. On performance of sparse fast Fourier transform algorithms using the flat window filter. IEEE Access 2020, 8, 79134–79146. [Google Scholar] [CrossRef]
  9. Hsieh, S.H.; Lu, C.S.; Pei, S.C. Sparse Fast Fourier Transform by downsampling. In Proceedings of the ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing—Proceedings, Vancouver, BC, Canada, 26–31 May 2013; pp. 5637–5641. [Google Scholar] [CrossRef] [Green Version]
  10. Hsieh, S.h.; Lu, C.s.; Pei, S.c. Sparse Fast Fourier Transform for Exactly and Generally. arXiv 2015, arXiv:1407.8315. [Google Scholar]
  11. Pawar, S.; Ramchandran, K. Computing a k-sparse n-length Discrete Fourier Transform using at most 4k samples and O(k log k) complexity. In Proceedings of the IEEE International Symposium on Information Theory—Proceedings, Istanbul, Turkey, 7–12 July 2013; pp. 464–468. [Google Scholar] [CrossRef] [Green Version]
  12. Pawar, S.; Ramchandran, K. FFAST: An algorithm for computing an exactly k-Sparse DFT in O(k log k) time. IEEE Trans. Inf. Theory 2018, 64, 429–450. [Google Scholar] [CrossRef]
  13. Pawar, S.; Ramchandran, K. A robust sub-linear time R-FFAST algorithm for computing a sparse DFT. arXiv 2015, arXiv:1501.00320. [Google Scholar] [CrossRef]
  14. Ong, F.; Heckel, R.; Ramchandran, K. A Fast and Robust Paradigm for Fourier Compressed Sensing Based on Coded Sampling. In Proceedings of the ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing—Proceedings, Brighton, UK, 12–17 May 2019; pp. 5117–5121. [Google Scholar] [CrossRef]
  15. Plonka, G.; Wannenwetsch, K.; Cuyt, A.; Lee, W.S. Deterministic sparse FFT for M-sparse vectors. Numer. Algorithms 2018. [Google Scholar] [CrossRef]
  16. Iwen, M.A. Combinatorial sublinear-time Fourier algorithms. Found. Comput. Math. 2010, 10, 303–338. [Google Scholar] [CrossRef]
  17. Iwen, M.A. Improved approximation guarantees for sublinear-time Fourier algorithms. Appl. Comput. Harmon. Anal. 2013, 34, 57–82. [Google Scholar] [CrossRef] [Green Version]
  18. Lawlor, D.; Wang, Y.; Christlieb, A. Adaptive Sub-Linear Time Fourier Algorithms. arXiv 2013, arXiv:1207.6368. [Google Scholar] [CrossRef]
  19. Christlieb, A.; Lawlor, D.; Wang, Y. A multiscale sub-linear time Fourier algorithm for noisy data. Appl. Comput. Harmon. Anal. 2016, 40, 553–574. [Google Scholar] [CrossRef] [Green Version]
  20. Merhi, S.; Zhang, R.; Iwen, M.A.; Christlieb, A. A New Class of Fully Discrete Sparse Fourier Transforms: Faster Stable Implementations with Guarantees. J. Fourier Anal. Appl. 2019, 25, 751–784. [Google Scholar] [CrossRef] [Green Version]
  21. Gilbert, A.C.; Indyk, P.; Iwen, M.; Schmidt, L. Recent Developments in the Sparse Fourier Transform. Signal Processing Mag. 2014, 31, 91–100. [Google Scholar] [CrossRef]
  22. Kapralov, M. Sample efficient estimation and recovery in sparse FFT via isolation on average. In Proceedings of the Annual Symposium on Foundations of Computer Science—Proceedings, Berkeley, CA, USA, 15–17 October 2017; pp. 651–662. [Google Scholar] [CrossRef] [Green Version]
  23. Indyk, P.; Kapralov, M.; Price, E. (Nearly) sample-optimal sparse Fourier transform. In Proceedings of the Annual ACM-SIAM Symposium on Discrete Algorithms, Portland, OR, USA, 5–7 January 2014; pp. 480–499. [Google Scholar] [CrossRef]
  24. López-Parrado, A.; Velasco Medina, J. Efficient Software Implementation of the Nearly Optimal Sparse Fast Fourier Transform for the Noisy Case. Ingeniería y Ciencia 2015, 11, 73–94. [Google Scholar] [CrossRef]
  25. Chen, G.L.; Tsai, S.H.; Yang, K.J. On Performance of Sparse Fast Fourier Transform and Enhancement Algorithm. IEEE Trans. Signal Process. 2017, 65, 5716–5729. [Google Scholar] [CrossRef]
  26. Schumacher, J.; Püschel, M. High-performance sparse fast Fourier transforms. In Proceedings of the IEEE Workshop on Signal Processing Systems, SiPS: Design and Implementation, Belfast, UK, 20–22 October 2014. [Google Scholar] [CrossRef] [Green Version]
  27. Wang, C. CusFFT: A High-Performance Sparse Fast Fourier Transform Algorithm on GPUs. In Proceedings of the Proceedings—2016 IEEE 30th International Parallel and Distributed Processing Symposium, IPDPS 2016, Chicago, IL, USA, 23–27 May 2016; pp. 963–972. [Google Scholar] [CrossRef]
  28. Abari, O.; Hamed, E.; Hassanieh, H.; Agarwal, A.; Katabi, D.; Chandrakasan, A.P.; Stojanovic, V. A 0.75-million-point Fourier-transform chip for frequency-sparse signals. In Proceedings of the Digest of Technical Papers—IEEE International Solid-State Circuits Conference, San Francisco, CA, USA, 9–13 February 2014. [Google Scholar] [CrossRef]
  29. Wang, S.; Patel, V.M.; Petropulu, A. Multidimensional Sparse Fourier Transform Based on the Fourier Projection-Slice Theorem. IEEE Trans. Signal Process. 2019, 67, 54–69. [Google Scholar] [CrossRef]
  30. Kapralov, M.; Velingker, A.; Zandieh, A. Dimension-independent sparse Fourier transform. In Proceedings of the Annual ACM-SIAM Symposium on Discrete Algorithms, San Diego, CA, USA, 6–9 January 2019. [Google Scholar] [CrossRef] [Green Version]
  31. Pang, C.; Liu, S.; Han, Y. High-speed target detection algorithm based on sparse Fourier transform. IEEE Access 2018, 6, 37828–37836. [Google Scholar] [CrossRef]
  32. Kumar, G.G.; Sahoo, S.K.; Meher, P.K. 50 Years of FFT Algorithms and Applications. Circuits Syst. Signal Process. 2019, 38, 5665–5698. [Google Scholar] [CrossRef]
  33. Plonka, G.; Wannenwetsch, K. A deterministic sparse FFT algorithm for vectors with small support. Numer. Algorithms 2016, 71, 889–905. [Google Scholar] [CrossRef] [Green Version]
  34. Plonka, G.; Wannenwetsch, K. A sparse fast Fourier algorithm for real non-negative vectors. J. Comput. Appl. Math. 2017, 321, 532–539. [Google Scholar] [CrossRef] [Green Version]
  35. Chien, R.T. Cyclic Decoding Procedures for Bose-Chaudhuri-Hocquenghem Codes. IEEE Trans. Inf. Theory 1964, 10, 357–363. [Google Scholar] [CrossRef]
  36. Bunse-Gerstner, A.; Gragg, W.B. Singular value decompositions of complex symmetric matrices. J. Comput. Appl. Math. 1988. [Google Scholar] [CrossRef]
  37. Hua, Y.; Sarkar, T.K. Matrix Pencil Method for Estimating Parameters of Exponentially Damped/Undamped Sinusoids in Noise. IEEE Trans. Acoust. Speech Signal Process. 1990. [Google Scholar] [CrossRef] [Green Version]
  38. Kay, S. A Fast and Accurate Single Frequency Estimator. IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 1987–1990. [Google Scholar] [CrossRef] [Green Version]
Figure 1. A system block diagram of the one-shot framework.
Figure 1. A system block diagram of the one-shot framework.
Electronics 10 01117 g001
Figure 2. A system block diagram of the peeling framework.
Figure 2. A system block diagram of the peeling framework.
Electronics 10 01117 g002
Figure 3. An example of a bipartite graph with N = 20 variable nodes on the left and N b = 9 parity check nodes on the right (blue square with “N” indicates a “zero-ton” bucket, green square with “S” indicates a “single-ton” bucket, red square with “M” indicates a “multi-ton” bucket).
Figure 3. An example of a bipartite graph with N = 20 variable nodes on the left and N b = 9 parity check nodes on the right (blue square with “N” indicates a “zero-ton” bucket, green square with “S” indicates a “single-ton” bucket, red square with “M” indicates a “multi-ton” bucket).
Electronics 10 01117 g003
Figure 4. Example of the process of the genie-assisted peeling decoder. In the first peeling, three “zero-ton”, three “single-ton”, and three “multi-ton” are distinguished. In the second peeling, two “single-ton” and one ”multi-ton” are distinguished. In the third peeling, one “zero-ton” is distinguished.
Figure 4. Example of the process of the genie-assisted peeling decoder. In the first peeling, three “zero-ton”, three “single-ton”, and three “multi-ton” are distinguished. In the second peeling, two “single-ton” and one ”multi-ton” are distinguished. In the third peeling, one “zero-ton” is distinguished.
Electronics 10 01117 g004
Figure 5. Example of the process of the sparse-graph decoder by the packet erasure channel method. In the process of expanding the subgraph of x ^ 2 , when the depth is equal to two, one “zero-ton” and one “single-ton” are distinguished; when the depth is equal to four, one “multi-ton” is distinguished; when the depth is equal to six, one “multi-ton” is distinguished. In the process of expanding the subgraph of x ^ 7 , one “zero-ton”, one “single-ton”, and two “multi-ton”s are distinguished.
Figure 5. Example of the process of the sparse-graph decoder by the packet erasure channel method. In the process of expanding the subgraph of x ^ 2 , when the depth is equal to two, one “zero-ton” and one “single-ton” are distinguished; when the depth is equal to four, one “multi-ton” is distinguished; when the depth is equal to six, one “multi-ton” is distinguished. In the process of expanding the subgraph of x ^ 7 , one “zero-ton”, one “single-ton”, and two “multi-ton”s are distinguished.
Electronics 10 01117 g005
Figure 6. Example of the iterative framework based on the binary tree search method. The number of the efficient nodes of the first, second, third, and fourth layers of the binary tree is equal to one, two, four, and five, respectively.
Figure 6. Example of the iterative framework based on the binary tree search method. The number of the efficient nodes of the first, second, third, and fourth layers of the binary tree is equal to one, two, four, and five, respectively.
Electronics 10 01117 g006
Figure 7. Run time of four algorithms using the aliasing filter, including the sFFT-DT2.0, sFFT-DT3.0, R-FFAST, and DSFFT algorithms in the general sparse case.
Figure 7. Run time of four algorithms using the aliasing filter, including the sFFT-DT2.0, sFFT-DT3.0, R-FFAST, and DSFFT algorithms in the general sparse case.
Electronics 10 01117 g007
Figure 8. Percentage of the signal sampled of three algorithms using the aliasing filter, including the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithms in the general sparse case.
Figure 8. Percentage of the signal sampled of three algorithms using the aliasing filter, including the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithms in the general sparse case.
Electronics 10 01117 g008
Figure 9. Run time and L 1 -error of three algorithms using the aliasing filter, including the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithms in the general sparse case vs. SNR.
Figure 9. Run time and L 1 -error of three algorithms using the aliasing filter, including the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithms in the general sparse case vs. SNR.
Electronics 10 01117 g009
Figure 10. Runtime of two typical algorithms using the aliasing filter and four other algorithms, including the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, fftw, and AAFFT0.9 algorithms in the general sparse case.
Figure 10. Runtime of two typical algorithms using the aliasing filter and four other algorithms, including the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, fftw, and AAFFT0.9 algorithms in the general sparse case.
Electronics 10 01117 g010
Figure 11. Percentage of the signal sampled of two typical algorithms using the aliasing filter and four other algorithms, including the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, fftw, and AAFFT0.9 algorithms in the general sparse case.
Figure 11. Percentage of the signal sampled of two typical algorithms using the aliasing filter and four other algorithms, including the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, fftw, and AAFFT0.9 algorithms in the general sparse case.
Electronics 10 01117 g011
Figure 12. Run time and L 1 -error of two typical algorithms using the aliasing filter and four other algorithms, including the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, fftw, and AAFFT0.9 algorithms in the general sparse case vs. SNR.
Figure 12. Run time and L 1 -error of two typical algorithms using the aliasing filter and four other algorithms, including the sFFT-DT3.0, R-FFAST, sFFT2.0, sFFT4.0, fftw, and AAFFT0.9 algorithms in the general sparse case vs. SNR.
Electronics 10 01117 g012
Table 1. The probability P B 0 , P B 1 , P B 2 , P B 3 of different K in the case of B from K to 8 K .
Table 1. The probability P B 0 , P B 1 , P B 2 , P B 3 of different K in the case of B from K to 8 K .
Probability K = 1 K = 2 K = 4 K = 6 K = 8
P B 0 1 1 2 = 0.5 3 ! 4 3 0.0938 5 ! 6 5 0.0154 7 ! 8 7 0.0024
P B 1 1 3 4 = 0.75 7 ! 8 3 4 ! 0.4101 11 ! 12 5 6 ! 0.2228 15 ! 16 7 8 ! 0.1208
P B 2 1 7 8 = 0.875 15 ! 16 3 12 ! 0.6665 23 ! 24 5 18 ! 0.5071 31 ! 32 7 24 ! 0.3857
P B 3 1 15 16 = 0.9375 31 ! 32 3 28 ! 0.8231 47 ! 48 5 42 ! 0.7224 63 ! 64 7 56 ! 0.6340
Table 2. The performance of fftw algorithm and sFFT algorithms in theory.
Table 2. The performance of fftw algorithm and sFFT algorithms in theory.
AlgorithmRuntime ComplexitySampling ComplexityRobustness
sFFT1.0 O ( K 1 2 N 1 2 log 3 2 N ) N 1 N w N log N medium
sFFT2.0 O ( K 2 3 N 1 3 log 4 3 N ) N 1 N w N log N medium
sFFT3.0 O ( K log N ) O ( K log N ) none
sFFT4.0 O ( K log N log 8 ( N / K ) ) O ( K log N log 8 ( N / K ) ) bad
MPFFT O ( K log N log 2 ( N / K ) ) O ( K log N log 2 ( N / K ) ) good
sFFT-DT1.0 O ( K log K ) O ( K ) none
sFFT-DT2.0 O ( K log K + N ) O ( K ) medium
sFFT-DT3.0 O ( K log K ) O ( K ) medium
FFAST O ( K log K ) O ( K ) none
R-FFAST O ( K log 7 / 3 N ) O ( K log 4 / 3 K ) good
DSFFT j = 0 j = log ( N / K ) ( P B j P B j 1 ) 2 j K log ( 2 j K ) j = 0 j = log ( N / K ) ( P B j P B j 1 ) 2 j K medium
AAFFT O ( K poly ( log N ) ) O ( K poly ( log N ) ) medium
fftw O ( N log N ) Ngood
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, B.; Jiang, Z.; Chen, J. On Performance of Sparse Fast Fourier Transform Algorithms Using the Aliasing Filter. Electronics 2021, 10, 1117. https://0-doi-org.brum.beds.ac.uk/10.3390/electronics10091117

AMA Style

Li B, Jiang Z, Chen J. On Performance of Sparse Fast Fourier Transform Algorithms Using the Aliasing Filter. Electronics. 2021; 10(9):1117. https://0-doi-org.brum.beds.ac.uk/10.3390/electronics10091117

Chicago/Turabian Style

Li, Bin, Zhikang Jiang, and Jie Chen. 2021. "On Performance of Sparse Fast Fourier Transform Algorithms Using the Aliasing Filter" Electronics 10, no. 9: 1117. https://0-doi-org.brum.beds.ac.uk/10.3390/electronics10091117

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop