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.

  • 7
    "DFT" doesn't refer to an algorithm: it refers to a mathematical operation. "FFT" refers to a class of methods for computing that operation. – Hurkyl Jun 9 '17 at 4:16
  • 1
    Just wanted to point out that the use of sudo in your code example could be confusing, as that is a well known command in the computer world. You probably meant psuedocode. – rwfeather Jun 9 '17 at 16:24
  • 1
    @nwfeather He probably meant 'pseudocode'. – user207421 Jun 9 '17 at 16:45

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).

  • would you be willing to list out what each of the variables stands for? I'm rather new to this, so W, j, X(), N and k don't yet have definitions for me. – Seph Reed Jun 8 '17 at 20:44
  • $W$ is already defined in my answer. I tried to better define some other notations. $k$ denotes index in the frequency domain and $n$ index in the time domain. – anpar Jun 8 '17 at 20:53

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

DFT, size 16

Diagram of the operations in a size-16 naïve DFT

FFT, size 16

Diagram of the operations in a size-16 radix-2 FFT

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.

  • 2
    nice pictoral depiction of the issue. :-) – robert bristow-johnson Jun 9 '17 at 18:12
  • 2
    Don't you love diagrams that are repeated everywhere and never actually explained anywhere :) – Mehrdad Jun 9 '17 at 23:03
  • 1
    I understood the picture after having just read anpar’s answer. – JDługosz Jun 11 '17 at 16:48

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.

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

FFT implementation

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.

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

  • prime FFT, various FFT. Using zero-pad is not the only option. Sorry, I just find zero-padding overused. One small question, I do not understand what you mean by "each layer with constant-time requirement per sample", if you could explain, it would be awesome. – Evil Jun 9 '17 at 22:30
  • 1
    Sorry I didn't mean to say zero padding was THE way, just wanted to point to further reading. And "layer" meaning a recursion, or a translation from a N DFT to 2 N/2 DFTs, with constant time per sample meaning this step is O(N). – kylefinn Jun 9 '17 at 23:20
  • So far, of all the descriptions, this one seems the closest to making a complex issue simple. The big thing it's missing, though, is an example of these matrix's. Would you happen to have one? – Seph Reed Jun 10 '17 at 18:30
  • Uploaded this, should help: whoiskylefinn.com/MatrixFFT.html – kylefinn Jun 13 '17 at 1:51

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.

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.

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.