Tuesday, 10 September 2013

A simple Fast Fourier Transform (FFT)

The Fast Fourier Transform is the most important spectral method in computing and lies at the heart of many common algorithms including compressed formats for image, video and sound.

The following is a simple 15-line implementation of the FFT written in F#:

open System.Numerics

let fft a =
  let rec loop n a =
    if n < 2 then Array.init n a else
      let e = loop (n/2) (fun i -> a(2*i))
      let o = loop (n/2) (fun i -> a(2*i + 1))
      let a = Array.create n Complex.Zero
      let f k k' =
        let t = 2. * System.Math.PI * float k / float n
        a.[k] <- complex="" cos="" e.="" k="" o.="" o:p="" sin="" t="">
      for k=0 to n/2 - 1 do
        f k k
      for k=n/2 to n - 1 do
        f k (k - n/2)
  loop (Array.length a) (Array.get a)

Note that this implementation works only for inputs with an integral-power-of-two number of samples. For a professional high-performance solution that handles any size input please buy our F# for Numerics software library.

The following array of complex numbers represents the spectrum of a low-frequency cosine wave:

let zs =
  [|0; 1; 0; 0; 0; 0; 0; 1|]
  |> Array.map (fun x -> Complex(float x, 0.0))

In order to make the output from F# Interactive comprehensible we shall provide our own custom pretty printer for complex numbers. Firstly, we wish to truncate small floating point numbers to zero:

let chop x =
  if abs x < 1e-5 then 0.0 else x

Now we can write a function to convert the chopped real and imaginary components of a complex number into a string, skipping any zeros:

let stringOfComplex (z: Complex) =
  match chop z.Real, chop z.Imaginary with
  | 0.0, 0.0 -> "0"
  | 0.0, y -> sprintf "%fi" y
  | x, 0.0 -> sprintf "%f" x
  | x, y -> sprintf "%f + %fi" x y

This function may be registered as the pretty printed for complex numbers with F# Interactive as follows:

fsi.AddPrinter stringOfComplex

The output of our fft function is then:

fft zs

val it : Complex [] =
  [|2.000000; 1.414214; 0; -1.414214; -2.000000; -1.414214; 0; 1.414214|]

This is one complete cycle of a cosine wave (the lowest possible frequency with this sampling) with amplitude two. Note the appearance of 2×cos(pi/4) = √2.


Unknown said...

hi John,

complex="" cos="" e.="" k="" o.="" o:p="" sin="" t="">
seems like a strange way of defining a complex number, and it doesn't seem to work with only System.Numerics

is there a typo, and are you relying on your FSharp.Numerics for some complex tpyes?


Jon Harrop said...

Hi Daniel,

Google Blogger has mangled my code. I assume I had a less than comparison and it interpreted the rest of the line as an HTML tag with attributes.

I'll try to fix it but Blogger is quite buggy so it may not be possible...