What are some of the differences between DFT and FFT that make FFT so fast?

Question

I'm trying to understand FFTs, here's what I have so far:

In order to find the magnitude of frequencies in a waveform, one must probe for them by multiplying the wave by the frequency they are searching for, in two different phases (sin and cos) and averaging each. The phase is found by its relation to the two, and the code for that is something like this:

//simple pseudocode
var wave = [...];                //an array of floats representing amplitude of wave
var numSamples = wave.length;
var spectrum = [1,2,3,4,5,6...]  //all frequencies being tested for.  

function getMagnitudesOfSpectrum() {
   var magnitudesOut = [];
   var phasesOut = [];

   for(freq in spectrum) {
       var magnitudeSin = 0;
       var magnitudeCos = 0;

       for(sample in numSamples) {
          magnitudeSin += amplitudeSinAt(sample, freq) * wave[sample];
          magnitudeCos += amplitudeCosAt(sample, freq) * wave[sample];
       }

       magnitudesOut[freq] = (magnitudeSin + magnitudeCos)/numSamples;
       phasesOut[freq] = //based off magnitudeSin and magnitudeCos
   }

   return magnitudesOut and phasesOut;
}

In order to do this for very many frequencies very quickly, FFTs use many tricks.

What are some of the tricks used to make FFTs so much faster than DFT?

P.S. I have tried looking at completed FFT algorithms on the web, but all the tricks tend to be condensed into one beautiful piece of code without much explanation. What I need first, before I can understand the entire thing, is some introduction to each of these efficient changes as concepts.

Thank you.

Answer 1

The naive implementation of an $N$-point DFT is basically a multiplication by a $N \times N$ matrix. This results in a complexity of $\mathcal{O}(N^2)$.

One of the most common Fast Fourier Transform (FFT) algorithm is the radix-2 Cooley-Tukey Decimation-in-Time FFT algorithm. This is a basic divide and conquer approach.

