Two Efﬁcient Sparse Fourier Algorithms Using the Matrix Pencil Method

: The research on efﬁcient computation of sparse signals by various Sparse Fast Fourier Transform (sFFT) algorithms has always been a hot topic in the direction of signal processing. The algorithms can decrease the sampling and running complexity by taking advantage of the signal’s inherent characteristics that a large number of signals are sparse in the frequency domain. The sFFT algorithm is generally divided into two stages: the ﬁrst stage is bucketization. The process is to divide N frequencies into B buckets through the ﬁlter. The main ﬁlters used are the ﬂat window ﬁlter and the aliasing ﬁlter. The second stage is the spectrum recovery. The process is to successfully locate the position of the large frequency in each bucket and successfully calculate the amplitude. Among these steps, the most difﬁcult and time-consuming step is to successfully locate the position of the large value frequency. For the sFFT algorithm based on a ﬂat window ﬁlter, the voting method used in sFFT1.0 and sFFT2.0 algorithms should use many rounds, so these two algorithms are time-consuming and indeterministic, while the phase estimation method used in sFFT3.0 and sFFT4.0 algorithms has medium robustness. For sFFT algorithms based on an aliasing ﬁlter, the Prony method used in sFFT-DT1.0 algorithm is only applicable to noiseless signals, while the enumeration method used in the sFFT-DT2.0 algorithm has high complexity and poor robustness. In view of the performance of the old methods, new and more efﬁcient methods are needed to achieve spectrum recovery. The spectrum restoration can be converted to estimating the complex amplitudes and attenuation coefﬁcient in the model of the sum of complex exponentials. The matrix pencil method is a standard technique for mode frequency identiﬁcation. Therefore, we propose the sFFT5.0 algorithm and sFFT-DT3.0 algorithm using the matrix pencil method to do spectrum recovery. These two algorithms are low computational complexity and strong robustness and have achieved good results in the actual comparative test.


Introduction
As we all know, the best method of the Discrete Fourier Transform (DFT) operation is the famous Fast Fourier Transform (FFT) [1], but the complexity of FFT is O(NlogN), which can not meet the needs of some real-time applications and big data requirments. Therefore, sFFT algorithms are proposed to recover only the large value frequency in the sparse signal by using the sparsity feature of the signal. Once proposed, the new algorithm has attracted wide attention and has been continuously developed. In 2012, it was rated as one of the top 10 breakthrough technologies of the year by MIT Technology Review.
The implementation of the sparse Fourier algorithm mainly consists of two steps: bucketization and spectrum recovery. According to different filters used in bucketization, the algorithms can be divided into three categories [2][3][4][5]. The first kind of sFFT algorithm is the randomized algorithm based on the Dirichlet kernel filter bank. The most typical algorithm is the class Ann Arbor Fast Fourier Transform (AAFFT) [6][7][8] algorithm proposed the two new algorithms and other sFFT algorithms are tested. It is proved that the new two algorithms are more efficient and robust.
Through the research of this paper, we can get the advantages and novelty of the new two sFFT algorithms compared with the old algorithms. In theory, the time complexity of the new algorithms is independent of N, so they are very suitable for processing large length signals. In addition, they have good robustness due to the use of the new matrix pencil method. These theoretical conclusions have been verified by experiments, and their general performance table has been obtained through this paper, which provides a reference for the use and research of these two new excellent sFFT algorithms.

Preliminaries
In this section, we introduce the preliminary knowledge and basic definitions of sparse Fourier transform.
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: The DFT of a vector x ∈ C N (consider a signal of size N) is a vectorx ∈ C N defined as follows:x If the spectrumx has exactly K non-zero frequency coefficients while the remaining N − K coefficients are zero, the signal is exactly K-sparse. If the spectrumx has K significant frequency coefficients while the remaining N − K coefficients are approximately equal to zero, the signal is general K-sparse. The goal of sFFT is to recover a K-sparse approximationx by locating K frequency positions f 0 , . . . , f K−1 and estimating K largest frequency coefficientsx f 0 , . . . ,x f K−1 .

