Fast Fourier transform: Difference between revisions

Content added Content deleted

Inline

Revision as of 12:37, 22 July 2011

The purpose of this task is to calculate the FFT (Fast Fourier Transform) of an input sequence. The most general case allows for complex numbers at the input and results in a sequence of equal length, again of complex numbers. If you need to restrict yourself to real numbers the output should be the magnitude (i.e. sqrt(re²+im²)) of the complex result. The classic version is the recursive Cooley–Tukey FFT. Wikipedia has pseudocode for that. Further optimizations are possible but not required.

Ada

<lang Ada> with Ada.Numerics.Generic_Complex_Arrays;

generic

  with package Complex_Arrays is
     new Ada.Numerics.Generic_Complex_Arrays (<>);
  use Complex_Arrays;

function Generic_FFT (X : Complex_Vector) return Complex_Vector;

</lang>

with Ada.Numerics; with Ada.Numerics.Generic_Complex_Elementary_Functions;

function Generic_FFT (X : Complex_Vector) return Complex_Vector is

  package Complex_Elementary_Functions is
     new Ada.Numerics.Generic_Complex_Elementary_Functions
       (Complex_Arrays.Complex_Types);

  use Ada.Numerics;
  use Complex_Elementary_Functions;
  use Complex_Arrays.Complex_Types;
  
  function FFT (X : Complex_Vector; N, S : Positive)
     return Complex_Vector is
  begin
     if N = 1 then
        return (1..1 => X (X'First));
     else
        declare
           F : constant Complex  := exp (Pi * j / Real_Arrays.Real (N/2));
           Even : Complex_Vector := FFT (X, N/2, 2*S);
           Odd  : Complex_Vector := FFT (X (X'First + S..X'Last), N/2, 2*S);
        begin
           for K in 0..N/2 - 1 loop
              declare
                 T : constant Complex := Odd (Odd'First + K) / F ** K;
              begin
                 Odd  (Odd'First  + K) := Even (Even'First + K) - T;
                 Even (Even'First + K) := Even (Even'First + K) + T;
              end;
           end loop;
           return Even & Odd;
        end;
     end if;
  end FFT;

begin

  return FFT (X, X'Length, 1);

end Generic_FFT; </lang>

with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays; with Ada.Complex_Text_IO; use Ada.Complex_Text_IO; with Ada.Text_IO; use Ada.Text_IO;

with Ada.Numerics.Complex_Elementary_Functions; with Generic_FFT;

procedure Example is

  function FFT is new Generic_FFT (Ada.Numerics.Complex_Arrays);
  X : Complex_Vector := (1..4 => (1.0, 0.0), 5..8 => (0.0, 0.0));
  Y : Complex_Vector := FFT (X);

begin

  Put_Line ("       X              FFT X ");
  for I in Y'Range loop
     Put (X (I - Y'First + X'First), Aft => 3, Exp => 0);
     Put (" ");
     Put (Y (I), Aft => 3, Exp => 0);
     New_Line;
  end loop;

end; </lang>

C

Inplace FFT with O(n) memory usage. Note: array size is assumed to be power of 2 and not checked by code; you can just pad it with 0 otherwise. Also, complex is C99 standard.<lang C>#include <stdio.h>

include <math.h>
include <complex.h>

double PI; typedef double complex cplx;

void _fft(cplx buf[], cplx out[], int n, int step) { if (step < n) { _fft(out, buf, n, step * 2); _fft(out + step, buf + step, n, step * 2);

for (int i = 0; i < n; i += 2 * step) { cplx t = cexp(-I * PI * i / n) * out[i + step]; buf[i / 2] = out[i] + t; buf[(i + n)/2] = out[i] - t; } } }

void fft(cplx buf[], int n) { cplx out[n]; for (int i = 0; i < n; i++) out[i] = buf[i];

_fft(buf, out, n, 1); }

int main() { PI = atan2(1, 1) * 4; cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};

void show(char * s) { printf(s); for (int i = 0; i < 8; i++) if (!cimag(buf[i])) printf("%g ", creal(buf[i])); else printf("(%g, %g) ", creal(buf[i]), cimag(buf[i])); }

show("Data: "); fft(buf, 8); show("\nFFT : ");

return 0; }</lang>Output<lang>Data: 1 1 1 1 0 0 0 0 FFT : 4 (1, -2.41421) 0 (1, -0.414214) 0 (1, 0.414214) 0 (1, 2.41421)</lang>

C++

<lang cpp>#include <complex>

include <iostream>
include <valarray>

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex; typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place) void fft(CArray& x) {

   const size_t N = x.size();
   if (N <= 1) return;

   // divide
   CArray even = x[std::slice(0, N/2, 2)];
   CArray  odd = x[std::slice(1, N/2, 2)];

   // conquer
   fft(even);
   fft(odd);

   // combine
   for (size_t k = 0; k < N/2; ++k)
   {
       Complex t = exp(Complex(0, -2 * PI * k / N)) * odd[k];
       x[k    ] = even[k] + t;
       x[k+N/2] = even[k] - t;
   }

}

int main() {

   const Complex test[] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
   CArray data(test, 8);

   fft(data);

   for (int i = 0; i < 8; ++i)
   {
       std::cout << data[i] << "\n";
   }

}</lang> Output:

(4,0)
(1,-2.41421)
(0,0)
(1,-0.414214)
(0,0)
(1,0.414214)
(0,0)
(1,2.41421)

D

<lang d>import std.stdio, std.math, std.algorithm, std.range;

auto fft(creal[] x) {

 int N = x.length;
 if (N <= 1) return x;
 auto ev = fft(array(stride(x, 2)));
 auto od = fft(array(stride(x[1 .. $], 2)));
 auto l = map!((k){return ev[k]+expi(-2*PI*k/N)*od[k];})(iota(N/2));
 auto r = map!((k){return ev[k]-expi(-2*PI*k/N)*od[k];})(iota(N/2));
 return array(chain(l, r));

}

void main() {

 writeln(fft([1.0L+0i, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]));

}</lang>

Output:

[4+0i, 1+-2.41421i, 0+0i, 1+-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

GAP

<lang gap># Here an implementation with no optimization (O(n^2)).

In GAP, E(n) = exp(2*i*pi/n), a primitive root of the unity.

Fourier := function(a) local n, z; n := Size(a); z := E(n); return List([0 .. n - 1], k -> Sum([0 .. n - 1], j -> a[j + 1]*z^(-k*j))); end;

InverseFourier := function(a) local n, z; n := Size(a); z := E(n); return List([0 .. n - 1], k -> Sum([0 .. n - 1], j -> a[j + 1]*z^(k*j)))/n; end;

Fourier([1, 1, 1, 1, 0, 0, 0, 0]);

[ 4, 1-E(8)-E(8)^2-E(8)^3, 0, 1-E(8)+E(8)^2-E(8)^3,
0, 1+E(8)-E(8)^2+E(8)^3, 0, 1+E(8)+E(8)^2+E(8)^3 ]

InverseFourier(last);

[ 1, 1, 1, 1, 0, 0, 0, 0 ]</lang>

Go

<lang go>package main

import (

   "fmt"
   "math"
   "cmath"

)

func ditfft2(x []float64, y []complex128, n, s int) {

   if n == 1 {
       y[0] = complex(x[0], 0)
       return
   }
   ditfft2(x, y, n/2, 2*s)
   ditfft2(x[s:], y[n/2:], n/2, 2*s)
   for k := 0; k < n/2; k++ {
       tf := cmath.Rect(1, -2*math.Pi*float64(k)/float64(n)) * y[k+n/2]
       y[k], y[k+n/2] = y[k]+tf, y[k]-tf
   }

}

func main() {

   x := []float64{1, 1, 1, 1, 0, 0, 0, 0}
   y := make([]complex128, len(x))
   ditfft2(x, y, len(x), 1)
   for _, c := range y {
       fmt.Printf("%8.4f\n", c)
   }

}</lang> Output:

(  4.0000 +0.0000i)
(  1.0000 -2.4142i)
(  0.0000 +0.0000i)
(  1.0000 -0.4142i)
(  0.0000 +0.0000i)
(  1.0000 +0.4142i)
(  0.0000 +0.0000i)
(  1.0000 +2.4142i)

Haskell

<lang haskell>import Control.Monad import Data.Array import Data.Complex

-- Cooley-Tukey fft [] = [] fft [x] = [x] fft xs = zipWith (+) ys ts ++ zipWith (-) ys ts

   where n = length xs
         xs' = listArray (1,n) xs
         ys = fft $ map (xs'!) [1,3..n]
         zs = fft $ map (xs'!) [2,4..n]
         ts = zipWith (\z k -> (exp' k n) * z) zs [0..(n `div` 2)-1]
         exp' k n = exp (-2 * pi * (0 :+ 1) * (fromIntegral k) / (fromIntegral n))

main =

   do let res = fft $ map (:+ 0) [1,1,1,1,0,0,0,0]
      forM_ res print

</lang>

And the output:

4.0 :+ 0.0
1.0 :+ (-2.414213562373095)
0.0 :+ 0.0
1.0 :+ (-0.4142135623730949)
0.0 :+ 0.0
0.9999999999999999 :+ 0.4142135623730949
0.0 :+ 0.0
0.9999999999999997 :+ 2.414213562373095

J

Based on j:Essays/FFT, with some simplifications, sacrificing accuracy, optimizations and convenience not visible here, for clarity:

<lang j>cube =: ($~ q:@#) :. , rou =: ^@j.@o.@(% #)@i.@-: NB. roots of unity floop =: 4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.' fft =: ] floop&.cube rou@#</lang>

Example:

<lang j> require'printf'

  fmt =: [:, sprintf~&'%7.3f'"0
  ('wave:',:'fft:'),.fmt"1 (,: |@fft) 1 o. 2p1*3r16 * i.16

wave: 0.000 0.924 0.707 -0.383 -1.000 -0.383 0.707 0.924 0.000 -0.924 -0.707 0.383 1.000 0.383 -0.707 -0.924 fft: 0.000 0.000 0.000 8.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 8.000 0.000 0.000</lang>

Note that sprintf does not support complex arguments, so we only display the magnitude of the fft here.

Mathematica

Mathematica has a built-in FFT function which uses a proprietary algorithm developed at Wolfram Research. It also has an option to tune the algorithm for specific applications. The options shown below, while not default, produce output that is consistent with most other FFT routines.

<lang Mathematica> Fourier[{1,1,1,1,0,0,0,0}, FourierParameters->{1,-1}] </lang>

Output:

{4. + 0. I, 1. - 2.4142136 I, 0. + 0. I, 1. - 0.41421356 I, 0. + 0. I, 1. + 0.41421356 I, 0. + 0. I, 1. + 2.4142136 I}

OCaml

This is a simple implementation of the Cooley-Tukey pseudo-code <lang OCaml>open Complex;;

let negtwopi = -4.0 *. acos 0.0;;

let fac k n = polar 1.0 (negtwopi*.(float k)/.(float n));;

let merge l r n =

  let f (k,t) x = (succ k, (mul (fac k n) x)::t) in
  let z = List.rev (snd (List.fold_left f (0,[]) r)) in
  (List.map2 add l z) @ (List.map2 sub l z);;

let fft lst =

  let rec ditfft2 a n s =
     if n = 1 then [List.nth lst a] else
     let o = ditfft2 a (n/2) (2*s) in
     let e = ditfft2 (a+s) (n/2) (2*s) in
     merge o e n in
  ditfft2 0 (List.length lst) 1;;

let show l = (List.iter (fun x -> Printf.printf "(%f %f) " x.re x.im) l; print_newline ()) in let indata = [one;one;one;one;zero;zero;zero;zero] in begin show indata; show (fft indata); end</lang>

Output:

(1.000000 0.000000) (1.000000 0.000000) (1.000000 0.000000) (1.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) 
(4.000000 0.000000) (1.000000 -2.414214) (0.000000 0.000000) (1.000000 -0.414214) (0.000000 0.000000) (1.000000 0.414214) (0.000000 0.000000) (1.000000 2.414214)

PARI/GP

Naive implementation, using the same testcase as Ada: <lang parigp>FFT(v)=my(t=-2*Pi*I/#v,tt);vector(#v,k,tt=t*(k-1);sum(n=0,#v-1,v[n+1]*exp(tt*n))); FFT([1,1,1,1,0,0,0,0])</lang> Output:

[4.0000000000000000000000000000000000000, 1.0000000000000000000000000000000000000 - 2.4142135623730950488016887242096980786*I, 0.E-37 + 0.E-38*I, 1.0000000000000000000000000000000000000 - 0.41421356237309504880168872420969807856*I, 0.E-38 + 0.E-37*I, 0.99999999999999999999999999999999999997 + 0.41421356237309504880168872420969807860*I, 4.701977403289150032 E-38 + 0.E-38*I, 0.99999999999999999999999999999999999991 + 2.4142135623730950488016887242096980785*I]

Perl 6

<lang perl6>sub fft {

   return @_ if @_ == 1;
   my @evn = fft( @_[0,2...^* >= @_] );
   my @odd = fft( @_[1,3...^* >= @_] );
   my $twd = 2i * pi / @_; # twiddle factor
   @odd  »*=« (^@odd).map(* * $twd)».exp;
   return @evn »+« @odd, @evn »-« @odd;

}

my @seq = ^16; my $cycles = 3; my @wave = (@seq »*» (2*pi / @seq * $cycles))».sin; say "wave: ", @wave.fmt("%7.3f");

say "fft: ", fft(@wave)».abs.fmt("%7.3f");</lang>

Output:

wave:   0.000   0.924   0.707  -0.383  -1.000  -0.383   0.707   0.924   0.000  -0.924  -0.707   0.383   1.000   0.383  -0.707  -0.924
fft:    0.000   0.000   0.000   8.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   8.000   0.000   0.000

Python

<lang python>from cmath import exp, pi

def fft(x):

   N = len(x)
   if N <= 1: return x
   even = fft(x[0::2])
   odd =  fft(x[1::2])
   return [even[k] + exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)] + \
          [even[k] - exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]

print fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])</lang>

Output:

[(4+0j), (1-2.4142135623730949j), 0j, (1-0.41421356237309492j), 0j, (0.99999999999999989+0.41421356237309492j), 0j, (0.99999999999999967+2.4142135623730949j)]

Using module numpy

<lang python>>>> from numpy.fft import fft >>> from numpy import array >>> a = array((0.0, 0.924, 0.707, -0.383, -1.0, -0.383, 0.707, 0.924, 0.0, -0.924, -0.707, 0.383, 1.0, 0.383, -0.707, -0.924)) >>> print( ' '.join("%5.3f" % abs(f) for f in fft(a)) ) 0.000 0.001 0.000 8.001 0.000 0.001 0.000 0.001 0.000 0.001 0.000 0.001 0.000 8.001 0.000 0.001</lang>

R

The function "fft" is readily available in R <lang R>fft(c(1,1,1,1,0,0,0,0))</lang> Output:

4+0.000000i 1-2.414214i 0+0.000000i 1-0.414214i 0+0.000000i 1+0.414214i 0+0.000000i 1+2.414214i

Tcl

Library: Tcllib (Package: math::constants)

Library: Tcllib (Package: math::fourier)

<lang tcl>package require math::constants package require math::fourier

math::constants::constants pi

Helper functions

proc wave {samples cycles} {

   global pi
   set wave {}
   set factor [expr {2*$pi * $cycles / $samples}]
   for {set i 0} {$i < $samples} {incr i} {

lappend wave [expr {sin($factor * $i)}]

   }
   return $wave

} proc printwave {waveName {format "%7.3f"}} {

   upvar 1 $waveName wave
   set out [format "%-6s" ${waveName}:]
   foreach value $wave {

append out [format $format $value]

   }
   puts $out

} proc waveMagnitude {wave} {

   set out {}
   foreach value $wave {

lassign $value re im lappend out [expr {hypot($re, $im)}]

   }
   return $out

}

set wave [wave 16 3] printwave wave

Uses FFT if input length is power of 2, and a less efficient algorithm otherwise

set fft [math::fourier::dft $wave]

Convert to magnitudes for printing

set fft2 [waveMagnitude $fft] printwave fft2</lang> Output:

wave:   0.000  0.924  0.707 -0.383 -1.000 -0.383  0.707  0.924  0.000 -0.924 -0.707  0.383  1.000  0.383 -0.707 -0.924
fft2:   0.000  0.000  0.000  8.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  8.000  0.000  0.000

Ursala

The fftw library is callable from Ursala using the syntax ..u_fw_dft for a one dimensional forward discrete Fourier transform operating on a list of complex numbers. Ordinarily the results are scaled so that the forward and reverse transforms are inverses of each other, but additional scaling can be performed as shown below to conform to convention. <lang ursala>#import nat

import flo

f = <1+0j,1+0j,1+0j,1+0j,0+0j,0+0j,0+0j,0+0j> # complex sequence of 4 1's and 4 0's

g = c..mul^*D(sqrt+ float+ length,..u_fw_dft) f # its fft

cast %jLW

t = (f,g)</lang> output:

(
   <
      1.000e+00+0.000e+00j,
      1.000e+00+0.000e+00j,
      1.000e+00+0.000e+00j,
      1.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j>,
   <
      4.000e+00+0.000e+00j,
      1.000e+00-2.414e+00j,
      0.000e+00+0.000e+00j,
      1.000e+00-4.142e-01j,
      0.000e+00+0.000e+00j,
      1.000e+00+4.142e-01j,
      0.000e+00+0.000e+00j,
      1.000e+00+2.414e+00j>)