For my research I needed to implement a discrete Fourier transform in the Clike 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 CooleyTukey radix2 decimationintime 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[(N1)s]):
if N = 1 then
X[0] ← x[0] trivial size1 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) outofplace by a radix2 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 highlevel language (Python) to make sure I get the algorithm itself right, and 2) reduce the use of highlevel features until it is trivial to port to a Clike language.
List handling in Python is always a charm, but the main highlevel feature
I wanted was the handling of complex numbers (notice the −2πi
s
above). My first highlevel 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 

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 oneminute 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 

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 

Testing still showed that the approaches disagreed:
1 2 3 4 5 

1 2 3 4 

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 infft3
.  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 tofft3
(rather than simply givinga
).  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 

by
1 2 3 4 5 6 

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 

However, the goal was using fewer highlevel 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 

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.