Techniques
In this section, we introduce five technologies related to the new sFFT algorithm. They are frequency bucketization by the flat window filter, frequency bucketization by the aliasing filter, estimating the number of effective signals of the "sum of complex exponentials" signal model, the traditional Prony method used to calculate the poles of the model, and the new MPP method used to calculate the poles of the model.

Frequency Bucketization by the Flat Window Filter
The first stage of sFFT is encoding by frequency bucketization. The process of frequency bucketization using the flat window filter is achieved through shift operation, scaling operation, flat filter operation, and subsampling operation.
The shift operation representing the original signal multiplied by matrix S τ . S τ ∈ R N×N (the offset parameter is denoted by τ ∈ R) is defined as Equation (3). The vector x = S τ x, such that x i = x (i−τ) . The time scaling operation representing the original signal multiplied by matrix P σ ∈ R N×N (the scaling parameter is denoted by σ ∈ R) is defined as Equation (4). A vector x = P σ x, such that x i = x σi . Let matrix Q L ∈ C N×N be a diagonal matrix whose diagonal entries represent flat window filter coefficients (filter G ∈ C N be an (L/N, L/2N, δ, w) flat window) in the time domain, defined as Equation (5). The last operation is the frequency subsampled operation, and the signal in the time domain is aliased such that the corresponding signal in the frequency domain is subsampled. Let matrix U L ∈ R B×N represent the aliasing operator as Equation (6): Through the Introduction, we have learned about the above four operations. After the above four operations, the signal can be divided into buckets based on the flat window filter. The specific process is as follows:ŷ L,τ,σ = F B U L Q L S τ P σ x. The results of bucketization in bucket i are shown in Equation (7). In addition, the specific contents can be seen in paper [3]:

Frequency Bucketization by Aliasing Filter
The process of frequency bucketization using the aliasing filter is achieved through shift operation and aliasing filter operation.
The first operation time shift operation is similar to that described in the previous section. The second operation 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. The subsampling factor is denoted by L, representing how many frequencies are aliasing in one bucket. The subsampling number is denoted by B(B = N/L). The subsampling operation represents the original signal multiplied by matrix D L . D L ∈ R B×N is defined as follows: After the above two operations, the signal can be divided into buckets based on the aliasing filter. The specific process is as follows:ŷ B,τ = F B D L S τ x, The results of bucketization in bucket i are shown in Equation (9). In addition, the specific contents can be seen in [4]:ŷ

Estimation of the Number of Effective Signals
This section describes the first step in solving the "sum of complex exponentials" signal model, estimating the number of effective signals.

Problem Statement
In general, the signal model of the observed late time of electromagnetic-energyscattered response from an object can be considered as the "sum of complex exponentials" signal model and be described as follows: In the above formula, m K represents the sampling signal at time k, p j represents a complex coefficient containing amplitude information, z j represents an attenuation index containing phase information and represents an exponential with an amplitude equal to 1. a m means that we assume that the sampled signal is composed of the sum of at most a m effective complex exponential signals. In addition, the sampled signal is actually composed of the sum of a effective complex exponential signals. The approximate equal sign indicates that the signal contains noise in addition to the effective signal. The first problem now is how to estimate the number of effective complex exponential signals.

Problem Solution
According to the maximum possible, the sampling signal is composed of at most a m number of signals, 2a m times is used for sampling. According to the sampled values, a matrix is constructed, SVD decomposition is performed on the matrix, and PCA analysis is performed to obtain the number of effective complex signals.
The specific process is as follows. First, 2a m samples are taken to obtain Equation (11). The matrix M a m is constructed as Equation (12) according to 2a m sampling results. In addition, the vectorẑ j ∈ C a m ×1 defined asẑ j = [z 0 j , z 1 j , . . . , z a m −1 j ] T is composed of each attenuation index: From the above three definitions, it can be concluded that the relationship between M a m , p j ,ẑ j satisfies the following Theorem 1: Proof of Theorem 1.
Based on the properties as mentioned above, we obtain Equation (13).
Equation (13) is similar to the symmetric singular value decomposition (SSVD). Nevertheless, there are some differences. Ref. [27] proved that, for the symmetric matrix, the Σ obtained from the SVD is equal to the Σ obtained from the SSVD. , σ 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 components p j 's. For a matrix that is not very large, the calculation amount of PCA is the third power of the matrix size, so the time complexity is a 3 m , and the sampling complexity is 2a m for the estimation of the number of effective complex exponential signals.

