# Scientific computing: Fast Fourier transform

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(N2) 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.$$

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
j <- 2 * j