ATSFFT: A Novel Sparse Fast Fourier Transform Enabled With Sparsity Detection

TL;DR

This paper introduces ATSFFT, an adaptive sparse Fourier transform algorithm that detects signal sparsity dynamically, outperforming existing methods and FFT libraries in efficiency and error control.

Contribution

The paper presents ATSFFT, a novel sparse FFT algorithm with sparsity detection and adaptive tuning, enhancing reliability and performance over previous SFFT implementations.

Findings

ATSFFT outperforms traditional SFFT in efficiency and error control.

ATSFFT achieves an order of magnitude faster performance than FFTW.

The method reliably detects sparsity without prior knowledge.

Abstract

The Fast Fourier Transform(FFT) is a classic signal processing algorithm that is utilized in a wide range of applications. For image processing, FFT computes on every pixel's value of an image, regardless of their properties in frequency domain. The Sparse Fast Fourier Transform (SFFT) is an innovative algorithm for discrete Fourier transforms on signals that possess characteristics of the sparsity in frequency domain. A reference implementation of the algorithm has been proven to be efficient than modern FFT library in cases of sufficient sparsity. However, the SFFT implementation has a critical drawback that it only works reliably for very specific input parameters, especially signal sparsity $k$, which hinders the extensive application of SFFT. In this paper, we propose an Adaptive Tuning Sparse Fast Fourier Transform (ATSFFT), which is a novel sparse fast fourier transform enabled with sparsity detection. In the case of unknown sparsity $k$, ATSFFT is capable of probing the sparsity $k$ via adaptive dynamic tuning and completing the sparse Fourier transform. Experimental results show that ATSFFT outperforms SFFT while it is able to control the computation error better than SFFT. Furthermore, ATSFFT achieves an order of magnitude of performance improvement than the state-of-the-art FFT library, FFTW.

Tables

Table 1. TABLE I: The speedup of ASFFT and SFFT over FFTW with fixed sparsity

$k=50$	speedup     \bigstrut
N	$2^8$	$2^9$	$2^{10}$	$2^{11}$	$2^{12}$	$2^{13}$ \bigstrut
ATSFFT/SFFT	$10.9\times$	$5.2\times$	$8.9\times$	$11.4\times$	$17.5\times$	$13.5\times$ \bigstrut
SFFT/FFTW	$0.2\times$	$1.2\times$	$2.3\times$	$3.3\times$	$5.8\times$	$11.6\times$\bigstrut
$k=100$	speedup     \bigstrut
N	$2^8$	$2^9$	$2^{10}$	$2^{11}$	$2^{12}$	$2^{13}$ \bigstrut
ATSFFT/SFFT	$8.3\times$	$11.5\times$	$11.9\times$	$15.3\times$	$13.1\times$	$14.1\times$ \bigstrut
SFFT/FFTW	$0.2\times$	$0.5\times$	$1.7\times$	$2.5\times$	$3.5\times$	$9.3\times$\bigstrut

Table 2. TABLE II: The speedup of ASFFT, SFFT and FFTW with fixed image size

$N=2^{11}$	Speedup     \bigstrut
$k$	$50$	$100$	$200$	$400$	$600$	$800$	$1000$ \bigstrut
ATSFFT/SFFT	$11.4\times$	$15.3\times$	$18.1\times$	$12.6\times$	$11.2\times$	$14.0\times$	$16.6\times$ \bigstrut
SFFT/FFTW	$3.3\times$	$2.5\times$	$1.8\times$	$1.2\times$	$1.1\times$	$0.8\times$	$0.6\times$\bigstrut
$N=2^{12}$	Speedup     \bigstrut
$k$	$50$	$100$	$200$	$400$	$600$	$800$	$1000$ \bigstrut
ATSFFT/SFFT	$17.5\times$	$13.1\times$	$14.2\times$	$16.7\times$	$18.0\times$	$30.3\times$	$32.2\times$ \bigstrut
SFFT/FFTW	$5.8\times$	$3.5\times$	$2.9\times$	$2.3\times$	$1.9\times$	$1.2\times$	$0.9\times$\bigstrut

Table 3. TABLE III: The Error of ASFFT and SFFT

$k=50$	Error     \bigstrut
N	$2^8$	$2^9$	$2^{10}$	$2^{11}$	$2^{12}$	$2^{13}$ \bigstrut
ATSFFT	$5.35\times {10}^{-03}$	$3.15\times {10}^{-05}$	$1.92\times {10}^{-05}$	$5.17\times {10}^{-10}$	$5.38\times {10}^{-10}$	$1.33\times {10}^{-10}$ \bigstrut
SFFT	$8.54\times {10}^{-03}$	$2.27\times {10}^{-04}$	$7.19\times {10}^{-05}$	$5.58\times {10}^{-10}$	$1.44\times {10}^{-09}$	$7.46\times {10}^{-11}$\bigstrut
$k=100$	Error     \bigstrut
N	$2^8$	$2^9$	$2^{10}$	$2^{11}$	$2^{12}$	$2^{13}$ \bigstrut
ATSFFT	$3.66\times {10}^{-04}$	$2.10\times {10}^{-04}$	$3.08\times {10}^{-07}$	$2.08\times {10}^{-07}$	$3.69\times {10}^{-10}$	$2.02\times {10}^{-09}$ \bigstrut
SFFT	$5.46\times {10}^{-04}$	$2.84\times {10}^{-04}$	$6.36\times {10}^{-06}$	$7.94\times {10}^{-07}$	$2.38\times {10}^{-10}$	$2.32\times {10}^{-10}$\bigstrut

Table 4. TABLE IV: The Error of ASFFT, SFFT and FFTW

$N=2^{11}$	Error     \bigstrut
N	$50$	$100$	$200$	$400$	$600$	$800$	$1000$ \bigstrut
ATSFFT	$5.17\times {10}^{-10}$	$2.08\times {10}^{-07}$	$3.43\times {10}^{-05}$	$9.23\times {10}^{-07}$	$8.02\times {10}^{-06}$	$1.45\times {10}^{-05}$	$1.65\times {10}^{-05}$ \bigstrut
SFFT	$5.58\times {10}^{-10}$	$7.94\times {10}^{-07}$	$8.72\times {10}^{-05}$	$1.49\times {10}^{-06}$	$2.57\times {10}^{-05}$	$3.70\times {10}^{-05}$	$3.18\times {10}^{-05}$\bigstrut
$N=2^{12}$	Error     \bigstrut
N	$50$	$100$	$200$	$400$	$600$	$800$	$1000$ \bigstrut
ATSFFT	$5.38\times {10}^{-10}$	$3.69\times {10}^{-10}$	$1.44\times {10}^{-09}$	$2.54\times {10}^{-05}$	$4.73\times {10}^{-05}$	$1.63\times {10}^{-05}$	$1.77\times {10}^{-05}$ \bigstrut
SFFT	$1.44\times {10}^{-9}$	$2.38\times {10}^{-10}$	$1.60\times {10}^{-09}$	$2.64\times {10}^{-06}$	$1.27\times {10}^{-05}$	$1.87\times {10}^{-05}$	$2.25\times {10}^{-05}$\bigstrut

Equations

$${\hat{x}}_{i,j}=\sum _{u\in [N]}\sum _{v\in [N]}{x}_{u,v}{\omega }^{iu+jv},(i,j)\in {\Omega }_N$$

$${\Omega }_N=\{(i,j)|0\le i\le N-1,0\le j\le N-1\}$$

$$(P_{{\sigma }_1,{\sigma }_2,{\tau }_1,{\tau }_2}x{)}_{i,j}=x_{{\sigma }_{1}i+{\tau }_1,{\sigma }_{2}j+{\tau }_2}$$

$${\widehat{(P_{{\sigma }_1,{\sigma }_2,{\tau }_1,{\tau }_2}x)}}_{{\sigma }_{1}i,{\sigma }_{2}j}={\hat{x}}_{i,j}{\omega }^{-({\tau }_{1}i+{\tau }_{2}j)},(i,j)\in {\Omega }_n$$

