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.
You are encouraged to solve this task according to the task description, using any language you may know.
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]
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)
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)) >>> 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