Thing practiced:computers

Tools used:Visual Studio 2012 on a Windows machine

(Quick background: The Fourier transform lets us translate any irregular signal into a combination of pure sine waves. The numerical version of this is the discrete Fourier transform (DFT), which is ubiquitous today – it’s responsible for jpegs and mp3s, among other things – thanks to the efficient fast Fourier transform (FFT) family of algorithms.)

Let’s try to implement the FFT in F#. First, here’s the naive DFT algorithm (a direct translation to code of the definition): $$ X_k = \sum_{n=0}^{N-1} {x_n}e^{\frac{-2{\pi}i}{N} kn} \tag{DFT} $$

```
let slow_dft xs =
let N = Array.length xs
[| for k in 0 .. N - 1 ->
let mutable x_k = Complex.Zero
for n = 0 to N - 1 do
let a = 2.0 * Math.PI / float N *
float k * float n
x_k <- x_k + xs.[n] *
Complex((cos a), (sin a))
x_k |]
```

This gives an exact answer, but with poor performance – this has O(N^{2}) time complexity. To improve performance, we first reduce problem scope by assuming that the number of input terms is an integral power of 2, letting us use the radix-2 DIT divide-and-conquer algorithm. This splits the DFT into two DFTs – a DFT of the even-numbered and a DFT of the odd-numbered terms:

$$ X_k = \sum_{m=0}^{\frac{N}{2} -1} {x_{2m}}e^{\frac{-2{\pi}i}{N} k(2m)} + \sum_{m=0}^{\frac{N}{2} -1} {x_{2m+1}}e^{\frac{-2{\pi}i}{N} k(2m+1)} $$

which gets rearranged as: $$ X_k = \left\{ \begin{array} {r@{\quad}l} E_k + e^{\frac{-2{\pi}i}{N} k}{O_k} & \mbox{if} & k < N/2 \\ E_{k-N/2} - e^{\frac{-2{\pi}i}{N} (k-N/2)}{O_{k-N/2}} & \mbox{if} & k \ge N/2 \end{array} \right. $$

This article explains further. Code:

```
let radix2_fft xs =
let bit_rev arr =
let N = Array.length arr
let mutable j = 0
let swap (a : _[]) i j =
let tmp = a.[j]
a.[j] <- a.[i]
a.[i] <- tmp
let rec helper m j =
let m = m / 2
let j = j ^^^ m
if j &&& m = 0 then helper m j else j
for i = 0 to N - 2 do
if i < j then swap arr i j
j <- helper N j
let radix2_helper (xs : Complex[]) N j m =
let t = Math.PI * float m / float j
let a = Complex((cos t), (sin t))
let mutable i = m
while i < N do
let x_i = xs.[i]
let t = a * xs.[i+j]
xs.[i] <- x_i + t
xs.[i+j] <- x_i - t
i <- i + 2 * j
let N = Array.length xs
bit_rev xs
let mutable j = 1
while j < N do
for m = 0 to j - 1 do
radix2_helper xs N j m
j <- 2 * j
```

This is fast – O(n log n) – but only handles signals with lengths of integral powers of 2. To use the radix-2 algorithm for a signal of arbitrary length, we need to increase length to the next smallest power of 2. We can’t zero-pad the signal itself, but if we express the transform as a convolution then pad the convolution we’ll get an exact answer. (My hour’s up, but if you’re interested the convolution algorithm is explained here.)