Fast Fourier transform
You are encouraged to solve this task according to the task description, using any language you may know.
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
Ada has build-in complex numbers of different accuracy (see also Arithmetic/Complex. The solution uses these types. It is thus generic, parametrized by the given complex type: <lang Ada> with Ada.Numerics.Generic_Complex_Arrays; with Ada.Numerics.Generic_Complex_Elementary_Functions;
generic
with package Arrays is new Ada.Numerics.Generic_Complex_Arrays (<>); use Arrays; with package Elementary_Functions is new Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Types);
function Generic_FFT (X : Complex_Vector) return Complex_Vector; </lang> Generic_FFT is a generic function which takes as the parameters instances of generic packages Ada.Numerics.Generic_Complex_Arrays (ARM G.3.2) for complex vectors and Ada.Numerics.Generic_Complex_Elementary_Functions (ARM G.1.2) for elementary functions like exp. The implementation is: <lang Ada> with Ada.Numerics; use Ada.Numerics;
function Generic_FFT (X : Complex_Vector) return Complex_Vector is
use Complex_Types; use Elementary_Functions; 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> The following is an example of use: <lang Ada> 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 Test_FFT is
function FFT is new Generic_FFT ( Ada.Numerics.Complex_Arrays, Ada.Numerics.Complex_Elementary_Functions ); 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 Test_FFT; </lang> Sample output:
X FFT X ( 1.000, 0.000) ( 4.000, 0.000) ( 1.000, 0.000) ( 1.000,-2.414) ( 1.000, 0.000) ( 0.000, 0.000) ( 1.000, 0.000) ( 1.000,-0.414) ( 0.000, 0.000) ( 0.000, 0.000) ( 0.000, 0.000) ( 1.000, 0.414) ( 0.000, 0.000) ( 0.000, 0.000) ( 0.000, 0.000) ( 1.000, 2.414)
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
<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>)