How do you sub-sample an FFT?

Question

I'm currently trying to understand this paper on a Sparse Fourier Transform.

On page four, there is a section on a sub-sampled FFT. The purpose of this section is to show that you can compute a sub-sampled FFT using less computations than the original FFT. So if I have a signal of size N, I could compute a subsampled FFT of size B which would require O(N + B log B) operations.

However, I'm confused as to how you would obtain this subsampled FFT. The idea is that we want the original FFT spectrum to be sampled at B locations. Obviously we could miss coefficients by only sampling the spectrum at B locations, however, it seems like the bigger problem is that subsampling the signal would result in aliasing frequencies.

Essentially I'm not sold on the fact that the described algorithm would result in a sub sampled FFT of the original spectrum. Instead, I think the spectrum could be aliased resulting in incorrect coefficients.

Answer 1

it seems like the bigger problem is that subsampling the signal would result in aliasing frequencies

No, subsampling in frequency domain corresponds to aliasing in time domain. So the idea here is to purposefully alias in time domain so that you get a sub-sampled FFT. That is exactly why 'Claim 3.7' of paper mentions $y_i=\sum_{j=0}^{n/B-1}x_{i+Bj}$. These are $n/B$ copies of $x[n]$ shifted by $B$ and overlapped. Each copy is positioned at $0,B,2B,\ldots,(\frac{n}{B}-1)B$. This aliased signal can be used to compute FFT ($B\log_2 B$ operations) to arrive at subsampled FFT. Total number of operations will the operations need to create the aliasing (summing of non-zero values of $x$ which equals $O(\text{supp}(x))$ in addition to the operation needed to compute FFT of aliased $x$ ($O(B\log_2 B)$).

For example, if a time domain signal is $1024$ samples long, ideally you need $1024$ point FFT. But if you sub sample by taking every 8th sample, $B=128$, $n/B=1024/128=8$, the time domain signal at $\tilde{x}[0]$ will be having addition from $x[0]$,$x[128]$,$x[256]$,$\ldots$,$x[896]$. Like this for $\tilde{x}[n]=\sum_0^{7}x[128k+n]$. You can now use this $\tilde{x}$ to compute $B$-point FFT to arrive at the Subsampled FFT.

Answer 2

It is important to understand that the "sub-sampling" mentioned is occurring in the frequency domain, not the time-domain-- so any aliasing that the OP is thinking of would all be time-domain aliasing, not aliasing in the frequency domain. To "sub-sample" the frequency spectrum, the time domain sampling rate is unchanged but the duration in time (number of samples in the DFT) is reduced. This results in each DFT bin being wider in frequency, which is a reduced frequency resolution.

Here is the general relationship in case this provides further insight:

Consider a sequence $x[n]$ of length $N$, whose DFT would be $X[k]$, and a shorter length $B$ which is equally divisible into $N$ ($N/B$ is an integer $D$).

To create the $B$-length DFT where the result is exactly every $D^\rm{th}$ sample of $X[k]$ you would need to do the following:

• Split $x[n]$ into $D$ $B$-length blocks.

• Take the $B$-length DFT of each of the $D$ blocks.

• Sum the complex results in each bin.

This would result in an exact match of the original and down-sampled DFTs.

This does result in a computational savings of $\log_2(D)$: Comparing the number of operations in the DFT of the $N$-length sequence ($N\log_2(N)$ vs the DFT of $D$ $B$-length sequences $D(B\log_2(B))$, since $D = N/B$:

$$D\left(B\log_2(B)\right) = \frac NB\left(M\log_2(B)\right) = N\log_2(B)$$

There would be $N\log_2(B)$ operations (plus the $D$ summations) to compute a down-sampled DFT spectrum in comparison to the $N\log_2(N)$ operations for the original $N$-length sequence.

Doing the above is not the best solution!

The operations above describe true down-sampling (select every $D$th sample, throw away the rest.) However, by simply taking the DFT of one $B$ length sequence, which would further reduce the number of operations required by $B$, would result in a truly decimated result. With down-sampling alone as the above process describes and as described in the paper, we can lose critical information: consider a case where every $D^\rm{th}$ sample is zero but there is spectral content in every other bin. The process above would result in all zeros! This is the result of the time-aliasing which is not a desired property.

Decimation is the combination of filtering and down-sampling. The result of a DFT over a single $B$-length sequence would not be the exact down-sampled value of the DFT as described above but would be the result of the average value over the adjacent $D$ bins. This would be equivalent to the combination of a moving average in frequency over $D$ samples of $X[k]$ followed by the down-sampling of every $D^\rm{th}$ sample--which would be then the decimation of the frequency spectrum. Further, the implied moving average in frequency eliminates the time-domain aliasing (to the same extent that doing a moving average in time before resampling would eliminate aliasing in the frequency domain).

Answer 3

Another way to look at this is as a windowed FFT with a variable size window different from the FFT length, and where the window can be both narrower or wider than the FFT.

If the window aperture is shorter than the FFT, then the FFT has to be zero-padded, and the frequency response of each result bin of the FFT is proportionally wider or fatter than if the FFT was fed enough actual data to span the full width of the FFT.

If the data window is longer than the FFT, then the data has to be circularly wrapped (and summed) to fit it all into the FFT, and the response of each FFT result bin is proportionally narrower, thus becoming a narrow "sample" of the wider "normal" full data length FFT Sinc response.

So you can vary the frequency response of each FFT bin, both narrower and wider than the default Sinc (actually Dirichlet) response by varying the window, not only in shape (Von Hann, et.al.), but in width.