The Sparse Fourier Transform. Haitham Hassanieh

of the flat window function and the standard window function are the same up to a constant factor.

      Figure 2.1 An example flat window function for n = 256. This is the sum of 31 adjacent (1/22, 10⁻⁸, 133) Dolph-Chebyshev window functions, giving a (0.11, 0.06, 2 × 10⁻⁹, 133) flat window function (although our proof only guarantees a tolerance δ = 29 × 10⁻⁸, the actual tolerance is better). The top row shows G and in a linear scale, and the bottom row shows the same with a log scale.

      Proof Let f = (∊ − ∊′)/2, and let F be an standard window function with minimal . We can assume ∊, ∊′ > 1/(2n) (because [−∊n, ∊n] = {0} otherwise), so . Let be the sum of 1 + 2(∊′ + f)n adjacent copies of , normalized to . That is, we define

      so by the shift theorem, in the time domain

      Since for |i| ≤ fn, the normalization factor is at least 1. For each i ∈ [−∊′n, ∊′n], the sum on top contains all the terms from the sum on bottom. The other 2∊′n terms in the top sum have magnitude at most δ/((∊′ + ∊)n) = δ/(2(∊′ + f)n), so . For |i| > ∊n, however, . Thus F′ is an (∊, ∊′, δ, w) flat window function, with the correct w. ■

2.2.2 Permutation of Spectra

      Following Gilbert et al. [2005a], we can permute the Fourier spectrum as follows by permuting the time domain:

      Definition 2.3 Suppose σ⁻¹ exists mod n. We define the permutation P_{σ, a, b} by

      We also define πσ, _b(i) = σ(i − b) mod n.

      Claim 2.3 .

       Proof

      Lemma 2.1 If j ≠ 0, n is a power of two, and σ is a uniformly random odd number in [n], then Pr[σj ∊ [−C, C]] ≤ 4C/n.

      Proof If j = m2^l for some odd m, then the distribution of σj as σ varies is uniform over m′2^l for all odd m′. There are thus 2 · round(C/2^l+1) < 4C/2^l+1 possible values in [−C, C] out of n/2^l+1 such elements in the orbit, for a chance of at most 4C/n. ■

      Note that for simplicity we will only analyze our algorithm when n is a power of two. For general n, the analog of Lemma 2.1 would lose an n/φ(n) = O(log log n) factor, where φ is Euler’s totient function. This will correspondingly increase the running time of the algorithm on general n.

      Claim 2.3 allows us to change the set of coefficients binned to a bucket by changing the permutation; Lemma 2.1 bounds the probability of non-zero coefficients falling into the same bucket.

       2.2.3 Subsampled FFT

      Suppose we have a vector x ∈ ℂⁿ and a parameter B dividing n, and would like to compute ŷi = x̂i_(n/B) for i ∈ [B].

      Claim 2.4 ŷ is the B-dimensional Fourier transform of . Therefore, ŷ can be computed in O(|supp(x)| + B log B) time.

       Proof

       2.2.4 2D Aliasing Filter

      The aliasing filter presented in Section 1.1.2. The filter generalizes to two dimensions as follows:

      Given a 2D matrix , and Br, Bc that divide , then for all (i, j) ∈ [Br] × [Bc] set

      Then, compute the 2D DFT ŷ of y. Observe that ŷ is a folded version of x̂:

       3

       Simple and Practical Algorithm

       3.1 Introduction

      In this chapter, we propose a new sublinear Sparse Fourier Transform algorithm over the complex field. The key feature of this algorithm is its simplicity: the algorithm has a simple structure, which leads to efficient runtime with low big-Oh constant. This was the first algorithm to outperform FFT in practice for a practical range of sparsity as we will show later in Chapter 6.

       3.1.1 Results

      We present an algorithm which has the runtime of

      where n is the size of the signal and k is the number of non-zero frequency coefficients. Thus, the algorithm is faster than FFT for k up to O(n/log n). In contrast, earlier algorithms required asymptotically smaller bounds on k. This asymptotic improvement is also reflected in empirical runtimes. For example, we show in Chapter 6 that for n = 2²², this algorithm outperforms FFT for k up to about 2,200, which is an order of magnitude higher than what was achieved