For my research I needed to implement a discrete Fourier transform in the C-like programming language of the UPPAAL model checker. One option would be to port an existing C library such as FFTW, but for a proof of concept, porting a large, powerful library would be much slower and less instructive than implementing a simpler algorithm myself. Programming rant coming up.

I went with the Cooley-Tukey radix-2 decimation-in-time fast Fourier transform, which, despite the name, is one of the simpler choices. Pseudocode adapted from Wikipedia:

X[0,...,N−1] ← ditfft2(x, N, s):             DFT of (x[0], x[s], x[2s], ..., x[(N-1)s]):
    if N = 1 then
        X[0] ← x[0]                                  trivial size-1 DFT base case
    else
        X[0,...,N/2−1] ← ditfft2(x, N/2, 2s)         DFT of (x[0], x[2s], x[4s], ...)
        X[N/2,...,N−1] ← ditfft2(x+s, N/2, 2s)       DFT of (x[s], x[s+2s], x[s+4s], ...)
        for k = 0 to N/2−1                           combine DFTs of two halves into full DFT:
            t ← X[k]
            X[k] ← t + exp(−2πi k/N) X[k+N/2]
            X[k+N/2] ← t − exp(−2πi k/N) X[k+N/2]
        endfor
    endif

Here, ditfft2(x,N,1), computes X=DFT(x) out-of-place by a radix-2 DIT FFT,
where N is an integer power of 2 and s=1 is the stride of the input
x array. x+s denotes the array starting with xs.

I decided to handle the problem in two stages: 1) Implement the algorithm in a high-level language (Python) to make sure I get the algorithm itself right, and 2) reduce the use of high-level features until it is trivial to port to a C-like language.

List handling in Python is always a charm, but the main high-level feature I wanted was the handling of complex numbers (notice the −2πis above). My first high-level version worked nicely — at least as compared with the fft function of R with a few sample inputs. Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

import math
EVEN = 0
ODD = 1

def fft1(a):
    n = len(a)
    if n == 1:
        return a
    c = fft1(a[EVEN::2]) + fft1(a[ODD::2])
    for k in range(n / 2):
        t = c[k]
        c[k]         = t + pow(math.e, -2j * math.pi * k / n) * c[k + n / 2]
        c[k + n / 2] = t - pow(math.e, -2j * math.pi * k / n) * c[k + n / 2]
    return c

Next step was to handle the complex numbers manually, so they can be represented by two doubles. At first I tried to handle the general case of raising a real to a complex, but that was harder than I thought it would be. Luckily my colleague Louise could derive it for me:

 a + bi      a                   a
x        =  x  cos(b ln(x)) + i x  sin(b ln(x))

(and I wish adding MathJax (LaTeX) to Octopress was a one-minute job, but it isn’t.)

However, if I had remembered that in my case x = e and a = 0, it would have been simpler. It even has a name, Euler’s formula, and it is mindboggling in all its simple beauty:

 bi
e    =  cos(b) + i sin(b)

At least we both learned something about complex numbers!

The next (buggy) version of my code was this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

def fft2(a_re, a_im):
    n = len(a_re)
    assert n == len(a_im)
    if n == 1:
        return (a_re, a_im)
    c_re, c_im = fft2(a_re[EVEN::2], a_im[EVEN::2])
    c_re2, c_im2 = fft2(a_re[ODD::2], a_im[ODD::2])
    c_re += c_re2
    c_im += c_im2
    for k in range(n / 2):
        t_re, t_im = c_re[k], c_im[k]

        c_re[k] = t_re + math.cos(-2 * math.pi * k / n) * c_re[k + n / 2]
        c_im[k] = t_im + math.sin(-2 * math.pi * k / n) * c_im[k + n / 2]
        c_re[k + n / 2] = t_re - math.cos(-2 * math.pi * k / n) * c_re[k + n / 2]
        c_im[k + n / 2] = t_im - math.sin(-2 * math.pi * k / n) * c_im[k + n / 2]

    return (c_re, c_im)

It didn’t give the same output as the first one, and I couldn’t see why, so I put both versions in one to debug them in lockstep:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

def fft3(a, a_re, a_im):
    n = len(a)
    assert n == len(a_re) == len(a_im)
    if n == 1:
        return (a, a_re, a_im)
    c, c_re, c_im = fft3(a[EVEN::2], a_re[EVEN::2], a_im[EVEN::2])
    c2, c_re2, c_im2 = fft3(a[ODD::2], a_re[ODD::2], a_im[ODD::2])
    c += c2
    del c2
    c_re += c_re2
    del c_re2
    c_im += c_im2
    del c_im2
    for k in range(n / 2):
        t, t_re, t_im = c[k], c_re[k], c_im[k]

        c[k]         = t + pow(math.e, -2j * math.pi * k / n) * c[k + n / 2]
        c[k + n / 2] = t - pow(math.e, -2j * math.pi * k / n) * c[k + n / 2]

        c_re[k]         = t_re + math.cos(-2 * math.pi * k / n) * c_re[k + n / 2]
        c_im[k]         = t_im + math.sin(-2 * math.pi * k / n) * c_im[k + n / 2]
        c_re[k + n / 2] = t_re - math.cos(-2 * math.pi * k / n) * c_re[k + n / 2]
        c_im[k + n / 2] = t_im - math.sin(-2 * math.pi * k / n) * c_im[k + n / 2]

    return (c, c_re, c_im)

