The Discrete Fourier Transform

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,

\displaystyle X_k = \sum_{n=0}^{N-1}{x_n e^{-2\pi i k n / N}},

and of its corresponding inverse transform,

\displaystyle x_n =\frac{1}{N} \sum_{k=0}^{N-1}{X_k e^{2\pi i k n / N}}.

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 N2, 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 N2 of them to calculate, we can take advantage of the fact that

\displaystyle e^{\varphi i} = \cos \varphi + i \sin \varphi,

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 N2, 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

\displaystyle X_0 = x_0 + x_1, \\ X_1 = x_0 - x_1,

and if N = 4, one can similarly arrive at

\displaystyle X_0 = x_0 + x_1  + x_2 + x_3, \\ X_1 = x_0 - x_2 - i (x_1 - x_3), \\ X_2 = x_0 + x_2 - (x_1+x_3), \\ X_3 = x_0 - x_2 + i (x_1 - x_3).

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…

Advertisements

4 Responses to The Discrete Fourier Transform

  1. Chitiz says:

    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?

  2. Jaime says:

    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

  3. Alessio Treglia says:

    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.

  4. HKS says:

    Hello,

    Very good tutorial.

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

    Thanks and regards,

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: