A deterministic sparse FFT algorithm for vectors with small support
AA deterministic sparse FFT algorithm for vectors withsmall support
Gerlind Plonka ∗ Katrin Wannenwetsch † June 12, 2018
Abstract
In this paper we consider the special case where a signal x ∈ C N is known to vanishoutside a support interval of length m < N . If the support length m of x or a goodbound of it is a-priori known we derive a sublinear deterministic algorithm to compute x from its discrete Fourier transform (cid:98) x ∈ C N . In case of exact Fourier measurementswe require only O ( m log m ) arithmetical operations. For noisy measurements, wepropose a stable O ( m log N ) algorithm. Key words. discrete Fourier transform, sparse Fourier reconstruction, sublinearsparse FFT
AMS Subject classifications.
It is well-known that FFT algorithms for the computation of the discrete Fouriertransform of length N require O ( N log N ) arithmetical operations. However, if theresulting vector is a-priori known to be sparse, i.e., contains only a small number ofnon-zero components, the question arises, whether we can do this computation in aneven faster way. In this paper we derive a new deterministic sparse FFT algorithmfor signals x ∈ C N that are a-priori known to vanish outside a support interval oflength m < N . Related work.
In recent years many sublinear algorithms for the sparse Fouriertransform have been proposed that focus on (approximately) m -sparse vectors orvectors being norm bounded.Most of the these methods are randomized poly-logarithmic sparse Fourier algo-rithms achieving e.g. a complexity of O ( m log N ) [5] for m -sparse signals, or even O ( m log m ), see e.g. [11, 12]. However, these algorithms succeed to compute the cor-rect result only with constant probability being smaller than 1. Moreover, there isno efficient method to check the correctness of the result. Another drawback is that ∗ University of G¨ottingen, Institute for Numerical and Applied Mathematics, Lotzestr. 16-18, 37083G¨ottingen, Germany. Email: [email protected] † University of G¨ottingen, Institute for Numerical and Applied Mathematics, Lotzestr. 16-18, 37083G¨ottingen, Germany. Email: [email protected] a r X i v : . [ m a t h . NA ] A p r amples need to be drawn randomly, this is significantly more complicated to realizeby hardware. Regarding the numerical results, runtimes start to be more efficientthan standard FFT for m = 50 and N > , see [3, 14] for most of the consideredalgorithms, and for m = 50 and N > for the algorithm in [6]. Another exper-iment shows for fixed N = 2 that several sparse FFT algorithms start to pay offfor m < m ≤ , [14]. Newestalgorithms are even faster, but the probability to provide the exact result decreaseswith larger m . For a nice overview of the techniques of randomized sparse Fouriertransforms we refer to [3].Deterministic sparse FFT algorithms proposed e.g. in [1, 2, 9, 10, 11] are appli-cable to (approximately) sparse or so-called norm bounded signals. Iwen [9, 10] con-sidered sparse m -term approximations of the discrete Fourier transform and achievedan (cid:96) , (cid:96) / √ m error bound with O ( m log N ) operations. Similar runtime is achievedfor the deterministic algorithm in [11], where accessibility not only to the usual signalsamples but also to time-shifted samples of the input function is assumed. The deter-ministic algorithm in [2] for signals being norm bounded by √ m , the evaluation of δ -approximations of all τ -significant frequencies, has polynomial costs in log N, m, /τ and 1 /δ .Finally, we mention sparse FFT algorithms based on Prony’s method that arecompletely deterministic and can be operated with O ( m ) operations (independentof the signal size N ), see e.g. [7, 13], but usually with an unstable numerical behav-ior. A better numerical performance can be achieved by randomization leading to acomplexity m / log N , see [7].Compared to the approaches above, we restrict ourselves to signal vectors possess-ing a small support interval. Vectors with small support appear in different applica-tions as e.g. in X-ray microscopy, where compact support is a frequently used a-prioricondition in phase retrieval, as well as in computer tomography reconstructions. Notations.
First, we fix some notations. For a vector x ∈ C N of length N ,let its discrete Fourier transform be defined by (cid:98) x = F N x , where the Fourier matrix F N ∈ C N × N is given by F N := (cid:16) ω jkN (cid:17) N − j,k =0 , ω N := e − π i N . Let the support length m = | supp x | of x ∈ C N be defined as the minimal integer m for which there exists a µ ∈ { , . . . , N − } such that the components x k of x vanish for all k / ∈ I := { ( µ + r )mod N, r = 0 , . . . , m − } . The index set I is called support interval of x . Obviously, the support length m of x is an upper bound forits sparsity (cid:107) x (cid:107) , i.e., the number of nonzero components of x , since there may bezero components of x inside the support interval I . However, we always have x µ (cid:54) = 0and x µ + m − (cid:54) = 0. Observe that the support length m of a vector x ∈ C N is alwaysuniquely defined while the support interval itself resp. the first support index µ needsnot to be unique. For example, considering the vector x = ( x k ) N − k =0 ∈ C N of evenlength N given by x = x N/ = 1 and x k = 0 for k ∈ { , . . . , N − } \ { , N/ } ,we find for x the support length m = N/ { , . . . , N/ } or as { N/ , . . . , N − , , } . However, if m ≤ N , thenthe support interval and hence also the first support index µ are uniquely determined. ur contribution. In this paper, we will derive new deterministic algorithmsto reconstruct a vector x ∈ C N with support length m < N from its discrete Fouriertransform (cid:98) x ∈ C N . In case of exact data, we will use less than 4 m Fourier samplesand O ( m log m ) arithmetical operations to recover x . Thus, the algorithm alreadystarts to pay off for m < N .In case of noisy Fourier measurements, we will use O ( m log N ) arithmetical op-erations. More exactly, we will employ one or several FFTs of size less than 4 m aswell as the computation of (cid:98) log Nm (cid:99) − m to recover thefull vector in a stable way.In both cases, for exact as well as for noisy measurements, the algorithms arebased on the idea that for N = 2 J the nonzero components of x can already becomputed by evaluating a periodization of x with length 2 L ≥ m . Thus, for thecomplete reconstruction of x , we only need to compute in a second step the correctsupport interval resp. the first support index of x . Organization of the paper.
In Section 2, we recall some properties of thediscrete Fourier transform that will be used in our approach. Section 3 is devoted tothe sparse FFT algorithm for m -sparse vectors in case of exact Fourier measurements.Section 4 considers a more stable variant of the different steps of the algorithm thatmake the method robust in presence of noise, and even improves the accuracy of theresulting vector in comparison to the usual FFT method while using only O ( m log N )arithmetical operations. Finally, we present our numerical experiments in Section 5. Throughout the paper let N := 2 J with some J > x ( j ) ∈ C j of x , x ( j ) = ( x ( j ) k ) j − k =0 = J − j − (cid:88) (cid:96) =0 x k +2 j (cid:96) j − k =0 (2.1)for j = 0 , . . . , J . Obviously, x (0) = (cid:80) N − k =0 x k is the sum of all components of x , x (1) =( (cid:80) N/ − k =0 x k , (cid:80) N/ − k =0 x k +1 ) T and x ( J ) = x . We recall the following relationship forthe discrete Fourier transform of the vectors x ( j ) , j = 0 , . . . , J , in terms of (cid:98) x . Lemma 2.1
For the vectors x ( j ) ∈ C j , j = 0 , . . . , J , in (2.1) , we have the discreteFourier transform (cid:98) x ( j ) := F j x ( j ) = ( (cid:98) x J − j k ) j − k =0 , where (cid:98) x = ( (cid:98) x k ) N − k =0 = F N x is the Fourier transform of x ∈ C N . Proof:
By definition, we find for the components (cid:98) x ( j ) k of (cid:98) x ( j ) (cid:98) x ( j ) k : = j − (cid:88) r =0 x ( j ) r ω rk j = j − (cid:88) r =0 2 J − j − (cid:88) (cid:96) =0 x r +2 j (cid:96) ω J − j rkN = N − (cid:88) n =0 x n ω J − j nkN = (cid:98) x J − j k , k = 0 , . . . , j − . hus the assertion follows.Lemma 2.1 tells us that for a given vector (cid:98) x the Fourier transform of the peri-odized vectors (cid:98) x ( j ) needs not to be computed but is just obtained by picking suitablecomponents of (cid:98) x . Further, we want to make use of the following simple observation. Lemma 2.2
Let x = ( x k ) N − k =0 ∈ C N , N = 2 J , and let y = ( y k ) N − k =0 ∈ C N be itsshifted version with components y k := x ( k +2 j ν )mod N , k = 0 , . . . , N − , for some j ∈ { , . . . , J − } and ν ∈ { , . . . , J − j − } . Then the components of the Fouriertransformed vectors (cid:98) x = F N x , (cid:98) y = F N y satisfy (cid:98) y l = ω − lν J − j (cid:98) x l , l = 0 , . . . , N − . In particular, for j = J − and ν = 1 we have y k = x ( k + N/ N , k = 0 , . . . , N − ,with (cid:98) x k = (cid:98) y k , (cid:98) x k +1 = − (cid:98) y k +1 , k = 0 , . . . , N − . Proof:
With (cid:98) x = F N x = ( (cid:98) x l ) N − l =0 and (cid:98) y = F N y = ( (cid:98) y l ) N − l =0 we obtain (cid:98) y l = N − (cid:88) k =0 y k ω lkN = N − (cid:88) k =0 x ( k +2 j ν )mod N ω lkN = N − (cid:88) k =0 x k ω l ( k − j ν ) N = ω − lν J − j (cid:98) x l . For j = J − ν = 1 the assertion follows with ω − lν J − j = ω − l = ( − l . We assume that the Fourier transformed vector (cid:98) x = F N x is given, and that thesupport length m of x or an upper bound m of it is known a-priori. Now, we applythe following simple procedure to recover x .Let 2 L − < m ≤ L . Then, by Lemma 1, (cid:98) x ( L +1) = ( (cid:98) x J − ( L +1) k ) L +1 − k =0 is theFourier transform of x ( L +1) . In a first step, we compute x ( L +1) using inverse FFT oflength 2 L +1 .Since | supp x | = m ≤ L by assumption, it follows that x ( L +1) has already thesame support length as x since for each k ∈ { , . . . , L +1 − } the sum in x ( L +1) k = J − L − − (cid:88) (cid:96) =0 x k +2 L +1 (cid:96) (3.1)contains at most one nonvanishing term. Moreover, also the support itself, or equiv-alently the first support index µ = µ ( L +1) of x ( L +1) , is uniquely determined.Thus, in order to recover x from x ( L +1) , we only need to compute the correct firstsupport index µ ( J ) of x , such that the components of x are determined by x ( µ ( J ) + k )mod N = (cid:40) x ( L +1)( µ ( L +1) + k )mod 2 L +1 k = 0 , . . . , m − , k = m, . . . , N − . From (3.1) we know that µ ( J ) is of the form µ ( J ) = µ ( L +1) + 2 L +1 ν for some ν ∈{ , , . . . , J − L − − } . heorem 3.1 Let x ∈ C N , N = 2 J , have support length m (or a support lengthbounded by m ) with L − < m ≤ L . For L < J − , let x ( L +1) be the L +1 -periodization of x . Then x can be uniquely recovered from x ( L +1) and one nonzerocomponent of the vector ( (cid:98) x k +1 ) N/ − k =0 . Proof:
Let x ( L +1) have the support interval { ( µ ( L +1) + r )mod 2 L +1 , r = 0 , . . . , m − } .Since m ≤ L , this support interval is uniquely determined. Further, we know thatby construction the first index µ ( J ) of the support interval of x has the form µ ( J ) = µ ( L +1) +2 L +1 ν for some ν ∈ { , . . . , J − L − − } . We now consider the vector u ∈ C N that is given by the components u µ ( L +1) + k = (cid:40) x ( L +1)( µ ( L +1) + k )mod2 L +1 k = 0 , . . . , m − , k = m, . . . , N − . The vector u is obtained as one possible vector in C N with support length m possess-ing the 2 L +1 -periodization x ( L +1) . Obviously the first index of the support intervalof u is µ ( L +1) . Therefore, all further vectors in C N with a support length m andpossessing the 2 L +1 -periodization x ( L +1) can be written as a shifted version of u ofthe form u ν := ( u νk ) N − k =0 with u νk = u k +2 L +1 ν )mod N , k = 0 , . . . , N − , for some ν ∈ { , . . . , J − L − − } . Thus, the desired vector x is contained in the set { u ν : ν = 0 , . . . , J − L − − } .It remains to find ν such that x = u ν . From Lemma 2.2 it follows for the compo-nents of the Fourier transform (cid:98) u ν = ( (cid:98) u νl ) N − l =0 that (cid:98) u νl = ω − lν J − L − (cid:98) u l .Choosing now an odd-indexed nonzero Fourier value (cid:98) x k +1 of the given vector (cid:98) x ,we can compare it to the corresponding component of (cid:98) u , (cid:98) u k +1 = m − (cid:88) r =0 x ( L +1)( µ ( L +1) + r )mod2 L +1 ω ( µ ( L +1) + r )(2 k +1) N and obtain the correct value ν from ω − (2 k +1) ν J − L − = (cid:98) x k +1 / (cid:98) u k +1 . We remark that ν ∈ { , . . . , J − L − − } is indeed uniquely determined by ω − (2 k +1) ν J − L − since gcd(2 k +1 , J − L − ) = 1, i.e., 2 k + 1 and 2 J − L − are mutually prime. Finally, the vector x isobtained as x = u ν .We summarize the algorithm to evaluate x from (cid:98) x for exact Fourier data as follows. Algorithm 3.2 (Sparse FFT for vectors with small support for exact Fourier data)
Input: (cid:98) x ∈ C N , N = 2 J , | supp x | ≤ m < N . • Compute L such that L − < m ≤ L , i.e., L := (cid:100) log m (cid:101) . • If L = J or L = J − , compute x = F − N (cid:98) x using an FFT algorithm of length N . • If L < J − : . Choose (cid:98) a := ( (cid:98) x J − ( L +1) k ) L +1 − k =0 and compute x ( L +1) := F − L +1 (cid:98) a using an FFT algorithm for the inverse discrete Fourier transform of length L +1 .2. Determine the first support index µ ( L +1) ∈ { , . . . , L +1 − } of x ( L +1) such that x ( L +1) µ ( L +1) (cid:54) = 0 and x ( L +1) k = 0 for k / ∈ { ( µ ( L +1) + r )mod 2 L +1 , r =0 , . . . , m − } .3. Choose a Fourier component (cid:98) x k +1 (cid:54) = 0 of (cid:98) x and compute the sum (cid:98) u k +1 := m − (cid:88) (cid:96) =0 x ( L +1)( µ ( L +1) + (cid:96) )mod 2 L +1 ω (2 k +1)( µ ( L +1) + (cid:96) ) N .
4. Compute the quotient b := (cid:98) x k +1 / (cid:98) u k +1 that is by construction of the form b = ω p J − L − for some p ∈ { , . . . , J − L − − } , and find ν ∈ { , . . . , J − L − − } such that (2 k + 1) ν = p mod 2 J − L − .5. Set µ ( J ) := µ ( L +1) + 2 L +1 ν , and x := ( x k ) N − k =0 with entries x ( µ ( J ) + (cid:96) )mod N := (cid:40) x ( L +1)( µ ( L +1) + (cid:96) )mod 2 L +1 (cid:96) = 0 , . . . , m − , (cid:96) = m, . . . , N − . Output: x . We simply observe that the proposed algorithm has an arithmetical complexityof O ( m log m ). In step 1 we employ an FFT algorithm of length 2 L +1 < m havingthis complexity. All other steps can be performed with O ( m ) or less operations.Particularly, the support interval (resp. the first support index µ ( L +1) of x ( L +1) ) instep 2 can e.g. be found by computing the local energies e k := (cid:80) m + k − (cid:96) = k | x ( L +1) (cid:96) mod 2 L +1 | for k = 0 , . . . , L +1 −
1, and taking µ ( L +1) := argmax k e k . Here, e can be computedby at most O ( m ) operations, and all further energies by using the recursion e k +1 = e k − | x k | + | x k + m | , where we need to keep in mind that there are only m nonzeroentries in x ( L +1) .The complete algorithm requires less then 4 m Fourier values for the overall eval-uation.Let us give some further remarks on Algorithm 3.2.
Remark 3.3
It is always possible to choose a nonzero Fourier component of the vector ( (cid:98) x k +1 ) N/ − k =0 . This can be seen as follows. Let supp x = { µ ( J ) , µ ( J ) +1 mod N, . . . , µ ( J ) + m − N } be the support of x , where m ≤ L ≤ N/ . Then the trigonometricpolynomial p ( ω ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N − (cid:88) k =0 x k e − i ωk (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) e − i ωµ ( J ) m − (cid:88) (cid:96) =0 x ( µ ( J ) + (cid:96) ) mod N e − i ω(cid:96) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (3.2) is real, non-negative and of degree ≤ m − . Hence it possesses at most m − pairwise different zeros, where all zeros are at least double zeros. Observing now that | (cid:98) x k | = p (cid:0) πkN (cid:1) , we can conclude that not all (cid:98) x k +1 , k = 0 , . . . , N/ − , can be zerosof p ( ω ) since N/ ≥ m . emark 3.4 As shown in the remark above, there can be at most m − zero components in thevector ( (cid:98) x k +1 ) N/ − k =0 . However, for a stable computation of the correct first supportindex µ ( J ) , the Fourier value (cid:98) x k +1 used in Algorithm should have large modulus.This may be ensured by taking k := argmax {| (cid:98) x k +1 | : k = 0 , . . . , N/ − } . Unfortunately, this procedure is quite costly. Instead, we may compute k ( L +1)0 := argmax {| (cid:98) x J − ( L +1) k | : k = 0 , . . . , L +1 − } . Then we have for the trigonometric polynomial in (3.2) , | (cid:98) x J − ( L +1) k ( L +1)0 | = p (cid:18) πk ( L +1)0 L +1 (cid:19) = max k =0 ,..., L +1 − p (cid:0) πk L +1 (cid:1) > , and it is likely that p (cid:18) πk ( L +1)0 L +1 (cid:19) is close to the global maximum (cid:107) p (cid:107) ∞ of p ( ω ) . Instep 3 of Algorithm we may now choose the one neighboring Fourier value withmaximal modulus, i.e., either (cid:98) x J − ( L +1) k ( L +1)0 +1 or (cid:98) x J − ( L +1) k ( L +1)0 − .If p ( ω ) attains its global maximum at some point ω ∈ (cid:104) π ( k ( L +1)0 − / L +1 , π ( k ( L +1)0 +1 / L +1 (cid:17) ,then it follows that (cid:12)(cid:12)(cid:12) ω − πk ( L +1)0 L +1 (cid:12)(cid:12)(cid:12) ≤ π L +1 . With the above choice of (cid:98) x k +1 with k =2 J − L − k ( L +1)0 or k = 2 J − L − k ( L +1)0 − , we can further assume that (cid:12)(cid:12)(cid:12) ω − π (2 k +1) N (cid:12)(cid:12)(cid:12) ≤ π L +1 , and applying the Lemma of Steˇckin (see e.g. [4]), we find | (cid:98) x k +1 | = p (cid:16) π (2 k +1) N (cid:17) ≥ (cid:107) p (cid:107) ∞ cos (cid:16) π ( m − L +1 (cid:17) > . Let us now assume that the given Fourier data are noisy, i.e., they are perturbed byuniform noise, (cid:98) y k = (cid:98) x k + (cid:15) k (4.1)with | (cid:15) k | ≤ δ . Similarly as before, we want to reconstruct x from the noisy Fouriervector (cid:98) y using the a-priori knowledge that x has a support interval of length m < N .For that purpose, we propose a stabilized variant of Algorithm 3.2. The stabiliza-tion regards the following essential steps in the algorithm, namely(1) the correct determination of the support interval of x ( L +1) ,(2) the correct determination of the support interval resp. the first support index µ ( J ) of the desired vector x , and(3) the evaluation of the nonzero components of x within the support interval thatmay be improved by employing more Fourier values (cid:98) y than in the case of exact data. (1) Stable determination of the support interval of x ( L +1) . For that purpose we use the following observation. heorem 4.1 Let x ∈ C N , N = 2 J , have support length m (or a support lengthbounded by m ) with L − < m ≤ L . For L < J − , let x ( L +1) = ( x ( L +1) k ) L +1 − k =0 be the L +1 -periodization of x and (cid:98) x ( L +1) = ( (cid:98) x J − L − k ) L +1 − k =0 its Fourier transform.Then, for each shifted vector of the form (cid:98) z ( κ ) := (cid:0)(cid:98) x J − L − k + κ (cid:1) L +1 − k =0 , κ = 0 , . . . , J − L − − , the inverse Fourier transform z ( κ ) = ( z ( κ ) (cid:96) ) L +1 − (cid:96) =0 = F − L +1 (cid:98) z ( κ ) satisfies | z ( κ ) (cid:96) | = | x ( L +1) (cid:96) | (cid:96) = 0 , . . . , L +1 − . Proof:
By definition, we obtain for the components z ( κ ) (cid:96) of z ( κ ) z ( κ ) (cid:96) = 12 L +1 2 L +1 − (cid:88) k =0 (cid:98) x J − L − k + κ ω − k(cid:96) L +1 = 12 L +1 2 L +1 − (cid:88) k =0 N − (cid:88) r =0 x r ω r (2 J − L − k + κ ) N ω − k(cid:96) L +1 = 12 L +1 N − (cid:88) r =0 x r ω rκN L +1 − (cid:88) k =0 ω k ( r − (cid:96) )2 L +1 = J − L − − (cid:88) j =0 x (cid:96) +2 L +1 j ω ( (cid:96) +2 L +1 j ) κN = ω (cid:96)κN J − L − − (cid:88) j =0 x (cid:96) +2 L +1 j ω L +1 jκN . Since | supp x | = m < L +1 , for each (cid:96) the above sum contains only one value, thus | z ( κ ) (cid:96) | = | x ( L +1) (cid:96) | .Obviously, we have z (0) = x ( L +1) by definition. We observe that all vectors z ( κ ) are constructed from different Fourier components of (cid:98) x . This circumstance will beused to stabilize the algorithm. Applying the noisy measurements (cid:98) y k in (4.1) insteadof (cid:98) x k , let (cid:98)(cid:101) z ( κ ) := ( (cid:98) y J − L − k + κ ) L +1 − k =0 . As before in Algorithm 3.2, we compute in afirst step the vector (cid:101) z (0) = y ( L +1) from the noisy measurements ( (cid:98) y J − ( L +1) k ) L +1 − k =0 .In order to determine the support interval of x ( L +1) from the vector (cid:101) z (0) we considerthe energies (cid:101) e (0) k := m + k − (cid:88) (cid:96) = k | y ( L +1) (cid:96) mod2 L +1 | = m + k − (cid:88) (cid:96) = k | (cid:101) z (0) (cid:96) mod 2 L +1 | , k = 0 , . . . , L +1 − , and take µ ( L +1)0 := argmax k (cid:101) e (0) k as the first estimate for µ ( L +1) . For higher noiselevels, we stabilize the computation of µ ( L +1) as follows: We compute also the vector (cid:101) z (2 J − L − ) by applying the inverse FFT to (cid:98)(cid:101) z (2 J − L − ) = ( (cid:98) y J − L − (2 k +1) ) L +1 − k =0 , computethe energies (cid:101) e (1) k := m + k − (cid:88) (cid:96) = k | (cid:101) z (2 J − L − ) (cid:96) mod 2 L +1 | , k = 0 , . . . , L +1 − , nd take µ ( L +1)1 := argmax k ( (cid:101) e (0) k + (cid:101) e (1) k ). If µ ( L +1)1 = µ ( L +1)0 , then this index is takenas the first support index of x ( L +1) . Otherwise, we compute a further vector, e.g., (cid:101) z (2 J − L − ) = F − L +1 ( (cid:98) y J − L − (4 k +1) ) L +1 − k =0 , evaluate the corresponding energies (cid:101) e (2) k and µ ( L +1)2 := argmax k ( (cid:101) e (0) k + (cid:101) e (1) k + (cid:101) e (2) k ), etc. (2) Stable determination of the support interval of x. In order to stabilize the computation of the first support index µ ( J ) of the completevector x = x ( J ) resp. the shift ν such that µ ( J ) = µ ( L +1) − L +1 ν , we employ aniterative procedure, where at each iteration step, the support interval of the periodizedvector x ( j +1) , j = L + 1 , . . . J −
1, is computed from the support interval of the vector x ( j ) of half length.At each iteration step we apply the following procedure. Observing that x ( j +1) has the same support length m as x ( j ) , and using the relation x ( j +1) k + x ( j +1) k +2 j = x ( j ) k , k = 0 , . . . , j − , it is sufficient to check whether µ ( j +1) = µ ( j ) or µ ( j +1) = µ ( j ) + 2 j , i.e., whether thefirst support index µ ( j ) of x ( j ) is equal to the first support index µ ( j +1) of x ( j +1) , orif the support has to be shifted by 2 j . We illustrate the two cases in detail in Figure1. First case: µ ( j +1) = µ ( j ) . x ( j ) = ⇒ x ( j +1) x ( j ) = ⇒ x ( j +1) Second case: µ ( j +1) = µ ( j ) + 2 j . x ( j ) = ⇒ x ( j +1) x ( j ) = ⇒ x ( j +1) Figure 1:
Possible support change in one iteration step.
Generally, in order to recover x ( j +1) from x ( j ) , we apply the following procedure. Theorem 4.2
Let x ∈ C N , N = 2 J , have support of length m (or a support lengthbounded by m ) with L − < m ≤ L . Further, let x ( j ) , L + 1 ≤ j ≤ J − , be theperiodizations of x = x ( J ) as given in (2.1) . Then, for each j = L + 1 , . . . , J − , thevector x ( j +1) can be uniquely recovered from x ( j ) and one nonzero component of thevector of Fourier values y := ( (cid:98) x J − ( j +1) (2 k +1) ) j − k =0 . Proof:
Let x ( j ) be the given vector of length 2 j with support supp x ( j ) = { ( µ ( j ) + r ) mod 2 j , r = 0 , . . . , n − } . In order to determine x ( j +1) , we only need to decide hether µ ( j +1) = µ ( j ) or µ ( j +1) = µ ( j ) + 2 j , or equivalently, whether x ( j +1) is givenby x ( j +1) µ ( j ) + (cid:96) = (cid:26) x ( µ ( j ) + (cid:96) ) mod 2 j (cid:96) = 0 , . . . , n − , x ( j +1)( µ ( j ) +2 j + (cid:96) ) mod 2 j +1 = (cid:26) x ( µ ( j ) + (cid:96) ) mod 2 j (cid:96) = 0 , . . . , n − , x ( j +1) only differ by a shift of all components by 2 j . Forsimplicity let us denote the two solutions by u (0) and u (1) . Then, according to Lemma2.2, the Fourier transformed vectors (cid:98) u (0) = ( (cid:98) u (0) k ) j +1 − k =0 and (cid:98) u (1) = ( (cid:98) u (1) k ) j +1 − k =0 satisfythe relationship (cid:98) u (0)2 k +1 = − (cid:98) u (1)2 k +1 , k = 0 , . . . , j − . Using now one nonzero Fourier value (cid:98) x J − ( j +1) (2 k +1) = (cid:98) x ( j +1)2 k +1 , we can determinethe correct vector x ( j +1) . For that purpose, we just compute (cid:98) u (0)2 k +1 = m − (cid:88) (cid:96) =0 x ( j )( µ ( j ) + (cid:96) ) mod 2 j ω (2 k +1)( (cid:96) + µ ( j ) )2 j +1 and compare it to the Fourier value (cid:98) x J − ( j +1) (2 k +1) . If | (cid:98) u (0)2 k +1 − (cid:98) x J − ( j +1) (2 k +1) | < | (cid:98) u (0)2 k +1 + (cid:98) x J − ( j +1) (2 k +1) | , then x ( j +1) = u (0) and µ ( j +1) = µ ( j ) , otherwise, we find x ( j +1) = u (1) and µ ( j +1) = µ ( j ) + 2 j . (3) Evaluation of the nonzero components of x. Finally, we can use the vectors (cid:101) z ( κ ) ∈ C L +1 computed in step (1) also for a moreexact evaluation of the nonzero components x k of x as follows. First, we employ ourinvestigations in Theorem 4.1 and observe that (cid:101) z ( κ ) (cid:96) mod 2 L+1 = ω (cid:96)κN ( x ( (cid:96) +2 L +1 ν ) mod N ω L +1 νκN ) , (cid:96) = µ ( L +1) , . . . , µ ( L +1) + m − . Using the equation µ ( J ) = µ ( L +1) + 2 L +1 ν , we can reformulate this as (cid:101) z ( κ )( µ ( L +1) + k ) mod 2 L+1 = x ( µ ( J ) + k ) mod N ω κ ( µ ( J ) + k ) N , k = 0 , . . . , m − . Therefore, if we have to compute the support entries x ( µ ( J ) + k ) mod N , k = 0 , . . . , m − (cid:98) y k ) N − k =0 of (cid:98) x , we can take the average x ( µ ( J ) + k ) mod N = 1 B + 1 B (cid:88) r =0 (cid:101) z ( κ r )( µ ( L +1) + k ) mod 2 L+1 ω − κ r ( µ ( J ) + k ) N , k = 0 , . . . , m − , where B + 1 is the number of the vectors (cid:101) z ( κ ) = F − L +1 ( (cid:98) y J − L − k + κ ) L +1 − k =0 that wewant to involve.The complete algorithm to compute x from its Fourier transform by a fast sparseFFT algorithm can now be summarized as follows. lgorithm 4.3 (Sparse FFT for vectors with small support for noisy Fourier data) Input: noisy measurement vector (cid:98) y ∈ C N , N = 2 J , | supp x | ≤ m < N . • Compute L such that L − < m ≤ L , i.e., L := (cid:100) log m (cid:101) . • If L = J or L = J − , compute x = F − N (cid:98) y by inverse FFT. • If L < J − :1. Choose (cid:98)(cid:101) z (0) := (cid:98) y ( L +1) = ( (cid:98) y J − ( L +1) k ) L +1 − k =0 and compute (cid:101) z (0) = y ( L +1) := F − L +1 (cid:98) y ( L +1) using an FFT algorithm for the inverse discrete Fourier trans-form.2. Determine the first support index µ ( L +1) ∈ { , . . . , L +1 − } of x ( L +1) usingthe following iteration: • Compute the energies (cid:101) e (0) k := m + k − (cid:88) (cid:96) = k | (cid:101) z (0) (cid:96) mod 2 L +1 | , k = 0 , . . . , L +1 − , and compute µ ( L +1)0 := argmax k (cid:101) e (0) k . • Compute (cid:101) z ( J − L − by IFFT from ( (cid:98) y J − L − (2 k +1) ) L +1 − k =0 , determine (cid:101) e (1) k := m + k − (cid:88) (cid:96) = k | (cid:101) z (2 J − L − ) (cid:96) mod 2 L +1 | , k = 0 , . . . , L +1 − , and take µ ( L +1)1 := argmax k ( (cid:101) e (0) k + (cid:101) e (1) k ) . • Set j := 0 .While µ ( L +1) j (cid:54) = µ ( L +1) j +1 proceed by computing for a further κ ∈ { , . . . , J − L − − } \ { J − L − } the vector (cid:101) z ( κ ) by IFFT from ( (cid:98) y J − L − k + κ ) ) L +1 − k =0 , the energies (cid:101) e ( j +2) k := m + k − (cid:88) (cid:96) = k | (cid:101) z ( κ ) (cid:96) mod 2 L +1 | , k = 0 , . . . , L +1 − , and take µ ( L +1) j +2 := argmax k j +3 (cid:80) j +2 r =0 (cid:101) e ( r ) k .Set µ ( L +1) := µ ( L +1) j +2 and j := j + 1 .End (while).Set x ( L +1) = ( x ( L +1) k ) L +1 − k =0 with x ( L +1)( k + µ ( L +1) ) mod 2 L +1 := (cid:40) (cid:101) z (0)( k + µ ( L +1) ) mod 2 L +1 k = 0 , . . . , m − , k = m, . . . , L +1 − .
3. For j = L + 1 , . . . , J − Choose a Fourier component (cid:98) y J − ( j +1) (2 k +1) = (cid:98) y ( j +1)2 k +1 (cid:54) = 0 and compute a j +1 := m − (cid:88) (cid:96) =0 x ( L +1)( µ ( L +1) + (cid:96) ) mod 2 L +1 ω (2 k +1)( µ ( j ) + (cid:96) )2 j +1 . f | a j +1 − (cid:98) y ( j +1)2 k +1 | < | a j +1 + (cid:98) y ( j +1)2 k +1 | , then set µ ( j +1) := µ ( j ) and x ( j +1) :=( x ( j +1) k ) j +1 − k =0 with entries x ( j +1) µ ( j ) + (cid:96) = (cid:40) x ( j )( µ ( j ) + (cid:96) ) mod 2 j (cid:96) = 0 , . . . , n − , . If | a j +1 − (cid:98) y ( j +1)2 k +1 | ≥ | a j +1 + (cid:98) y ( j +1)2 k +1 | , then set µ ( j +1) := µ ( j ) +2 j and x ( j +1) :=( x ( j +1) k ) j +1 − k =0 with entries x ( j +1)( µ ( j ) +2 j + (cid:96) )mod 2 j +1 = (cid:40) x ( j )( µ ( j ) + (cid:96) )mod 2 j (cid:96) = 0 , . . . , n − , .
4. Assuming that (cid:101) z ( κ r ) ∈ C L +1 , r = 0 , . . . , B , have been evaluated already instep 2, we compute x ( (cid:96) +2 L +1 ν ) mod N = 1 B + 1 B (cid:88) r =0 (cid:101) z ( κ r ) (cid:96) ω − κ r ( (cid:96) +2 L +1 ν ) N ,(cid:96) = µ ( L +1) , . . . , µ ( L +1) + m − . Output: x ( J ) = x . Let us shortly summarize the arithmetical complexity of the algorithm. Step 1and 2 together require two or more inverse FFTs of length 2 L +1 < n , i.e. O ( m log m )arithmetical operations. The number of involved vectors (cid:101) z ( κ ) depends on the noiselevel. The evaluation of energies can be done at each iteration step by O ( m ) opera-tions.In step 3, J − ( L + 1) = log N − (cid:100) log m (cid:101) − < log ( N/m ) iterations are needed,where at each iteration a scalar product of length m has to be computed and tobe compared to a given Fourier value. To find a nonzero Fourier value needs lessthan m comparisons, see also Remark 3.4. Hence only O ( m log m + m log( N/m )) = O ( m log N ) arithmetical operations are needed to recover x in the case of noisy data. Remark 4.4
For the computation of x it is sufficient to compute x ( L +1) as well as the firstsupport index µ ( J ) , where µ ( j +1) , j = L + 1 , . . . , J − , is iteratively computed from µ ( j ) . The values a j +1 can be obtained directly from x ( L +1) and µ ( j ) . Hence, there isno reason to compute the intermediate vectors x ( j +1) in step 3 of Algorithm 4.3, theyare only given for better illustration. In this section, we discuss the behavior of the algorithm for noisy input data. Wereconstruct randomly generated vectors x from disturbed Fourier data (cid:98) y = (cid:98) x + ε ,where we assume uniform noise ε = ( ε k ) N − k =0 with | ε k | ≤ δ at different noise levels.As a noise measure, we use the SNR valueSNR = 20 · log (cid:107) (cid:98) x (cid:107) (cid:107) ε (cid:107) . et us first illustrate the algorithm for a vector x of length N = 2 = 256 and supportlength m = 6, i.e., we have J = 8 and L = 3. The nonzero entries of x are x = 8, x = − x = − x = 2. We add a noise vector ε of SNR = 20 to (cid:98) x and reconstruct x from (cid:98) y = (cid:98) x + ε . For the noise vector ε in this example, itholds that (cid:107) ε (cid:107) ∞ = 1 .
721 and (cid:107) ε (cid:107) /
256 = 0 . x , and wecompare the reconstruction results of our algorithm to the result of the inverse Fouriertransform applied to (cid:98) y . The reconstruction x (cid:48) computed by our deterministic sparseFFT algorithm has nonzero components x (cid:48) = 7 . .
131 i, x (cid:48) = − . − .
149 i, x (cid:48) = − .
100 + 0 .
171 i, x (cid:48) = − .
955 + 0 .
094 i, x (cid:48) = − .
105 + 0 .
022 i, x (cid:48) =1 . − .
076 i, and an error (cid:107) x − x (cid:48) (cid:107) /
256 = 0 . (cid:107) x − F − (cid:98) y (cid:107) /
256 = 0 . µ ( L +1) in step 2 ofAlgorithm 4.3, i.e., we only use the two vectors (cid:101) z (0) and (cid:101) z ( J − L − = (cid:101) z (3) of length2 L +1 = 16 in this step. Thus, we have taken 32 + 4 = 36 Fourier values to recover x . Let us first illustrate the algorithm for a vector x of length N = 2 = 256 and supportlength m = 6. The nonzero entries of x are x = 8, x = x = x = 2. We add a noise vector " of SNR = 20 to b x and reconstruct x from b y = b x + " .For the noise vector " in this example, it holds that k " k k = 1 .
721 and k " k k /
256 =0 . x as well as we compare the reconstruction resultsfor our algorithm and the inverse Fourier transform on b y . The reconstruction x computed by our deterministic sparse FFT algorithm has nonzero components x =7 . .
131 i, x = . .
149 i, x = . .
171 i, x = . .
094 i, x = .
105 + 0 .
022 i, x = 1 . .
076 i, and an error k x x k /
256 = 0 . k x F b y k /
256 = 0 . µ ( L +1) in step 2 of Algorithm 4.3, i.e., we only use two vectors of length 2 L +1 inthis step. − x − Reconstruction of x by deterministic sparse FFT algorithm − Reconstruction by inverse FFT
Figure 2:
Reconstruction of a vector x of length N = 256:results of the deterministic sparse FFT algorithm and a direct inverse FFT. In Figure 3 we show results for vectors x of length N = 2 and support length m = 50 where the vectors are randomly generated with | Re( x k ) | | Im( x k ) | k = 0 , . . . ,
1. We consider SNR values between 0 and 50 and compute thereconstruction x for 100 vectors x at each noise level. The reconstructed vectors x are evaluated by computing the norm k x x k /N . Additionally, we comparethe results of the reconstruction by our algorithm to the result of an inverse Fouriertransform applied to the noisy vector b y . Up to a noise level of SNR = 10, ouralgorithm makes use of additional vectors in step 2 of Algorithm 4.3 in order toimprove the identification of the first support index µ ( L +1) of x ( L +1) . Here, thealgorithm evaluates at most five additional vectors (SNR = 0), hence a total numberof seven vectors is used to determine µ ( L +1) in this case. (a) Let us first illustrate the algorithm for a vector x of length N = 2 = 256 and supportlength m = 6. The nonzero entries of x are x = 8, x = x = x = 2. We add a noise vector " of SNR = 20 to b x and reconstruct x from b y = b x + " .For the noise vector " in this example, it holds that k " k k = 1 .
721 and k " k k /
256 =0 . x as well as we compare the reconstruction resultsfor our algorithm and the inverse Fourier transform on b y . The reconstruction x computed by our deterministic sparse FFT algorithm has nonzero components x =7 . .
131 i, x = . .
149 i, x = . .
171 i, x = . .
094 i, x = .
105 + 0 .
022 i, x = 1 . .
076 i, and an error k x x k /
256 = 0 . k x F b y k /
256 = 0 . µ ( L +1) in step 2 of Algorithm 4.3, i.e., we only use two vectors of length 2 L +1 inthis step. − x − Reconstruction of x by deterministic sparse FFT algorithm − Reconstruction by inverse FFT
Figure 2:
Reconstruction of a vector x of length N = 256:results of the deterministic sparse FFT algorithm and a direct inverse FFT. In Figure 3 we show results for vectors x of length N = 2 and support length m = 50 where the vectors are randomly generated with | Re( x k ) | | Im( x k ) | k = 0 , . . . ,
1. We consider SNR values between 0 and 50 and compute thereconstruction x for 100 vectors x at each noise level. The reconstructed vectors x are evaluated by computing the norm k x x k /N . Additionally, we comparethe results of the reconstruction by our algorithm to the result of an inverse Fouriertransform applied to the noisy vector b y . Up to a noise level of SNR = 10, ouralgorithm makes use of additional vectors in step 2 of Algorithm 4.3 in order toimprove the identification of the first support index µ ( L +1) of x ( L +1) . Here, thealgorithm evaluates at most five additional vectors (SNR = 0), hence a total numberof seven vectors is used to determine µ ( L +1) in this case. (b) Let us first illustrate the algorithm for a vector x of length N = 2 = 256 and supportlength m = 6. The nonzero entries of x are x = 8, x = x = x = 2. We add a noise vector " of SNR = 20 to b x and reconstruct x from b y = b x + " .For the noise vector " in this example, it holds that k " k k = 1 .
721 and k " k k /
256 =0 . x as well as we compare the reconstruction resultsfor our algorithm and the inverse Fourier transform on b y . The reconstruction x computed by our deterministic sparse FFT algorithm has nonzero components x =7 . .
131 i, x = . .
149 i, x = . .
171 i, x = . .
094 i, x = .
105 + 0 .
022 i, x = 1 . .
076 i, and an error k x x k /
256 = 0 . k x F b y k /
256 = 0 . µ ( L +1) in step 2 of Algorithm 4.3, i.e., we only use two vectors of length 2 L +1 inthis step. − x − Reconstruction of x by deterministic sparse FFT algorithm − Reconstruction by inverse FFT
Figure 2:
Reconstruction of a vector x of length N = 256:results of the deterministic sparse FFT algorithm and a direct inverse FFT. In Figure 3 we show results for vectors x of length N = 2 and support length m = 50 where the vectors are randomly generated with | Re( x k ) | | Im( x k ) | k = 0 , . . . ,
1. We consider SNR values between 0 and 50 and compute thereconstruction x for 100 vectors x at each noise level. The reconstructed vectors x are evaluated by computing the norm k x x k /N . Additionally, we comparethe results of the reconstruction by our algorithm to the result of an inverse Fouriertransform applied to the noisy vector b y . Up to a noise level of SNR = 10, ouralgorithm makes use of additional vectors in step 2 of Algorithm 4.3 in order toimprove the identification of the first support index µ ( L +1) of x ( L +1) . Here, thealgorithm evaluates at most five additional vectors (SNR = 0), hence a total numberof seven vectors is used to determine µ ( L +1) in this case. (c) Figure 2: (a) Original vector x of length N = 256; (b) Reconstruction of x using the sparseFFT Algorithm 4.3; (c) Reconstruction of x using the inverse FFT. In Figure 3 we show the reconstruction error using the sparse FFT Algorithm 4.3for vectors x of length N = 2 and support length m = 50 where the vectors arerandomly generated with | Re( x k ) | ≤ | Im( x k ) | ≤
10 for k in the support interval.We consider noisy Fourier data with SNR values between 0 and 50 and computethe reconstruction x (cid:48) for 100 vectors x at different noise levels. The quality of thereconstructed vectors x (cid:48) is evaluated by computing the norm (cid:107) x − x (cid:48) (cid:107) /N . Figure3 shows the average error norm over all 100 considered vectors. We compare thereconstruction results of our algorithm to the results of an inverse Fourier transformapplied to the noisy vectors (cid:98) y . For noise levels SNR ≤
10, our algorithm made use of dditional vectors (cid:101) z ( κ ) in step 2 of Algorithm 4.3 in order to improve the identificationof the first support index µ ( L +1) of x ( L +1) . Here, the algorithm evaluated at most fiveadditional vectors (for SNR = 0), hence a total number of up to seven vectors (cid:101) z ( κ ) has been used to determine µ ( L +1) in some cases. − − − − SNR E rr o r sparse FFTIFFT Figure 3:
Average reconstruction error (cid:107) x − x (cid:48) (cid:107) /N for different levels of uniform noise,comparing our deterministic sparse FFT algorithm and inverse FFT. Determining the first index of the support of x (i.e., finding the correct µ = µ ( J ) )is one of the crucial points of the algorithm for noisy input data. If µ ( L +1) is notidentified correctly, the correct support interval cannot be found anymore, even if allshifts in step 3 of Algorithm 4.3 are correct.In a third experiment we again take randomly generated vectors x of length N =2 with support length m = 50 or m = 2 , where | Re( x k ) | ≤ | Im( x k ) | ≤
10 for k in the support interval. SNR N = 2 , m = 50 N = 2 , m = 2 correctly (cid:107) ε (cid:107) ∞ (cid:107) ε (cid:107) /N correctly (cid:107) ε (cid:107) ∞ (cid:107) ε (cid:107) /N identified µ identified µ Table 1:
Percentage of correctly identified µ and average norm of noise in 100 randomly chosenvectors for different noise levels, dependent on length N and support length m of x . able 1 shows the number of cases in which the first support index µ has beendetermined correctly by Algorithm 4.3 for 100 randomly chosen vectors for each noiselevel. Additionally, we give an average for the norm of the noise vectors ε at eachnoise level.The results show that the algorithm succeeds for very short support intervals aswell as for long support intervals compared to the full vector length. In cases where µ could not be determined correctly, the error originates from step 2 of the algorithmwhere µ ( L +1) has to be determined. However, the deviation from the correct supportwas small ( ≤
6) in any case.We can conclude that even for high noise level, the support of the reconstructedvector x is correctly found in most of the cases. The support is always correct fornoise levels with SNR ≥
15. Table 1 also shows that the absolute noise ε k at eachcomponent can be considerably large. It is essentially larger for vectors with largersupport since in this case also the signal energy grows accordingly. Acknowledgement
The research in this paper is partially funded by the project PL 170/16-1 of theGerman Research Foundation (DFG). This is gratefully acknowledged.
References [1] A. Akavia, Deterministic sparse Fourier approximation via fooling arithmeticprogressions, in Proc. 23rd COLT, 2010, pp. 381–393.[2] A. Akavia, Deterministic sparse Fourier approximation via approximating arith-metic progressions, IEEE Trans. Inform. Theory (3) (2014), 1733–1741.[3] A. Gilbert, P. Indyk, M.A. Iwen, and L. Schmidt, Recent developments in thesparse Fourier transform, IEEE Signal Processing Magazine (5) (2014), 91–100.[4] J.J. Green, Calculating the maximum modulus of a polynomial using Steˇckin’sLemma, SIAM J. Numer. Anal. (4) (1999), 1022-1029.[5] H. Hassanieh, P. Indyk, D. Katabi, and E. Price, Near-optimal algorithm forsparse Fourier transform, Proc. 44th annual ACM symposium on Theory of Com-puting, 2012, pp. 563–578.[6] H. Hassanieh, P. Indyk, D. Katabi, and E. Price, Simple and practical algorithmfor sparse Fourier transform, Proc. 23th Annual ACM-SIAM Symposium onDiscrete Algorithms (SODA ’12), 2012, pp. 1183–1194.[7] S. Heider, S. Kunis, D. Potts, and M. Veit, A sparse Prony FFT, Proc. 10th In-ternational Conference on Sampling Theory and Applications (SAMPTA), 2013,pp. 572–575.[8] P. Indyk, M. Kapralov, and E. Price, (Nearly) sample-optimal Fourier transform,Proc. 25th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 14),2014, pp. 480–499.[9] M.A. Iwen, Combinatorial sublinear-time Fourier algorithms, Found. Comput.Math. (2010), 303–338.
10] M.A. Iwen, Improved approximation guarantees for sublinear-time Fourier algo-rithms, Appl. Comput. Harmon. Anal. (2013), 57–82.[11] D. Lawlor, Y. Wang, and A. Christlieb, Adaptive sub-linear time Fourier algo-rithms, Advances in Adaptive Data Analysis (1) (2013), 1350003 (25 pages).[12] S. Pawar and K. Ramchandran, Computing a k -sparse n -length discrete Fouriertransform using at most 4 k samples and O ( k log k ) complexity, IEEE Interna-tional Symposium on Information Theory (2013), pp. 464–468.[13] G. Plonka and M. Tasche, Prony methods for recovery of structured functions,GAMM-Mitt. (2) (2014), 239–258.[14] J. Schumacher, High performance sparse fast Fourier transform, Master Thesis,Computer Science. ETH Z¨urich, Switzerland, 2013.(2) (2014), 239–258.[14] J. Schumacher, High performance sparse fast Fourier transform, Master Thesis,Computer Science. ETH Z¨urich, Switzerland, 2013.