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, x[s], x[2s], ..., x[(N-1)s]): if N = 1 then X ← x trivial size-1 DFT base case else X[0,...,N/2−1] ← ditfft2(x, N/2, 2s) DFT of (x, 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
above). My first high-level version worked nicely — at least as compared with
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 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
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
- 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
- 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
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 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
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.
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.