Fast Fourier transform

[andrewdavidmackenzie]

[andrewdavidmackenzie]

from https://www.cs.cmu.edu/~scandal/nesl/alg-numerical.html#fft

This is a simple parallel version of the standard sequential fast Fourier transform algorithm. The algorithm does O(n log n) work and has O(log n) depth. In the code we first give a general FFT that works with any commutative ring that has primitive nth roots of unity. In the code we first give a general FFT that works with any commutative ring. As well as the data, this FFT takes as arguments the n different nth-roots of unity (an array w) and the +, * functions on the ring (add, mult). We then specialize the algorithm to a particular commutative ring---the complex numbers.

function fft(a,w,add,mult) = if #a == 1 then a else let r = {fft(b, even_elts(w), add, mult): b in [even_elts(a),odd_elts(a)]} in {add(a, mult(b, w)): a in r[0] ++ r[0]; b in r[1] ++ r[1]; w in w};

function complex_fft(a) = let c = 2.pi/float(#a); w = {cos(cfloat(i)),sin(cfloat(i)) : i in [0:#a]}; add = ((ar,ai),(br,bi)) => (ar+br,ai+bi); mult = ((ar,ai),(br,bi)) => (arbr-aibi,arbi+ai*br); in fft(a,w,add,mult);

complex_fft([(2.,0.),(-1.,1.),(0.,0.),(-1.,-1.)]);