The Sparse Fourier Transform. Haitham Hassanieh
The best k-sparse approximation of x̂ can be obtained by setting all but the largest k coefficients of x̂ to 0. The goal is to compute an approximation x̂′ in which the error in approximating x̂ is bounded by the error on the best k-sparse approximation. Formally, x̂ has to satisfy the following ℓ2/ℓ2 guarantee:
where C is some approximation factor and the minimization is over exactly k-sparse signals.
In the remainder of this section, we will describe the algorithmic framework and techniques in terms of exactly sparse signals. However, the full details and extensions to the general case of approximately sparse signals can be found in Chapters 3, 4, and 5.
1.1.2 Algorithmic Framework
The Sparse Fourier Transform has three main components: Frequency Bucketization, Frequency Estimation, and Collision Resolution.
1.1.2.1 Frequency Bucketization
The Sparse Fourier Transform starts by hashing the frequency coefficients of x̂ into buckets such that the value of the bucket is the sum of the values of the frequency coefficients that hash into the bucket. Since x̂ is sparse, many buckets will be empty and can be simply discarded. The algorithm then focuses on the non-empty buckets and computes the positions and values of the large frequency coefficients in those buckets in what we call the frequency estimation step.
The process of frequency bucketization is achieved through the use of filters. A filter suppresses and zeroes out frequency coefficients that hash outside the bucket while passing through frequency coefficients that hash into the bucket. The simplest example of this is the aliasing filter. Recall the following basic property of the Fourier transform: subsampling in the time domain causes aliasing in the frequency domain. Formally, let b be a subsampled version of x, i.e., b(i) = x(i · p) where p is a subsampling factor that divides n. Then, b̂, the Fourier transform of b is an aliased version of x̂, i.e.,
Thus, an aliasing filter is a form of bucketization in which frequencies equally spaced by an interval B = n/p hash to the same bucket and there are B such buckets, as shown in Figure 1.1. The hashing function resulting from this bucketization can be written as: h(f) = f mod n/p. Further, the value in each bucket is the sum of the values of only the frequency coefficients that hash to the bucket as can be seen from Equation 1.3.
For the above aliasing filter, the buckets can be computed efficiently using a B-point FFT which takes O(B log B) time. We set B = O(k) and hence bucketization takes only O(k log k) time and uses only O(B) = O(k) of the input samples of x. In Section 1.1.3, we describe additional types of filters that are used by our Sparse Fourier Transform algorithms to perform frequency bucketization.
Figure 1.1 Bucketization using aliasing filter. Subsampling a signal by 3× in the time domain results in the spectrum aliasing. Specifically, the 12 frequency will alias into 4 buckets. Frequencies that are equally spaced by 4 (shown with the same color) end up in the same bucket.
1.1.2.2 Frequency Estimation
In this step, the Sparse Fourier Transform estimates the positions and values of the non-zero frequency coefficients which created the energy in each of the non-empty buckets. Since x̂ is sparse, many of the non-empty buckets will likely have a single non-zero frequency coefficient hashing into them, and only a small number will have a collision of multiple non-zero coefficients. We first focus on buckets with a single non-zero frequency coefficients and estimate the value and the position of this non-zero frequency, i.e., x̂(f) and the corresponding f.
In the absence of a collision, the value of the non-zero frequency coefficient is the value of the bucket it hashes to since all other frequencies that hash into the bucket have zero values. Hence, we can easily find the value of the non-zero frequency coefficient in a bucket. However, we still do not know its frequency position f, since frequency bucketization mapped multiple frequencies to the same bucket. The simplest way to compute f is to leverage the phase-rotation property of the Fourier transform, which states that a shift in time domain translates into phase rotation in the frequency domain [Lyons 1996]. Specifically, we perform the process of bucketization again, after a circular shift of x by τ samples. Since a shift in time translates into a phase rotation in the frequency domain, the value of the bucket changes from b̂(i) = x̂(f) to b̂(τ)(i) = x̂(f) · ejΔϕ where the phase rotation is
Hence, using the change in the phase of the bucket, we can estimate the position of the non-zero frequency coefficient in the bucket. Note that the phase wraps around every 2π and so the shift τ should be 1 to avoid the phase wrapping for large values of f.1 Since there are k non-zero frequency coefficients, this frequency estimation can be done efficiently using at most O(k) computations. In Section 1.1.3, we describe additional techniques that are used by our Sparse Fourier Transform algorithms to estimate the values and positions of non-zero frequency coefficients.
1.1.2.3 Collision Resolution
Non-zero frequency coefficients that are isolated in their own bucket can be properly estimated, as described above. However, when non-zero frequencies collide in the same bucket, we are unable to estimate them correctly. Hence, to recover the full frequency spectrum, we need to resolve the collisions.
To resolve collision, we need to repeat the frequency bucketization in a manner that ensures that the same non-zero frequencies do not collide with each other every time. The manner in which we achieve this depends on the type of filter used for bucketization. For example, with the aliasing filters described above, we can bucketize the spectrum multiple times using aliasing filters with co-prime sampling rates. This changes the hashing function from h(f) = f mod n/p to h′(f) = f mod n/p′ where p and p′ are co-prime. Co-prime aliasing filters guarantee that any two frequencies that collide in one bucketization will not collide in the other bucketization. To better understand this point, consider the example in Figure 1.2. The first time we bucketize, we use an aliasing filter that subsamples the time signal by a factor of 3. In this case, the two frequencies labeled in red and blue collide in a bucket whereas the frequency labeled in green does not collide, as shown in the figure. The second time we bucketize, we use an aliasing filter that subsamples by 4. This time the blue and green frequencies collide whereas the red frequency does not collide. Now we can resolve collisions by iterating between the two bucketizations. For example, we can estimate the green frequency from the first