Problem Statement
After knowing the effective number a, the original Equation (12) becomes Equation (14) and the original matrix M a m becomes a new matrix M a . In the next step, it is necessary to estimate each complex coefficient p j and each attenuation index z j according to the m j 's:

Problem Solution by the Prony Method
The Prony method is a method for solving polynomials, so the following polynomial of order a is defined as follows: The Prony method has three steps. The first step is to approximate the coefficient C of the polynomial through the formula, the second step is to approximate z j through the above approximate polynomial, and the third step is to calculate p j .
Let the orthogonal polynomial formula P(z) be defined as Equation (16) and P(z) ≈ 0. Let Matrix M a ∈ C a×a be defined as Equation (15). 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 2a−1 ] T . The moments' formula satisfies Theorem 2. Through Theorem 2, we can obtain C ≈ (M −1 a )M s . Thus, the first step of the Prony method has been completed.
The first element of M s has been proven and other elements of M s can also be proven.
The second step of the Prony method is to calculate the coefficients of the polynomial and find the approximate root of the polynomial. When the polynomial is an equation, directly use the method of solving the equation, so the a number of roots of P(z) = 0 is the solution of z j 's, such as if a = 1, through P(z) = z + c 0 = 0, then z 0 = −c 0 . If a = 2, through If the polynomial is an approximate equation, z j 's can be solved by solving the eigenvalue of the matrix. The third step of the Prony method is to calculate p j 's. Obviously, p j 's can be obtained by the least square method when Equation (14) is known.
In the first step, the complexity of calculate polynomial coefficients C ≈ (M −1 a )M s is a 3 m . In the second step, the complexity of calculate a approximate roots of the polynomial with known coefficients P(z) = z a + c a−1 z a−1 + · · · + c 1 z + c 0 is a 3 m . In the third step, the computational complexity of the least square method is a 3 m . When dealing with noisy data, the Prony method solves the following of polynomials, so the poles obtained are extremely inaccurate. There is a small deviation in resonance frequency and amplitude, about 10 percent to 20 percent, which is very sensitive to noise. A large number of proportional examples are used in the paper [28,29] to prove the low noise immunity of the Prony method.

Calculate the Poles by the MPP Method
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. This method only needs one step to directly convert the pole solving problem into solving the generalized characteristics of the prediction matrix. Ref. [24] points out that this method is efficient and has better accuracy. We first assume that the approximate equation sign becomes the equation sign, that is, there is no noise in the sampled signal.
For our problem, let the Toeplitz matrix Y be defined as Equation (18), let Y 1 be Y with the rightmost column removed and be defined as Equation (19), and let Y 2 be Y with the leftmost column removed and be defined as Equation (20). Let the Vandermondee martix U a+1 ∈ C (a+1)×a be defined as Equation (21): It there is no noise in the sampled signal. so the rank of the matrix Y, Y 1 , Y 2 is equal to a. Then, Equation (22) can be obtained through the Vandermonde decomposi- . Then, Equation (23) can be obtained through the Vandermonde decomposition matrix Y 1 . For ex- . Then, Equation (24) can be obtained through the Vandermonde decomposition matrix Y 2 . For example, a = 2, Through Equations (22)-(24), we can obtain , so the set of generalized eigenvalues of Y 2 − λY 1 are the z j 's we seek. 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 afterward. For details, please refer to paper [24,30]. The third step to calculate p j 's of the MPP method is the least square method similar to the way the Prony method does.
In the first step, the complexity of computing the SVD of the Y, Y = VΣV * is a 3 m . In the second step, the complexity of calculating the set of nonzero 2 ordinary eigenvalues of (Y 1 + )Y 2 is a 3 m . In the third step, the computational complexity of the least square method is a 3 m .

sFFT5.0 Algorithm
The new sFFT5.0 algorithm is a sparse Fourier algorithm based on a flat window filter using the matrix pencil method. The block diagram of this algorithm is shown in Figure 1.

Problem Statement
As shown in Equation (7), the filtered signal of each bucket after the flat window filter is obtained by summing the linear changes of L original signals. In order to recover the spectrum, we need to successfully recover the effective a large value frequency, so Equation (7) is converted to the following equation:

Old Problem Solutions
The old method assumes that there is only one valid signal at most for each bucket of sampling, so the formula can be converted as follows: The sFFT1.0 [9] and sFFT2.0 [9] algorithms use the voting method to recover the spectrum of the large value frequency in the bucket. The method is as follows: ifŷ L,τ,σ [i]'s are large, it means that positions of the effective spectrum may appear in the frequency elements that are hashed by these buckets, then the corresponding position is voted and scored. Correspondingly, ifŷ L,τ,σ [i]'s are small, it means that positions of the effective spectrum do not exist in the frequency elements that are hashed by these buckets, then the corresponding position is voted and without scoring. After multiple rounds of voting, the position with more vote scores is judged as the position of the effective frequency point. This method is relatively simple. However, due to the more rounds of voting, the operation efficiency will be affected. Moreover, this algorithm is a statistical algorithm, a probabilistic algorithm, not a deterministic algorithm.
The sFFT3.0 [10] and sFFT4.0 [10] algorithms use the phase offset method to recover the spectrum of the large value frequency in the bucket. The method is as follows: for the exactly sparse case, in the first round, we set τ 1 = 0, and, in the second round, we set τ 2 = 1; then, suppose the bucket i contains only one large frequency, so we obtain y L,1,σ [i] so we can locate the position f = (σ −1 u)modN. In a general sparse case, we can estimate the possible range of u, narrow the range through multiple iterations, such as binary search or multi-scale search, and finally determine u and f = (σ −1 u)modN. This method has three limitations. One limitation is that the efficiency is very low. The scope needs to be continuously reduced each time in one round, and a total of many rounds are required. The second limitation is that the search has an assumption that only at most one large value in the bucket can be estimated correctly. The third limitation is the general robustness of the algorithm.

New Problem Solution
The sFFT5.0 algorithm, as shown in Equation (23), can convert the spectrum problem in one bucket into the solution of the "sum of complex exponentials" signal model. In the solution process, for specified bucket i, the number of large value frequencies is defined as a, the maximum number of large value frequencies is defined as a m ,Ĝ iL−u jxσ −1 u j = p j represents a complex coefficient, ω u j N = z j represents attenuation index, τ = k represents time coefficient,ŷ L,τ,σ [i] = m k represents sampling signal, and u j ∈ [ (2i−1)L 2 , (2i+1)L 2 − 1] represents the possible range. After such a transformation, the original problem of the spectrum restoration in a bucket can be transformed into the problem of solving the "sum of complex exponentials" signal model.
In the actual solution, we assume that there are at most two coefficients in one bucket that means a m = 2. The detail process is as follows: (1) It is known that the maximum number of effective signals is a m = 2. Calculate the number of effective signals a by SVD decomposition and PCA analysis. (2) After the number of signals a is obtained, z j 's are obtained by using the matrix pencil method. (3) p j 's are obtained by using the least square method. After a effective signals are obtained, it can be concluded that if other signals except these a effective signals still have a lot of energy, that is, the maximum number of effective signals representing this bucket is greater than two. Therefore, the bucket division is inconsistent with our preset assumption, so the calculation results are not adopted. These effective signals will be obtained in the case of bucket division in subsequent iterations. (4) ThroughĜ iL−u jxσ −1 u j = p j , ω τu j N = z j , f = (σ −1 u)modN, obtain every f j andx f j in one bucket. After the spectrum recovery of all buckets is processed, the spectrum recovery of the signal is completed as well. Secondly, we analyze the robustness: from the above analysis, the old method assumes that the sampling of each bucket has at most only one effective signal, so it can not deal with aliasing. The new method can handle the case that there are two effective signals in the bucket, and the matrix pencil method ensures that the position and amplitude of one or two large frequency values can be normally estimated even if the bucket length L is large, so the robustness of the new method is better.

