Simple OCaml implementation of Fast Fourier Transform.

fft.ml

type complex = {re: float; im: float}

let complex_add a b = {re = a.re +. b.re; im = a.im +. b.im}
let complex_sub a b = {re = a.re -. b.re; im = a.im -. b.im}
let complex_mul a b = {re = a.re *. b.re -. a.im *. b.im; im = a.re *. b.im +. a.im *. b.re}
let complex_exp theta = {re = cos theta; im = sin theta}

let rec fft x =
  let n = Array.length x in
  if n <= 1 then x
  else
    let even = Array.init (n / 2) (fun i -> x.(2 * i))
    and odd = Array.init (n / 2) (fun i -> x.(2 * i + 1)) in
    let even = fft even
    and odd = fft odd in
    let t = ref {re = 0.0; im = 0.0} in
    let result = Array.make n {re = 0.0; im = 0.0} in
    for k = 0 to n / 2 - 1 do
      t := complex_mul (complex_exp (-. 2.0 *. Float.pi *. float_of_int k /. float_of_int n)) odd.(k);
      result.(k) <- complex_add even.(k) !t;
      result.(k + n / 2) <- complex_sub even.(k) !t;
    done;
    result

(* example usage *)
let () =
  let signal = [|{re = 1.0; im = 0.0}; {re = 2.0; im = 0.0}; {re = 3.0; im = 0.0}; {re = 4.0; im = 0.0}|] in
  let result = fft signal in
  Array.iter (fun c -> Printf.printf "(%f, %f)\n" c.re c.im) result