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. $$

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.)