$${\widehat{(P_{{\sigma }_1,{\sigma }_2,{\tau }_1,{\tau }_2}x)}}_{i,j}=\sum _{u\in [n]}\sum _{v\in [n]}{x}_{{\sigma }_{1}u+{\tau }_1,{\sigma }_{2}v+{\tau }_2}{\omega }^{iu+jv}$$

$${\widehat{(P_{{\sigma }_1,{\sigma }_2,{\tau }_1,{\tau }_2}x)}}_{i,j}$$

$$=\sum _{a_1\in [n]}\sum _{a_2\in [n]}{x}_{a_1,a_2}{\omega }^{\frac{(a_1-{\tau }_1)}{{\sigma }_1}i+\frac{(a_2-{\tau }_2)}{{\sigma }_2}j}$$

$$={\omega }^{-(\frac{{\tau }_1}{{\sigma }_1}i+\frac{{\tau }_2}{{\sigma }_2}j)}\sum _{a_1\in [n]}\sum _{a_2\in [n]}{x}_{a_1,a_2}{\omega }^{(\frac{{\tau }_1}{{\sigma }_1}i+\frac{{\tau }_2}{{\sigma }_2}j)}$$

$$={\omega }^{-({\tau }_1{{\sigma }_1}^{-1}i+{\tau }_2{{\sigma }_2}^{-1}j)}{\hat{x}}_{{{\sigma }_1}^{-1}i,{{\sigma }_2}^{-1}j}$$

