Fast Fourier transform

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.

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));
 creal[] l, r;
 foreach (k; 0 .. N/2) {
   l ~= ev[k] + expi(-2 * PI * k / N) * od[k];
   r ~= ev[k] - expi(-2 * PI * k / N) * od[k];
 }
 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]

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.

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)) >>> fft(a) array([ 0.00000000e+00 +0.00000000e+00j,

       -7.77156117e-16 +1.28750059e-03j,
        0.00000000e+00 +0.00000000e+00j,
       -7.66053887e-15 -8.00062775e+00j,
        0.00000000e+00 +0.00000000e+00j,
       -8.88178420e-16 -1.23179335e-03j,
        0.00000000e+00 +0.00000000e+00j,
        2.33146835e-15 +6.83454981e-04j,
        0.00000000e+00 +0.00000000e+00j,
       -7.77156117e-16 -6.83454981e-04j,
        0.00000000e+00 +0.00000000e+00j,
        2.99760217e-15 +1.23179335e-03j,
        0.00000000e+00 +0.00000000e+00j,
        2.44249065e-15 +8.00062775e+00j,
        0.00000000e+00 +0.00000000e+00j,   2.33146835e-15 -1.28750059e-03j])</lang>

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