I’m currently a little fed up with number theory, so its time to change topics completely. Specially since the post on basic integer factorization completes what I believe is a sufficient toolkit to tackle a very cool subject: the fast Fourier transform (FFT).

I have some mixed feelings about how does Fourier analysis qualify for the “uncomplicated complexity” rule I imposed on myself when starting this blog. There certainly is a lot of very elaborate math behind it, but I think that if you skip all the mathematical subtleties and dive head first into a “the Fourier transform gives you the frequencies present in the signal” type of intuitive description, the journey from the more naive implementations of the discrete Fourier transform (DFT) to the many flavors of the FFT is a very enjoyable one, which requires only a limited amount of complication, but gentle doses of complexity. Which is exactly what I’m after…

So I will take the definition of the DFT for granted, and concentrate on the algorithms to compute it efficiently. There will have to be a few leaps of faith, but I will try to mark them clearly, so you can whisper “amen” and keep going without a second thought. Although things may unfold differently, I think there is material for about half a dozen posts. Towards the end I may devote some time to optimizing the code, with things such as avoiding recursion, or doing the calculations in place, but be warned that most of my efforts will be devoted to optimizing the math behind the code, not the code itself.

**The Discrete Fourier Transform**

Without further explanation, we will begin by writing down the analytical expression of the DFT,

and of its corresponding inverse transform,

With the built-in support for complex arithmetic, there really isn’t much mistery in turning these two formulas into python functions, or as I have chosen, one with an `inverse`

switch:

from __future__ import division import math import time def dft(x, inverse = False, verbose = False) : t = time.clock() N = len(x) inv = -1 if not inverse else 1 X =[0] * N for k in xrange(N) : for n in xrange(N) : X[k] += x[n] * math.e**(inv * 2j * math.pi * k * n / N) if inverse : X[k] /= N t = time.clock() - t if verbose : print "Computed","an inverse" if inverse else "a","DFT of size",N, print "in",t,"sec." return X

Of course, having two nested loops of the size of the input means that computation time will grow with N^{2}, which is extremely inefficient.

**Twiddle Factors**

Apart from the FFT algorithms we will start with in the next post, there are some pretty straightforward optimizations at hand. While some texts reserve the name for similar factors appearing in FFT algorithms, the complex exponentials coming up in the DFT formula are often known as twiddle factors. While in principle there are N^{2} of them to calculate, we can take advantage of the fact that

,

and owing to the periodic nature of trigonometric functions, it isn’t hard to see that n·k can be replaced by n·k mod N without altering the result, which means only N different ones do exist.

A python implementation with pre-calculated twiddle factors could look something like this…

from __future__ import division import math import time def dft(x, inverse = False, verbose = False) : t = time.clock() N = len(x) inv = -1 if not inverse else 1 twiddles = [math.e**(inv*2j*math.pi*k/N) for k in xrange(N)] X =[0] * N for k in xrange(N) : for n in xrange(N) : X[k] += x[n] * twiddles[n * k % N] if inverse : X[k] /= N t = time.clock() - t if verbose : print "Computed","an inverse" if inverse else "a","DFT of size",N, print "in",t,"sec." return X

On my machine this code is about 5 times faster than the original one, although it still scales with time as N^{2}, so we are only postponing the inevitable. Do notice as well that the increase in speed comes at the cost of an additional array of complex of size N memory requirement.

**Particular Cases**

Without being touched with insightful inspiration, the only other route to optimization is by working out the math on paper for particular cases, and code explicitly the resulting formulas to avoid doing costly exponentiations.

As examples, if N = 2, it is easy to show, that

and if N = 4, one can similarly arrive at

Adding these particular cases to the python code could be done as follows…

def dft(x, inverse = False, verbose = False) : t = time.clock() N = len(x) X =[0] * N inv = -1 if not inverse else 1 if N == 1 : X[0] = x[0] elif N == 2 : X[0] = x[0] + x[1] X[1] = x[0] - x[1] elif N == 4 : a = x[0] + x[2] b = x[0] - x[2] c = x[1] + x[3] d = x[1] - x[3] X[0] = a + c X[1] = b + inv * 1j * d X[2] = a - c X[3] = b - inv * 1j * d if inverse : X = [k/N for k in X] else : twiddles = [math.e**(inv*2j*math.pi*k/N) for k in xrange(N//2+1)] for k in xrange(N//2+1,N) : twiddles += [twiddles[N-k].conjugate()] for k in xrange(N) : for n in xrange(N) : X[k] += x[n] * twiddles[n * k % N] if inverse : X[k] /= N t = time.clock() - t if verbose : print "Computed","an inverse" if inverse else "a","DFT of size",N, print "in",t,"sec." return X

Overall, and for this particular cases only, an additional factor of 3 increase in speed is achieved in this way. But this solution is not scalable to arbitrary input sizes, and so wiser approaches are needed to tackle even moderately large inputs of few tens of thousands samples…

Hello, Nice tutorial!

I just wanted to ask you one thing though. Is it possible to get a value of X(k) where k>N? I mean if the input signal has 5 values, would be possible to obtain such a value as X[6]. The mathematics seems to allow, but what would that mean?

With zero based arrays, if k > N, then X(k) = X(k % N).

You may want to read this, as it elaborates on the effects of windowing and sampling of a real world signal:

http://en.wikipedia.org/wiki/A_derivation_of_the_discrete_Fourier_transform

Hi,

thanks for the great tutorial!

I’d like to use and adapt this code in my thesis project so it would be good to know under which license it is released.

Thanks in advance.

Hello,

Very good tutorial.

I wanted to know how to do this for a pdb (protein data bank) file.

Thanks and regards,