Testing still showed that the approaches disagreed:

1
2
3
4
5

a = [1, 2, 1, 0]
x1, x2_re, x2_im = fft3(map(complex, a), a, [0] * len(a))
x2 = [complex(re, im) for re, im in zip(x2_re, x2_im)]
for x, y in zip(x1, x2):
    print x, y, x - y, '!' * 10 if x - y != 0 else ''

1
2
3
4

(4+0j) (4+0j) 0j
(1.22464679915e-16-2j) (1.22464679915e-16+0j) -2j !!!!!!!!!!
0j 0j 0j
(-1.22464679915e-16+2j) (-1.22464679915e-16+0j) 2j !!!!!!!!!!

Stepping through the code in Winpdb, I found three things:

• Debugging is easier if you clean up temp variables in the namespace, hence the del statements in fft3.
• Debugging is easier if the complex numbers are represented as complex numbers all the time, hence the test code giving map(complex, a) as the first argument to fft3 (rather than simply giving a).
• Line 21 was the first to show the bug (for the given input): I had been so focused on complex exponentiation that I had forgot all about complex multiplication! Which is btw
  • (a + bi) (c + di) = (ac - bd) + (bc + ad) i

Now I (naively hoped I) was done. Replacing

1
2
3
4

c_re[k]         = t_re + math.cos(-2 * math.pi * k / n) * c_re[k + n / 2]
c_im[k]         = t_im + math.sin(-2 * math.pi * k / n) * c_im[k + n / 2]
c_re[k + n / 2] = t_re - math.cos(-2 * math.pi * k / n) * c_re[k + n / 2]
c_im[k + n / 2] = t_im - math.sin(-2 * math.pi * k / n) * c_im[k + n / 2]

by

1
2
3
4
5
6

p_re = math.cos(-2 * math.pi * k / n)
p_im = math.sin(-2 * math.pi * k / n)
c_re[k]         = t_re +  (p_re * c_re[k + n / 2]) - (p_im * c_im[k + n / 2])
c_im[k]         = t_im +  (p_im * c_re[k + n / 2]) + (p_re * c_im[k + n / 2])
c_re[k + n / 2] = t_re - ((p_re * c_re[k + n / 2]) - (p_im * c_im[k + n / 2]))
c_im[k + n / 2] = t_im - ((p_im * c_re[k + n / 2]) + (p_re * c_im[k + n / 2]))

still didn’t work. But that is just a silly bug: Reading c_re[k + n / 2] in line 6 after it was overwritten in line 5. Python has a really nice solution for problems like this: Do both assignments at once! Instead of lines 5 and 6:

1
2
3
4

c_re[k + n / 2], c_im[k + n / 2] = (
    t_re - ((p_re * c_re[k + n / 2]) - (p_im * c_im[k + n / 2])),
    t_im - ((p_im * c_re[k + n / 2]) + (p_re * c_im[k + n / 2])),
)

However, the goal was using fewer high-level features, so my final solution factors out the reused parts of the computation, which probably makes the code simpler in this case anyway:

1
2
3
4
5
6
7
8

p_re = math.cos(-2 * math.pi * k / n)
p_im = math.sin(-2 * math.pi * k / n)
mult_re = (p_re * c_re[k + n / 2]) - (p_im * c_im[k + n / 2])
mult_im = (p_im * c_re[k + n / 2]) + (p_re * c_im[k + n / 2])
c_re[k]         = t_re + mult_re
c_im[k]         = t_im + mult_im
c_re[k + n / 2] = t_re - mult_re
c_im[k + n / 2] = t_im - mult_im

So, after a few detours I got my “manually” computed complex numbers, and learned a bit more about how complex numbers and FFTs work. Also, debugging with a debugger (rather than print statements) doesn’t mean you don’t need to modify your code – in my case, merging both functions made the debugging a lot easier.

Finally, the Python math library doesn’t support complex numbers! But there is a good reason for that, as explained in the cmath (complex math) library:

The reason for having two modules is that some users aren’t interested in complex numbers, and perhaps don’t even know what they are. They would rather have math.sqrt(-1) raise an exception than return a complex number.