sFFT-DT3.0 Algorithm
The new sFFT-DT3.0 algorithm is a sparse Fourier algorithm based on the aliasing filter using the matrix pencil method. The block diagram of this algorithm is shown in Figure 1.

Problem Statement
As shown in Equation (9), the filtered signal of each bucket after the aliasing window filter is obtained by summing the linear changes of L original signals. In order to recover the spectrum, we need to successfully recover the effective a large value frequency, so Equation (9) is converted to the following equation:

Old Problem Solutions
In the solution process, for a specified bucket i, the number of large value frequencies is defined as a, the maximum number of large value frequencies is defined as a m ,x f j = p j represents a complex coefficient, ω f j N = z j represents an attenuation index, τ = k represents time coefficient,ŷ B,τ [i] = m k represents a sampling signal, f j ∈ {i, B + i, 2B + i, . . . , and(L − 1)B + i} represents the possible range. After such a transformation, the original problem of the spectrum restoration in a bucket can be transformed into the problem of solving the "sum of complex exponentials" signal model.
After the number of effective signals a is obtained, the polynomial Prony method is used to continue the solution in the old algorithm. The sFFT-DT1.0 [14,15] algorithm defaults that the polynomial is an equation and uses the polynomial root method to solve it. This is only suitable for the case without noise. The sFFT-DT2.0 [14,15] algorithm tries every possible position z j = ω (16) P(z) = z a + c a−1 z a−1 + · · · + c 1 z + c 0 , the first a number of the smallest |P(z)| is the solution of z j 's. L attempts are needed in one solution.

New Problem Solution
In the actual solution, we assume that there are at most four coefficients in one bucket that means a m = 4. The detailed process is as follows: (1) It is known that the maximum number of effective signals is a m = 4. Calculate the number of effective signals a by SVD decomposition and PCA analysis. (2) After the number of signals a is obtained, z j 's are obtained by using the matrix pencil method. (3) p j 's are obtained by using the least square method. After a effective signals obtained, it can be concluded that, if other signals except these a effective signals still have a lot of energy, that is, the maximum number of effective signals representing this bucket is greater than four. Therefore, the bucket division is inconsistent with our preset assumption, so the calculation results are not adopted. These effective signals will be obtained in the case of bucket division in subsequent iterations. (4) Throughx f j = p j , ω f j N = z j , get every f j andx f j into one bucket. After the spectrum recovery of all buckets is processed, the spectrum recovery of the signal is completed as well.

Performance Analysis in Theory
The sFFT-DT1.0 algorithm is only suitable for ideal signals by using the polynomial root method, and is not suitable for general situations. The sFFT-DT2.0 algorithm uses the enumeration method. The time complexity of calculating a bucket is O(L + a m ), which is relatively high and has poor robustness. The new sFFT-DT3.0 algorithm uses the matrix pencil method. The time complexity of calculating a bucket is O(a m ), which realizes the optimization of performance and is more robust than the traditional Prony method.

