Good Times

by James Classen

The BBC has done it again. Sherlock is the best modern day Sherlock Holmes story, Benedict Cumberbatch the best Holmes since Jeremy Brett. Affectations, mannerisms, they’re “spot on.” I would very much suspect that the greatest consulting detective of all time would do much the same as this one, given today’s technology. I’ve only seen “A Study In Pink” and “The Blind Banker” so far, and while they’ve done the sensational stories first, I hope they delve deeper into development of the characters as the series progresses. It’s been too long since I’ve read the books, but I recognized so much of these episodes, and they remind me of “A Study In Scarlet” (obviously), “The Sign of the Four”, and “The Five Orange Pips” along with a hint of “The Adventure of the Dancing Men”. Had her name been Mary Morstan…but perhaps she’ll show up later.

The discrete Fourier transform. Simple enough, sure:

\mathcal{F}_n=\sum_{k=0}^{N-1}{f_ke^{-2\pi ink/N}}

I have no problems with this. Or the conversion to a matrix formula:

\mathcal{F}_n=\begin{bmatrix}f_0 & f_1 & \cdots & f_{N-1}\end{bmatrix}\begin{bmatrix}W^{(0)(0)} & W^{(1)(0)} & \cdots & W^{(N-1)(0)} \\ W^{(0)(1)} & W^{(1)(1)} & \cdots & W^{(N-1)(1)} \\ \vdots & \vdots & \ddots & \vdots \\ W^{(0)(0)} & W^{(1)(N-1)} & \cdots & W^{(N-1)(N-1)}\end{bmatrix}

where W=e^{-2\pi in/N} (so W^N=1, these are roots of unity or de Moivre numbers).

Eventually I found the Danielson-Lanczos lemma, which gives

\mathcal{F}_n=\sum_{k=0}^{N-1}{f_ke^{-2\pi ink/N}}

=\sum_{k=0}^{N/2-1}{e^{-2\pi ikn/(N/2)}f_{2k}}+W^n\sum{k=0}^{N/2-1}{e^{-2\pi ikn/(N/2)}f_{2k+1}}

=\mathcal{F}_n^e+W^n\mathcal{F}_n^o

This is all fine and good. Then it comes to working the FFT, using the Danielson-Lanczos lemma and the Cooley-Tukey algorithm, I’m suddenly stumped. Not the mechanics, but the programming bit. I know I’ve done this successfully before, certainly in college and at least once since then—I had a Python script that used this, first to compute DFTs, then to multiply numbers. Obviously unnecessary with Python supporting arbitrary precision arithmetic built in, but the purpose was (and is) to bolster my own understanding of one of the principles behind APM.

So, University at Buffalo to the rescue. Found some C++ code that gave me what I needed. So here’s (hopefully) a more-or-less plain English version of the method (with an example):

Given a vector (we’ll start small, a size of 8) f=\begin{bmatrix}7 & 3 & 9 & 1 & 3 & 7 & 3 & 5\end{bmatrix}.

If the number of elements in the vector is even, split the odd and even parts of the vector into separate vectors (so I cheated and chose a nice number. But the method doesn’t help if the above isn’t true). We’ll call those f_o and f_e, respectively:

f_o=\begin{bmatrix}3 & 1 & 7 & 5\end{bmatrix}

f_e=\begin{bmatrix}7 & 9 & 3 & 3\end{bmatrix}

These are still divisible by two, so repeat:

f_{oo}=\begin{bmatrix}1 & 5\end{bmatrix}

f_{oe}=\begin{bmatrix}3 & 7\end{bmatrix}

f_{eo}=\begin{bmatrix}9 & 3\end{bmatrix}

f_{ee}=\begin{bmatrix}7 & 3\end{bmatrix}

Still even!

f_{ooo}=\begin{bmatrix}5\end{bmatrix}

f_{ooe}=\begin{bmatrix}1\end{bmatrix}

f_{oeo}=\begin{bmatrix}7\end{bmatrix}

f_{oee}=\begin{bmatrix}3\end{bmatrix}

f_{eoo}=\begin{bmatrix}3\end{bmatrix}

f_{eoe}=\begin{bmatrix}9\end{bmatrix}

f_{eeo}=\begin{bmatrix}3\end{bmatrix}

f_{eee}=\begin{bmatrix}7\end{bmatrix}

Now the size of these vectors is odd, so the only choice is to do a regular discrete transform. But this turns out to be nice and simple. Using superscripts to denote which of the above vectors I’m refering to, since

\mathcal{F}^{ooo}=\sum_{k=0}^{N-1}{f_ke^{-2\pi ink/N}}

we can plug in each of the above into this formula:

\mathcal{F}^{ooo}=\sum_{k=0}^{0}{f_0e^{-2\pi i\times 0 \times k/1}}

which nicely reduces to

\mathcal{F}^{ooo}=f_0

But this is only the discrete transform of the one-element vector—these need to be recombined to get the DFT of the entire thing.

So if you reduce the vectors to size 1, the DFT of that is trivial. Then we need to combine them. So let’s start that combination one level higher, with the two-element vectors.

f_{oo}=\begin{bmatrix}1 & 5\end{bmatrix}

From the Danielson-Lanczos lemma,

\mathcal{F}_n^{oo}=\mathcal{F}_{n\mod(N/2)}^{ooe}+W^n\mathcal{F}_{n\mod(N/2)}^{ooo}

Therefore,

\mathcal{F}_0^{oo}=\mathcal{F}_0^{ooe}+W^0\mathcal{F}_0^{ooo}

\mathcal{F}_1^{oo}=\mathcal{F}_0^{ooe}+W^1\mathcal{F}_0^{ooo}

In this case, W=e^{-2\pi i/2}=e^{-\pi i}=-1, so

\mathcal{F}^{oo}=\begin{bmatrix}1+5 & 1-5\end{bmatrix}=\begin{bmatrix}6 & -4\end{bmatrix}

So for a two-element vector, \mathcal{F}=\begin{bmatrix}\mathcal{F}^{e}+\mathcal{F}^{o} & \mathcal{F}^{e}-\mathcal{F}^{o}\end{bmatrix}

“Solving” the others in a similar manner, we get

\mathcal{F}^{oo}=\begin{bmatrix}6 & -4\end{bmatrix}

\mathcal{F}^{oe}=\begin{bmatrix}10 & -4\end{bmatrix}

\mathcal{F}^{eo}=\begin{bmatrix}12 & 6\end{bmatrix}

\mathcal{F}^{ee}=\begin{bmatrix}10 & 4\end{bmatrix}

Of course, each step gets a bit more difficult. Thanks to the DLL, we know

\mathcal{F}_n^{o}=\mathcal{F}_{n\mod(N/2)}^{oe}+W^n\mathcal{F}_{n\mod(N/2)}^{oo}

so let’s take it one step at a time.

First, W=e^{-2\pi i/4}=-i

So

\mathcal{F}_0^{o}=\mathcal{F}_0^{oe}+W^0\mathcal{F}_0^{oo}

=10+6=16

\mathcal{F}_1^{o}=\mathcal{F}_1^{oe}+W^1\mathcal{F}_1^{oo}

=-4+-i\times -4=-4+4i

\mathcal{F}_2^{o}=\mathcal{F}_0^{oe}+W^2\mathcal{F}_0^{oo}

=10+-1\times 6=4

\mathcal{F}_3^{o}=\mathcal{F}_1^{oe}+W^3\mathcal{F}_1^{oo}

=-4+i\times -4=-4-4i

Therefore,

\mathcal{F}^{o}=\begin{bmatrix}16 & -4+4i & 4 & -4-4i\end{bmatrix}

Do the same thing to find \mathcal{F}^{e}:

\mathcal{F}^{e}=\begin{bmatrix}22 & 4-6i & -2 & 4+6i\end{bmatrix}

Putting it all together,

W=e^{-2\pi i/8}={{1}\over{\sqrt{2}}}(1-i)

\mathcal{F}_n=\mathcal{F}_{n\mod(N/4)}^{e}+W^n\mathcal{F}_{n\mod(N/4)}^{o}

\mathcal{F}_0=\mathcal{F}_0^e+W^0\mathcal{F}_0^o

=22+1\times 16=38

\mathcal{F}_1=\mathcal{F}_1^e+W^1\mathcal{F}_1^o

=(-4+4i)+{{1}\over{\sqrt{2}}}(1-i)(4-6i)=(-4-\sqrt{2})+(4-5\sqrt{2})i

\mathcal{F}_2=\mathcal{F}_2^e+W^2\mathcal{F}_2^o

=-2+-i\times4=-2-4i

\mathcal{F}_3=\mathcal{F}_3^e+W^3\mathcal{F}_3^o

=(-4-4i)+{{1}\over{\sqrt{2}}}(-1-i)(4+6i)=(-4+\sqrt{2})+(-4-5\sqrt{2})i

\mathcal{F}_4=\mathcal{F}_0^e+W^4\mathcal{F}_0^o

=22+-1\times 16=6

\mathcal{F}_5=\mathcal{F}_1^e+W^5\mathcal{F}_1^o

=(-4+4i)+{{1}\over{\sqrt{2}}}(-1+i)(4-6i)=(-4+\sqrt{2})+(4+5\sqrt{2})i

\mathcal{F}_6=\mathcal{F}_2^e+W^6\mathcal{F}_2^o

=-2+i\times4=-2+4i

\mathcal{F}_7=\mathcal{F}_3^e+W^7\mathcal{F}_3^o

=(-4-4i)+{{1}\over{\sqrt{2}}}(1+i)(4+6i)=(-4-\sqrt{2})+(-4+5\sqrt{2})i

Thus

\mathcal{F}=\begin{bmatrix}38 & (-4-\sqrt{2})+(4-5\sqrt{2})i & -2-4i & (-4+\sqrt{2})+(-4-5\sqrt{2})i & 6 & (-4+\sqrt{2})+(4+5\sqrt{2})i & -2+4i & (-4-\sqrt{2})+(-4+5\sqrt{2})i\end{bmatrix}

Seems like a lot of work, I know. But here’s the thing: normally, using the straight definition, you’d be doing 64 multiplications, 64 sums, and 64 exponentiations (O(N^2). We did 24 of each (O(N\log_2{N}). And that’s without using a trick: if all f_k are real, \mathcal{F}_{N-n}=\overline{\mathcal{F}_{n}}

Yes, it’s still a lot of work, but as the vectors get longer, these savings adds up fast:

Vector Size Naive Operations FFT operations
1 1 0
2 4 2
4 16 8
8 64 24
16 256 64
32 1024 160

And once you put it in a computer program (like the following listing), massive data sets become easy to analyze.