First define the "twiddle factor" as: $$W_N \triangleq e^{-j\frac{2\pi}{N}}$$ where $j \triangleq \sqrt{-1}$ is the imaginary unit, then the DFT $X[k]$ of $x[n]$ is given by $$X[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{kn}$$ If $N$ is even (and $\tfrac{N}{2}$ is an integer), the sum can then be divided in two sums as follows $$X[k] = \sum_{n=0}^{N/2-1} x[2n]W_N^{2kn} + \sum_{n=0}^{N/2-1} x[2n+1]W_N^{k(2n+1)}$$ where the first summation deals with the even samples of $x[n]$ and the second with the odd samples of $x[n]$. Defining $x_e[n] \triangleq x[2n]$ and $x_o[n] \triangleq x[2n+1]$ and using the fact that

1. $W_N^{k(2n+1)} = W_N^{2kn}W_N^k$ and
2. $W_N^{2kn} = W_{N/2}^{kn}$

this can be re-written as $$ \begin{align} X[k] &= \sum_{n=0}^{N/2-1} x_e[n] W_{N/2}^{kn} + W_N^k\sum_{n=0}^{N/2-1} x_o[n]W_{N/2}^{kn} \\ & = X_e[k] + W_N^k X_o[k] \end{align} $$ where $X_e[k]$ and $X_o[k]$ are the $\tfrac{N}{2}$-point DFT transforms of the even and odd samples of $x[n]$ respectively. So we just transformed a single $N$-point DFT into two smaller $\tfrac{N}{2}$-point DFTs. This reduces computational cost because $$ 2 \left( \frac{N}{2} \right)^2 + N < N^2 $$ when $N > 2$.

We can then re-iterate the same process on those two smaller DFTs. This divide-and-conquer approach allows to reach complexity of $\mathcal{O}(N\log N)$, which is way better than the $\mathcal{O}(N^2)$ we had with the naive DFT implementation (as is greatly illustrated by leftaroundabout's answer).

Answer 2

http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629

DFT, size 16

FFT, size 16

The difference in complexity is pretty evident from that, isn't it?

---

Here's how I understand FFT.

First off, I would always think about Fourier transforms foremostly as transforms of continuous functions, i.e. a bijective mapping $\operatorname{FT} : \mathcal{L}^2(\mathbb{R}) \to \mathcal{L}^2(\mathbb{R})$. In that light it's clear that it can't be really be necessary to go to the “deepest level” and loop over individual elements, because the “individual elements” are single points on the real line, of which there are over-countably infinite.

So how comes this transform is still well-defined? Well, it's crucial that it operates not on the general function space $\mathbb{R}\to\mathbb{C}$ but only on the space of (Lebesgue-, square-) integrable functions. Now, this integrability is not a very strong property (much weaker than differentiability etc.), but it does demand that the function becomes “locally discribable with countable information”. Such a discription is given by the coefficients of a short-time Fourier Transform.^† The simplest case is that your function is continuous and you divide it in so small regions that it's basically constant in each of them. Then each of the STFTs has most strongly a zeroth term. If you ignore the (anyways decaying) other coefficients then each domain is just a single data point. Of all these short-time–LF-limit coefficients, you could take a discrete Fourier transform. In fact, that's exactly what you do when performing any FT on measured real-world data!

The measured data needn't necessarily correspond to a fundamental physical quantity, though. For instance, when you measure some light intensity, you're really just measuring the amplitude of an electromagnetic wave whose frequency is itself too high to be sampled with an ADC. But clearly you can also compute the DFT of a sampled light-intensity signal, and cheaply so, despite the insane frequency of the light-wave.

This could be understood as the most important reason FFT is cheap:

Don't bother trying to see the individual oscillation cycles from the highest level. Instead, transform only somewhat high-level information that's already been preprocessed locally.

That's not all there is to it, though. The great thing about FFT is that is still gives you all the information a complete DFT would give. I.e. all the information you'd also get when sampling the exact electromagnetic wave of a light-beam. Can this be accomplished by transforming a photodiode signal? – can you measure the exact light-frequency from that?

Well, the answer is no you can't. That is, unless you apply extra tricks.
First of all, you need to at least roughly measure the frequency in the short time blocks. Well, that's possible with a spectrograph. But it's only possible up to a precision of $\Delta \nu = 1/{\Delta t}$, a typical uncertainty relation^‡.

By having overall a longer time span, we should also be able to narrow down the frequency uncertainty. And this is indeed possible, if you measure locally not just the rough-frequency but also the phase of the wave. You know that a 1000 Hz signal will have exactly the same phase if you look at it one second later. Whereas a 1000.5 Hz signal, while being indistinguishable on the short-scale, will have inverted phase one second later.

Luckily, that phase information can very well be stored in a single complex number. And that's how FFT works! It starts out with lots of small, local transformations. These are cheap – for one thing obviously because they only use a small amount of data, but secondly because they know that, due to the short time span, they can't resolve the frequency very precisely anyway – so it's still affordable even though you do a whole lot of such transformations.

These do, however, record also the phase, and from that you can then make the frequency resolution more exact on the top-level. The required transformation is again cheap, because it doesn't itself bother with any high-frequency oscillations but only with the pre-processed low-frequency data.

---

^†_{Yup, my argumentation is a bit circular at this point. Let's just call it recursive and we're fine...}

^‡_{This relation is not quantum mechanical, but the Heisenberg uncertainty has actually the same fundamental reason.}

Answer 3

Here is a picture to add to Robert's good answer demonstrating the "re-use" of operations, in this case for an 8 point DFT. The "Twiddle Factors" are represented in the diagram using the notation $W_N^{nk}$ which is equal to $e^{j2\pi \frac{nk}{N}}$

Note the path shown and the equation underneath shows the result for the frequency bin X(1), as given by Robert's equation copied below:

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{j 2 \pi \frac{nk}{N}} $$

Dashed lines are no different than solid lines just to make clear where the summation joins are.

Answer 4

essentially, in computing the naive DFT directly from the summation:

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} $$

there are $N$ table lookups for the twiddle factor $ e^{-j 2 \pi \frac{nk}{N}} $, $N$ complex multiplications, and $N-1$ additions. and that's just for one value of $X[k]$ and one instance of $k$. then the naive DFT throws away all of that intermediate data away and goes through all of it again for $X[k+1]$.

1. so the FFT holds on to some intermediate data.
2. the FFT will also make use of factoring the twiddle factor a bit so that the same factor can be used for an intermediate combination of data.

Answer 5

I am a visual person. I prefer to imagine the FFT as a matrix trick rather than as a summation trick.

To explain at a high level:

A naive DFT computes each output sample independently and uses every input sample in each computation (classic N² algorithm).

A common FFT uses symmetries and patterns in the DFT definition to do the computation in "layers" (log N layers), each layer with constant-time requirement per sample creating an N log N algorithm.

More specifics:

One way to visualize these symmetries is to look at the DFT as a 1×N matrix input multiplied by an NxN matrix of all your complex exponentials. Let's start with the "radix 2" case. We're going to split out the even and odd rows of the matrix (corresponding to the even and odd input samples) and consider them as two separate matrix multiplications which add together to get the same final result.

Now look at these matrices: in the first one the left half is identical to the right half. In the other, the right half is the left half x −1. This means we only really have to use the left half of these matrices for multiplication and create the right half cheaply by multiplying by 1 or −1. Next, observe that the second matrix differs from the first matrix by factors that are the same in each column, so we can factor that out and multiply it into the input so now both even and odd samples use the same matrix, but require a multiplier first. And the final step is observing that this resulting N/2 × N/2 matrix is identical to an N/2 DFT matrix and we can do this again and again until we reach a 1×1 matrix where the DFT is an identity function.

To generalize beyond radix 2, you can look at splitting every third row and looking at three chunks of columns, or every 4th etc.

In the event of prime sized inputs, there exists a method to properly zero-pad, FFT, and truncate, but that is beyond the scope of this answer.

See: http://whoiskylefinn.com/MatrixFFT.html

Answer 6

The DFT does a brute force N^2 matrix multiply.

FFTs does clever tricks, exploiting properties of the matrix (degeneralizing the matrix multiply) in order to reduce computational cost.

Let us first look at a small DFT:

W=fft(eye(4));

x = rand(4,1)+1j*rand(4,1);

X_ref = fft(x);

X = W*x;

assert(max(abs(X-X_ref)) < 1e-7)

Great so we are able to substitute MATLABs call to the FFTW library by a small 4x4 ( complex) matrix multiplication by filling a matrix from the FFT function. So what does this matrix look like?

N=4,

Wn=exp(-1j*2*pi/N),

f=((0:N-1)'*(0:N-1))

f =

 0     0     0     0
 0     1     2     3
 0     2     4     6
 0     3     6     9

W=Wn.^f

W =

1 1 1 1

1 -i -1 i

1 -1 1 -1

1 i -1 -i

Each element is either +1, -1, +1j or -1j. Obviously, this means that we can avoid full complex multiplications. Further, the first column is identical, meaning that we are multiplying the first element of x over and over by the same factor.

It turns out that Kronecker tensor products, "twiddle factors" and a permutation matrix where the index is changed according to the binary represantation flipped is both compact and gives an alternate perspective on how FFTs are computed as a set of sparse matrix operations.

The lines below is a simple Decimation in Frequency (DIF) radix 2 forward FFT. While the steps may seem cumbersome, it is convenient to reuse for forward/inverse FFT, radix4/split-radix or decimation-in-time, while being a fair representation of how in-place FFTs tends to be implemented in the real world, I believe.

N = 4;

x = randn(N, 1) +1j*randn(N, 1);

T1 = exp(-1j*2*pi*([zeros(1, N/2), 0:(N/2-1)]).'/N),

M0 =kron(eye(2), fft(eye(2))),

M1 = kron(fft(eye(2)), eye(2)),

X=bitrevorder(x.'*M1*diag(T1)*M0),

X_ref=fft(x)

assert(max(abs(X(:)-X_ref(:)))<1e-6)

C F Van Loan has a great book on this subject.

Answer 7

If you want to drink from the Firehose of Wisdom, I suggest :

"Fast Transforms - Algorithms, Analyses, Applications" by Douglas F. Elliott, K. Ramamohan Rao

It covers FFT, Hartley, Winograd and applications.

One strong point is that is shows how the FFT is a set of sparse matrix factorizations with bit reversal ordering.