Algorithms Analysis in Practice
In this section, we evaluate the performance of five sFFT algorithms using the flat filter: the sFFT1.0, sFFT2.0, sFFT3.0, sFFT4.0, and sFFT5.0 algorithms. In addition, we evaluate the performance of five sFFT algorithms using the aliasing filter: the sFFT-DT1.0, sFFT-DT2.0, sFFT-DT3.0, R-FFAST, and DSFFT algorithms as well. These old algorithms are obtained through open source, and the new algorithms are implemented according to the design. There are three main test items, namely, the complexity of the algorithm under different lengths, the complexity of the algorithm under different sparsity, and the L 1 error of the algorithms. All experiments were run on a Linux CentOS computer with a 4 Intel(R) Core(TM) i7-4570 3 GHz CPU and 16 GB of RAM.

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. Each test record contains the run time and the L 0 , L 1 , L 2 error between the calculation result and the best result through the algorithm library. Detailed reports can be found on the github website (https://github.com/zkjiang/-/tree/master/docs, accessed on 10 June 2022).

Comparison Experiments of Two New Algorithms Using Different Parameters
In order to better compare and verify the algorithms, this section tests the different parameters of the two proposed algorithms, and determines the selected parameters through the experimental results. Since both algorithms use the matrix pencil method, the most important parameter of the algorithms is the assumed number of maximum effective signals a m after bucketization. It is obvious that, when a m is large, the rank of the matrix to be processed is large, so the time complexity will be greatly increased, but, at the same time, the robustness will be improved because it considers the worst possible condition. From Figures 2 and 3, based on the balance of the two performances, the parameter used in our proposed sFFT5.0 algorithm is a m = 2, and the parameter used in sFFT-DT3.0 algorithm is a m = 4.

Comparison Experiment of sFFT5.0 and Other Algorithms Using the Flat Filter
We plot Figure 4a,b representing the run time vs. signal size and vs. signal sparsity for the sFFT1.0, sFFT2.0, sFFT4.0, and sFFT5.0 algorithms in the general sparse case (The sFFT3.0 algorithm is ignored because it is not suitable for a general sparse case; it is only suitable for an exactly sparse case). Thus, from Figure 4, we can see the following: (1) The run time of these four algorithms were approximately linear in the log scale as a function of N and nonlinear as a function of K. (2) When N is large, the result of ranking run time of the four algorithms is sFFT5.0 > sFFT2.0 > sFFT4.0 > sFFT1.0. From this ranking, it can be seen that, when N is large, the sFFT5.0 algorithm using the matrix pencil method is the best because the complexity of the method is independent of N, and other algorithms are related to N. (3) When K is large, the result of ranking run time of the four algorithms is sFFT4.0 > sFFT2.0 > sFFT1.0 > sFFT5.0. From this ranking, we can see that, when k is large, the sFFT5.0 algorithm using the matrix pencil method has no advantage. Once there are many buckets to be processed, the matrix transformation involved in each bucket is time-consuming. In addition, through the intersection of these curves, we can also find that, with the increase of sparsity, the time complexity of sFFT5.0 algorithm is better than other algorithms when N is greater than or equal to 2 22 and, with the increase of sparsity, the time complexity of sFFT4.0 algorithm is always better than other algorithms when K is greater than or equal to 100.   3) The result of ranking robust of the four algorithms is sFFT5.0 > sFFT2.0 > sFFT1.0 > sFFT4.0. The reason is that, in the positioning methods, the matrix pencil method is better than the voting method, and the voting method is better than the phase method, which leads to the advantages and disadvantages of the algorithms using these methods.