$$\displaystyle r(x,y)=\left\{\begin{array}[]{lll}1,&(x,y)\in D\\ 0,&(x,y)\in D^{{}^{\prime}}\end{array}\right.$$

$$f(x,y)=Aexp[-(\frac{(x-x_0{)}^2}{2{\sigma }_x^2}+\frac{(y-y_0{)}^2}{2{\sigma }_y^2})]$$

$$h_{{\sigma }_1,{\sigma }_2}(i,j)=round({\sigma }_1{\sigma }_{2}ij\frac{B^2}{N^2}),i\in [N],j\in [N]$$

$${\hat{x}}_{i,j}^{{}^{\prime }}={\hat{Z}}_{h_{{\sigma }_1,{\sigma }_2}(i,j)}{\omega }^{{\tau }_{1}i+{\tau }_{2}j}/{\hat{G}}_{o_{{\sigma }_1,{\sigma }_2}(i,j)}$$

$$r_i=|1-\frac{k_i}{k_{i-1}}|.$$

$$\displaystyle B_{i+1}=\left\{\begin{array}[]{lll}{\varepsilon_{1}}B_{i},&0{\leq}{r_{i}}<{\delta}_{1}\\ B_{i},&{\delta}_{1}{\leq}{r_{i}}<{\delta}_{2}\\ (1+{\varepsilon_{2}})B_{i},&{\delta}_{2}{\leq}{r_{i}}\\ \end{array}\right.$$

$$\displaystyle B_{i+1}=\left\{\begin{array}[]{lll}{1/2}B_{i},&0{\leq}{r_{i}}<{2\%}\\ B_{i},&{2\%}{\leq}{r_{i}}<{5\%}\\ 2B_{i},&{5\%}{\leq}{r_{i}}\\ \end{array}\right.$$

$$Error=\frac{1}{k}\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}|{\hat{x}}_{i,j}^{{}^{\prime }}-{\hat{x}}_{i,j}|$$

Full text

I Introduction

Nowadays, the development of information technology has reached unprecedent level. The fast growing computing power stimulate emerging technologies that promote the development of human society. Magnetic Resonance Imaging (MRI)[1], Light Field Photography [2], Radio Astronomy [3] and etc. are the applications of image processing that have wide impact on health care, technology and science. There is tremendous demand of high efficient signal processing techniques for drasticlly increasing amount of images to process. The discrete Fourier transform (DFT) is one of the most fundamental and important numerical algorithms which plays the central role in signal processing[4][5], communications, and audio/image/video compression [6]. The Fast Fourier Transform (FFT) [7] that computes the DFT of an $n$-size signal in $O(n\log n)$ time greatly simplifies the complexity of the DFT and boosts the performance substantially, thus is utilized by a broad range of applications.

The general algorithms for computing the exact DFT necessitate the time that is at least proportional to its size $n$. However, it is well known that most image signals posses sparsity in frequency domain. That is, the image signals have naturally sparse representations with respect to fixed Fourier basis. This property has been widely used in various applications including High Efficiency Video Coding (HEVC)[8][9][10][11], compressed sensing [12][13][14] and radio astronomy[reference]. Therefore, for sparse image signals, the lower bound $\Omega(n)$ of the DFT complexity no longer applies. It is crucial to study the new strategy of the Fourier transform based on image sparsity. In 2012, Hassanieh et al proposed one-dimensional sparse fast Fourier transform [15] [16] which is faster than traditional FFT, needless to say, the algorithm demonstrates a promising approach.

However, the SFFT implementation has the drawback that it only works reliably for very specific input parameters, especially signal sparsity $k$. This drawback hinders the extensive applications of SFFT. In addition, two-dimensional sparse Fourier transform can not simply be implemented by utilizing two separate one-dimensional Sparse Fourier transform. Since two-dimensional transform for image signals are more widely used in practival applications, we propose a new two-dimensional Fourier transform that takes advantage of images sparsity(2D-SFFT) [17][18]. Furthermore, we propose an Adaptive Tuning Sparse Fast Fourier Transform (ATSFFT) to improve the accuracy and robustness of SFFT. With adaptive tuning, ATSFFT is able to probe the sparsity $k$ automatically and obtain the Fourier coefficients of signals. Experimental results show that ATSFFT not only can control the error better than SFFT, it could outperform SFFT with an order of magnitude of speedup in some cases.

The remainder of this paper is organized as follows. Section II presents details of the proposed algorithm of two-dimensional Sparse Fourier Transform (SFFT). Section III describes the Adaptive Tuning Sparse Fast Fourier Transform (ATSFFT) algorithm. Experimental results are shown in Section IV, and at last we conclude the paper in Section II-B.

II Two-dimensional Sparse Fast Fourier Transform

First, we lay out several conventions and notations that are used in this paper. A space-domain image is represented as a tow-dimensional matrix $x\in{C^{N\times N}}$, the Fourier spectrum of the image is represented as $\hat{x}$. We assume that $N$ is a power of 2, the notation $[N]$ is defined as the set $\{0,1,...,N-1\}$, and $[N]\times[N]={[N]}^{2}$ denotes the $N\times N$ grid $\{(i,j):i\in[N],j\in[N]\}$. The image support is denoted by $supp(x)\subseteq[N]\times[N]$. All matrix indices are the calculated modulo of the matrix size, e.g. $x_{i,j}$ of image $x$ is actually $x_{i~{}\mathbf{mod}~{}n,j~{}\mathbf{mod}~{}n}$. A set of matrix elements can be written as a matrix subscripted with a set of indices, for example $x_{I,J}=\{x_{i,j}|i\in I,j\in J\}$. In addition, we assume that $B$ is a power of 2, and $N$ can be divisible by $B$.

We define $\omega=e^{-2\pi{i}/N}$ to be a primitive $N$-th root of unity. In the following sections, we will use the following definition of the 2D-DFT without the constant scaling factor:

[TABLE]

This makes some of the proof easier, but it is not considered relevant in practical implementations.

II-A Hash Function

The 2D-SFFT algorithm firstly constructs and utilizes a hash function to extract useful information of an image. The hash function consists of random spectrum permutation, filtering and subsampling in frequency domain.

II-A1 Random Spectrum Permutation

Normally, we do not have access to the input images’ Fourier spectrum since it would involve performing DFT. The spectrum permutation is the primary component of the 2D-SFFT, which is defined in Definition 1. It aims to tear apart the nearby coefficients to reorder the image’s frequency-domain $\hat{x}$:

Definition 1

Let $\sigma_{1}$ and $\sigma_{2}$ be invertible modulo $n$, i.e. $\gcd(\sigma_{1},n)=1$, $\gcd(\sigma_{2},n)=1$, and $\tau_{1}\in[n]$, $\tau_{2}\in[n]$. Then, $i\rightarrow\sigma_{1}i+\tau_{1}~{}\mathrm{mod}~{}n$ and $j\rightarrow\sigma_{2}j+\tau_{2}~{}\mathrm{mod}~{}n$ are permutations on $[n]$. The associated permutation $P_{\sigma_{1},\sigma_{2},\tau_{1},\tau_{2}}$ on a matrix x is then given by

[TABLE]

When a permutation is applied to an image $x$ in space domain, the image’s frequency domain $\hat{x}$ is also permuted. This interesting property is derived in Lemma 1.

Lemma 1

Let $P_{\sigma_{1},\sigma_{2},\tau_{1},\tau_{2}}$ be a permutation and $x$ be an two-dimensional vector. Then

[TABLE]

$Proof.$ For $(i,j)\in{\Omega_{n}}$,

[TABLE]

with $a_{1}=\sigma_{1}u+\tau_{1}$, $a_{2}=\sigma_{2}v+\tau_{2}$

[TABLE]

The Lemma follows by substituting $i={\sigma_{1}}i$, $j={\sigma_{2}}j$. Note that $\omega^{-({\tau_{1}}i+{\tau_{2}}j)}$ changes the phase, but does not change the magnitude of $\hat{x}_{i,j}$.

The permutation in the 2D-SFFT algorithm allows to permute the image’s Fourier spectrum by modifying the image’s space-domain $x$.

II-A2 Window Function

In order to achieve substantial performance improvement, the 2D-SFFT only uses partial input image for computation. The standard window function acts like a filter, it supplies the sparse Fourier transform algorithm with a subset of the Fourier coefficients. Ideally, however, we would like the pass region of the filter to be as flat as possible to avoid spectral leakage. Specifically, two-dimensional flat Guassian window functions are used in 2D-SFFT.

The two-dimensional flat Guassian window function can be obtained from a 2D Gaussian standard window function which is shown in Figure 1 by convolving it with a two-dimensional ”box car” window function which can be presented as:

[TABLE]

where $D=\{(x,y)|-\frac{b}{2}\leq x\leq\frac{b}{2},-\frac{b}{2}\leq y\leq\frac{b}{2}\}$.

The 2D Gaussian window function is defined as

[TABLE]

By applying convolution of (8) and (9), we have the 2D Gaussian flat window function $G$. Figure 2 shows the function in time domain and frequency domain.

Using 2D Gaussian flat window function $G$, part of size $|supp(G)|$ can be extracted out from $P_{\sigma_{1},\sigma_{2},\tau_{1},\tau_{2}}x$ by multiplying $G$ with $x$ and neglecting the coefficients with value of zero. According to the convolution theorem, the multiplication is equivalent to a convolution of $\hat{G}$ and $\hat{x}$. The filtering process can expand the area of non-zero coefficients. This step is to prepare for the subsequent sub-sampling and reverse steps, and further increase the probability of detection of non-zero coefficients.

II-A3 Fast Subsampling and DFT

Lemma 2

Let $B\in N$ divide $n$, $x$ be an $N\times N$ two-dimensional matrix and $y$ be a $B\times B$ two-dimensional matrix with $y_{i,j}=\sum_{u\in[n]}\sum_{v\in[n]}x_{i+Bu,j+Bv}$ for $i=1,...,B$, $j=1,...,B$. Then, $\hat{y}_{i,j}=\hat{x}_{i(n/B),j(n/B)}$

Lemma 2 effectively reduces dimension by subsampling image in time domain and summing up the result. Since image is sparse in frequency domain, dimension reduction can reduce the complexity of position searching and amplitudes of non-zero elements.

An example of an image in frequency domain with sparsity $(k=2)$ is show in Figure 3(a). The process of random spectrum permutation, filtering and subsampling are shown in Figure 3. Permutation can separate nearby coefficients so that the non-zero coefficients can be approximately uniform distributed. Filtering process can expand the area of non-zero coefficients to increase the detection probability. Subsampling effectively reduces complexity.

Random spectrum permutation, filtering and subsampling define the hash function $h_{\sigma_{1},\sigma_{2}}:[N\times N]\to[B\times B]$

[TABLE]

Hash function $h_{\sigma_{1},\sigma_{2}}$ maps each of the $N\times N$ coordinates of the input image to one of $B\times B$ bins.

II-B Two-dimensional Sparse Fourier Transform Algorithm

The 2D-SFFT consists multiple executions of two kinds of operations: seek location and estimate coefficient. The seek location operation is to generate a list of candidate coordinates which have certain probability of being indices of non-zero coefficients in frequency domain. While the estimate coefficient operation is used to exactly determine the frequency coefficients. The implementation of the estimate coefficient also uses hash function. Given a set of coordinates $I$, $\hat{x}_{i,j}$ can be estimated by

[TABLE]

which basically removes the phase change due to the permutation and the filtering.

A simplified workflow diagram of 2D-SFFT is shown in Figure 5. After running multiple iterations of the seek location operation we only keep coordinates emerge in at least half of the seek location loops. For the coordinates $I^{{}^{\prime}}$, the median of the corresponding outputs of the $L$ rounds of the estimate coefficient operation is set to be the frequency coefficient.

III Adaptive Tuning Sparse Fast Fourier Transform

The drawback of the SFFT algorithm is that it only works reliably with prior knowledge of signal sparsity $k$. This drawback prevent it to be utilized by the wide range of applications. Although it is well known that most signals possess the characteristics of sparsity, it is almost impossible to foreknow the sparsity $k$ of the signals. Therefore, we propose a innovative Adaptive Tuning Sparse Fast Fourier Transform (ATSFFT) algorithm. By adaptive dynamic tuning, the ATSFFT can predict the sparsity $k$ and execute the Fourier transform. Experimental results show that the ATSFFT outperforms the SFFT, it also can have a better control of the error.

III-A Adaptive Tuning Iteration

The basic idea behind hash function is to hash the $N\times N$ coefficients of the input image into a small number of $B\times B$ bins. In SFFT, $B$ is determined by $k$, which is set to $\sqrt{Nk}$. From these bins, SFFT selects and only keeps coordinates of top $dk$ points with the largest magnitudes. The actual locations of the Fourier coefficients in the frequency domain are approximated based on these coordinates. We find that after the execution of the hash function, the number of local maximum points is the same with the number of original image’s non-zero coefficients in frequency domain. Specifically, the hash function hashes $2^{11}\times{2^{11}}$ points of the image shown in Figure 3(a) into $2^{6}\times{2^{6}}$ bins. The number of local maximum points is $2$ which is the same with the image sparsity($k=2$). Therefore, by finding the local maximum points, we can avoid prior knowledge sparsity $k$.

ATSFFT initially sets a small value to $B_{0}$ and $k_{0}$. After hashing, we can get $k_{1}$ by counting the number of local maximum values in $B\times B$ matrix. Next, $B_{1}$ can be updated by analyzing the relationship between $k_{0}$ and $k_{1}$. For the convenience of description, the change ratio factor $r$ of $k$ is defined as follows,

[TABLE]

The iterative process is shown as follows

[TABLE]

where $0<{\delta_{1}}<{\delta_{2}}<1$ and $0<{\varepsilon_{1}}<{\varepsilon_{2}}<1$. When $k_{i}$ approaches $k_{i-1}$ within a very small range of ${\delta}_{1}$, we consider $B_{i}$ adapts to the full scatter point set where non-zero elements can adequately distributed, so $B_{i+1}$ need to be decreased to find the better size. When $k_{i}$ approaches $k_{i-1}$ within a appropriate range between ${\delta}_{1}$ and ${\delta}_{2}$, we consider $B_{i}$ adapts to the appropriate scatter point set. If we further reduce the size of $B$, it will lead to large area overlap. Therefore, $B_{i+1}$ is kept the same with $B_{i}$. When $k_{i}$ fluctuates within a large range of $k_{i-1}$, we consider $B_{i}$ is too small to lead large area overlap, so we need to continue to increase the value of $B$. Since $B$ is a power of 2, and $N$ can be divisible by $B$. ATSFFT sets ${\varepsilon_{1}}=1/2$, ${\varepsilon_{2}}=1$, ${\delta_{1}}=2\%$, ${\delta_{2}}=5\%$.

[TABLE]

III-B Adaptive Tuning Sparse Fast Fourier Transform Algorithm

A simplified workflow diagram of ATSFFT is shown in Figure 6. Similar to the SFFT, ATSFFT consists two kinds of operations: seek location and estimate coefficient. In the operation of seek location, by running multiple adaptive tuning iterations, it is able to identify the sparsity $k$ and locate candidate coordinates with a high probability of being one of the $k$ non-zero coordinates. Given the set of candidate coordinates , the ATSFFT can use coefficient estimation to precisely determine the frequency coefficients, which basically remove the phase change due to the permutation and the effect of the filtering.

IV Numerical Experiments

In this section, we present numerical experimental results that compare runtime, speedup and error of ATSFFT, SFFT and FFTW, which is a FFT algorithms implementation known to be one of the fastest FFT libraries [19]. For the experiments, $k$ frequencies are randomly selected and assigned with magnitude of 1, the rest frequencies are set to 0. Each data point in the result graphs is the average over 100 runs with different images. All experiments were carried out on the system that is equiped with 8 Intel Xeon E7-8830 2.13GHz CPU total of 64 cores and 1TB memory.

IV-A Runtime

IV-A1 Runtime vs. Signal Size

The sparsity parameter $k$ is fixed to constant $k=50$ and $k=100$, there are 6 different image matrix size from $2^{8}\times 2^{8}$ to $2^{13}\times 2^{13}$. The average runtime are shown in Figure 7. We can see that SFFT runs faster than FFTW in most cases. However, ATSFFT outperforms its competitors by an order of magnitude. It is no surprise to see that the runtime of the three algorithms are approximately linear in the log scale with respect to problem size. Furthermore, the runtime of ATSFFT grows with the smallest slope than the others as the problem size increases.

In Figure 7(a), it shows that SFFT is faster than FFTW while recovering the exact 50 non-zero coefficients for bigger problem size, and $N={2^{9}}$ is the breaking dimension size. When signal size is less than${2^{9}}\times{2^{9}}$, SFFT is slower than FFTW. ATSFFT is the fastest among the three algorithms for all cases. Similarly, Figure 7(b) shows that for signal dimension size $N>{2^{10}}$ SFFT is faster than FFTW while recovering the exact 100 non-zero coefficients. When signal size is less than ${2^{10}}\times{2^{10}}$, SFFT is slower than FFTW. Again, ATSFFT is the fastest for all cases. The results show that the ATSFFT significantly extends the range of applications of sparse FFT where sparse approximation is applicable.

IV-A2 Runtime vs. Number of Non-zero Frequency

The image matrix size is fixed to constant ${2^{11}}\times{2^{11}}$ and ${2^{12}}\times{2^{12}}$ and the runtime of the comparing algorithms for sparsity $k(k=50,100,200,400,600,800,1000)$ are shown in Figure 8. For each value of $k$, the experiment is repeated 100 times. As the sparsity $k$ increases, ATSFFT and SFFT take more time to complete, while the runtime of FFTW which depends on image size but not sparsity $k$ is essentially constant.

Figure 8(a) shows that SFFT runs faster than FFTW when sparsity $k<700$, and SFFT presents disadvantage when sparsity $k>=800$. While the runtime of ATSFFT for all cases keeps under 0.1 second, ATSFFT outperforms SFFT and FFTW. Similarly, Figure 8(b) shows SFFT is faster than FFTW when sparsity $k<900$, and runs slower otherwise. ATSFFT appears to be the fastest where the runtimes for all cases are under 0.1 seconds as well. The above experimental results show that ATSFFT significantly extends the range of applications for which sparse approximation of SFFT is practical.

IV-B Speedup

An important question we would like to address is that, performance wise, how does ATSFFT compare with SFFT and FFT.

IV-B1 Speedup vs. Signal Size

Table I shows the speedup of ATSFFT to SFFT and SFFT to FFTW where the sparsity is set to be $(k=50)$ and $(k=100)$. For sparsity $(k=50)$, except that SFFT runs slower than FFTW for the case of $N={2^{8}}$, SSFT achieves more speedup when image size increases. However, ATSFFT achieves better performance across the board, and by average a magnitude of speedup than SFFT. For sparsity $(k=100)$, similar observation is obtained. In conclusion, ATSFFT exhibits strong performance and stability over its peers.

IV-B2 Speedup vs. Number of Non-zero Frequency

In this section, we compare the performance with regard to different number of non-zero frequencies when the input image size is fixed. Table II shows the speedup of ATSFFT to SFFT and SFFT to FFTW for image matrix size of ${2^{11}}\times{2^{11}}$ and ${2^{12}}\times{2^{12}}$. For both image sizes, we can see that SFFT performs the best when the sparsity is the smallest, and the speedup over FFTW reduces when the sparsity number $K$ increasing. However, ATSFFT could maintain performance superiority more consistently. For image size ${2^{12}}\times{2^{12}}$, when sparsity is 1000, ATSFFT is about 29 times faster than FFTW while SFFT runs slightly slower than FFTW. Again, the ATSFFT demonstrates robustness over sparsity variation.

IV-C Error

We compute the error metric per as the average $L_{1}$ error between the output $k-$sparse approximation $\hat{x}^{{}^{\prime}}$ and the fft of $x$ referred to as $\hat{x}$.

[TABLE]

IV-C1 Error vs. Signal Size

Tab. III shows the error of the compared algorithms ATSFFT and SFFT where the sparsity parameter $k$ is fixed to constant $(k=50)$ and $(k=100)$ and the image matrix size is from $2^{8}\times 2^{8}$ to $2^{13}\times 2^{13}$. Smaller is better. For sparsity $(k=50)$, the error of ATSFFT and SFFT decreased with the image size increasing. The error of ATSFFT is smaller than SFFT. Similar results appear in the experiments for sparsity $(k=100)$. This shows the sparser the image signal, the smaller the error of ATSFFT and SFFT. Morever, the error of ATSFFT is smaller than SFFT, which indicate ATSFFT can control the error better than SFFT.

IV-C2 Error vs. Number of Non-zero Frequency

In this section, we compare the error with regard to different number of non-zero frequencies and different image size. Tab. IV shows the error of ATSFFT and SFFT increased with the image sparsity increasing. The error of ATSFFT is smaller than SFFT. We can draw the conclusion the sparser the image signal, the smaller the error of ATSFFT and SFFT. Morever, ATSFFT can control the error better than SFFT.

V Conclusion

The Sparse Fast Fourier Transform (SFFT) is a novel algorithm for discrete Fourier transforms on signals with the sparsity in frequency domain. However, the SFFT implementation has the drawback that it only works reliably for very specific input parameters, especially signal sparsity $k$. This drawback hinders the extensive applications of SFFT. we propose an Adaptive Tuning Sparse Fast Fourier Transform (ATSFFT), which is a novel sparse fast fourier transform enabled with sparsity detection. In the case of unknown sparsity $k$, ATSFFT can probe the sparsity $k$ via adaptive dynamic tuning technology and then complete the Fourier transform of signal. We present some numerical experiments comparing runtime, speedup and error of ATSFFT, SFFT and FFT in FFTW algorithms library. Experimental results show that ATSFFT not only can control the error better than SFFT, but also performs faster than SFFT, which computes more efficiently than the state-of-the-art FFTW.