The Sparse Fourier Transform. Haitham Hassanieh
it does not collide.2 We subtract the green frequency from the colliding bucket in the second bucketization to obtain the blue frequency. We then go back to the first bucketization and subtract the blue frequency from the bucket where it collides to obtain the red frequency.
Figure 1.2 Resolving collisions with co-prime aliasing. Using two co-prime aliasing filters, we ensure the frequencies that collide in one filter will not collide in the second. For example, frequencies 5 and 9 collide in the first filter. But frequency 5 does not collide in the second which allows us to estimate it and subtract it.
Iterating between the different bucketizations by estimating frequencies from buckets where they do not collide and subtracting them from buckets where they do collide, ensures that each non-zero frequency will be isolated in its own bucket during some iteration of the algorithm. This allows us to estimate each non-zero frequency correctly. Thus, at the end of the of the collision resolution step, we have recovered all non-zero frequencies and hence have successfully computed the Fourier transform of the signal.
1.1.3 Algorithmic Techniques
The previous section established a general framework for computing the Sparse Fourier Transform and gave one example of a technique that can be used in each step of this framework. In this section, we describe a more comprehensive list of techniques that are used by different Sparse Fourier Transform algorithms.
Figure 1.3 Filters used for frequency bucketization shown in the time (upper row) and the frequency (lower row) domain.
1.1.3.1 Frequency Bucketization Techniques
As described earlier, bucketization is done using filters. The choice of the filter can severely affect the running time of a Sparse Fourier Transform algorithm. Ideally, we would like a filter that uses a small number of input time samples to hash the frequency coefficients into buckets. For example, the rectangular or boxcar filter shown in Figure 1.3(a), uses only B time samples to hash the frequency coefficients into B buckets which is ideal in time domain. However, in the frequency domain, it is equal to the sinc function,3 which decays polynomially as shown in Figure 1.3(a). This polynomial decay means that the frequency coefficients “leak” between buckets, i.e., the value of the bucket is no longer the sum of n/B coefficients that hash to the bucket. It is a weighted sum of all the n frequency coefficients. Hence, a nonzero frequency coefficient can never be isolated in a bucket and estimated correctly. One the other hand, a rectangular filter is ideal in the frequency domain since it has no leakage, as shown in Figure 1.3(b). However, it is sinc function in the time domain and hence requires using all n input samples which take at least Ω(n) time to process.
In this book, we identify several efficient filters that use a small number of samples in time domain and have minimal or no leakage in the frequency domain and as a result can be used to perform fast bucketization.
Flat Window Filter. This filter looks very similar to a rectangle or box in the frequency domain while still using a small number of time samples. An example of such filter is a Gaussian function multiplied by a sinc function in the time domain which is shown in Figure 1.3(c). Since the Gaussian function decays exponentially fast both in time and frequency, the leakage between buckets in this filter is negligible and can be ignored. Similarly, the filter is concentrated in time and hence uses only a small number of time samples. The resulting hash function of such filter can be written as h(f) = ⌈f/(n/B)⌉. Gaussian is only one example of such functions. One can potentially use a Dolph-Chebyshev or a Kaiser-Bessel function as we describe in more detail in Chapter 2.
Aliasing Filter. We presented this filter in the previous section. It is simply a spike-train of period p in time domain since it subsamples the time domain signal by a factor of p, as shown in Figure 1.3(d). It is also a spike-train in the frequency domain since it sums up frequency coefficients that are equally spaced by n/p. The aliasing filter is ideal both in time and in frequency since it only uses B time samples and has zero leakage between buckets. Unfortunately, we will later show that the aliasing filter does not lend itself to powerful randomization techniques and, as a result, can only be shown to work for average case input signals as opposed to worst case.
Fourier Projection Filter. This filter is a generalization of the aliasing filter to higher dimensions. It is a direct result of the Fourier Slice Projection Theorem which states that taking the Fourier transform of samples along a slice (e.g., a 1D line in 2D time signal or a 2D plane in a 3D signal) results in the orthogonal projection of the frequency spectrum onto the corresponding slice in the Fourier domain. For example, if we sample a row in a 2D time domain signal and then take its 1D Fourier transform, each output point of the Fourier transform will be the sum of the coefficients along one column, as shown in Figure 1.4(a). Alternatively, if we sample a column, each output point will be the sum of the coefficients along one row, as shown in Figure 1.4(b). This also holds for any discrete line, as shown in Figure 1.4(c,d). Thus, if we sample a line with slope equal to 1 or 2, we get a projection in the frequency domain along lines with slopes equal to −1or −1/2. Note that discrete lines wrap around and hence can generate non-uniform sampling, as shown in Figure 1.4(d).
Figure 1.4 Fourier projection filters. The top row of figures shows the sampled lines of the time signal and the bottom row of figures shows how the spectrum is projected. Frequencies of the same color are projected onto the same point: (a) row projection, (b) column projection, (c) diagonal projection, and (d) line with slope = 2.
1.1.3.2 Frequency Estimation Techniques
Recall that the goal of the frequency estimation step is to compute the positions f and values x̂(f) of the non-zero frequency coefficients that have been hashed to non-empty buckets. In what follows, we will present the techniques used by the Sparse Fourier Transform algorithms to perform frequency estimation.
Time-Shift/Phase-Rotation Approach. In this approach, which we have previously described in Section 1.1.2, we repeat the bucketization after shifting the time signal x by a circular shift of τ samples. For buckets with a single isolated non-zero frequency coefficient, this results in a phase rotation Δϕ of the complex value of the coefficient. Δϕ = 2πfτ/n which we can use to compute the frequency position f. Since the frequency is isolated, its value can be immediately computed from the value of the bucket as described in the previous section. This is an extremely efficient approach since it requires constant O(1) time to estimate each non-zero frequency. For noisy signals, we repeat this process for multiple different time shifts to average the noise. The details of how we average the noise to ensure robust recovery can be found in Chapter 4.
Voting