Comparison Experiment of sFFT-DT3.0 and Other Algorithms Using the Aliasing Filter
We plot Figure 6a,b (The DSFFT algorithm under the condition of large sparsity is ignored because it cannot deal with large sparse data) representing the run time vs. signal size and vs. signal sparsity for the sFFT-DT2.0, sFFT-DT3.0, DSFFT, and R-FFAST algorithm in the general sparse case (The FFAST algorithm and sFFT-DT1.0 algorithm are ignored because they are not suitable for a general sparse case; it is only suitable for an exactly sparse case). Thus, from Figure 6, we can see the following: (1) The run time of these four algorithms were approximately linear in the log scale as a function of N and nonlinear as a function of K. (2) The result of ranking run time of the four algorithms is sFFT-DT3.0 > sFFT-DT2.0 > DSFFT > R-FFAST. (3) From this ranking, it can be seen that, in terms of time complexity, the sFFT-DT algorithm using the "sum of complex exponentials" model is the best. Other FFAST algorithms using the peeling framework and DSFFT algorithm using a binary tree framework are more complex and time-consuming. (4) Among the sFFT-DT algorithms, the sFFT-DT3.0 algorithm using the matrix pencil method is better than the sFFT-DT2.0 algorithm using the Prony method because it requires less matrix processing times, and it is independent of N, while other algorithms are related to N. We plot Figure 7a,b representing the run time and L 1 -error vs. SNR. for the sFFT-DT2.0, sFFT-DT3.0, and R-FFAST algorithm in the general sparse case. Thus, from Figure 7, we can see the following: (1) The runtime is approximately equal vs. SNR except the R-FFAST algorithm. (2) To a certain extent, these three algorithms are all robust. (3) The result of ranking robust of the four algorithms is R-FFAST > sFFT-DT3.0 > sFFT-DT2.0. From this ranking, it can be seen that the R-FFAST algorithm using the peeling framework and MMSE estimation method has the best robustness, but among the sFFT-DT algorithms using the "sum of complex exponentials" model, the sFFT-DT3.0 algorithm is better than the sFFT-DT2.0 algorithm using the Prony method, which is determined by the characteristics of the method used.

Summary of the Performance of the Proposed Algorithms and the Related Algorithms
Through theoretical analysis and experimental results, we get a summary table as shown in Table 1.

Conclusions
At the beginning, the paper that introduces the filtered signal is approximately equal to the linear addition of signals in the bucket after aliasing filters or after flat window filters. Based on this situation, the problem of frequency recovery of large values in the bucket is transformed into the problem of solving the amplitude and phase of the "sum of complex exponentials" model. This paper then introduces two methods to calculate the amplitude and phase of the model, the Prony method and MPP method. After comparison, it can be concluded that the MPP method has higher efficiency and stronger robustness.
Based on these, this paper proposes the sFFT5.0 algorithm based on the flat filter and sFFT-DT3.0 algorithm based on the aliasing filter using the matrix pencil method. Compared with the previous sFFT1.0 and sFFT2.0 algorithms using the voting method, the sFFT5.0 algorithm reduces rounds and improves the algorithm efficiency. Compared with the previous sFFT3.0 and sFFT4.0 algorithms using the phase search method, the sFFT5.0 algorithm also reduces rounds, improves the algorithm efficiency, and improves the robustness of the algorithm. Compared with the former sFFT-DT1.0 algorithm using the Prony method, the sFFT-DT3.0 algorithm can process general sparse signals. Compared with the former sFFT-DT2.0 using the enumeration method, the sFFT-DT3.0 algorithm reduces the number of calculations in the bucket, improves the efficiency of the algorithm, and improves the robustness of the algorithm.
Finally, the time complexity and robustness of the new and old four sparse Fourier algorithms based on flat window filter, and the new and old four sparse Fourier algorithms based on aliasing filter are tested. It can be seen from the test results that the two new algorithms using the matrix pencil method have the best time complexity in the case of large length because the complexity of the method using the matrix pencil method is independent of N, and other algorithms are related to N. However, in the case of large sparsity, the time complexity of new algorithms is ordinary because, once there are many buckets to be processed, the matrix transformation involved in each bucket of this method is time-consuming. In addition, in terms of robustness, the two new algorithms are also optimal or suboptimal, which is determined by the anti-interference characteristics of the matrix pencil method. These conclusions also prove the superiority of the new algorithm and lay a foundation for the future use of this algorithm.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. We also made it public on the website https://github.com/zkjiang/-/tree/master/docs. (accessed on 10 June 2022).

Conflicts of Interest:
The authors declare no conflict of interest.