Task

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(re2 + im2))   of the complex result.

The classic version is the recursive Cooley–Tukey FFT. Wikipedia has pseudo-code for that.
Further optimizations are possible but not required.




import <Utilities/Complex.sl>;
import <Utilities/Math.sl>;
import <Utilities/Sequence.sl>;
 
fft(x(1)) :=
    let
        n := size(x);
 
        top := fft(x[range(1,n-1,2)]);
        bottom := fft(x[range(2,n,2)]);
 
        d[i] := makeComplex(cos(2.0*pi*i/n), -sin(2.0*pi*i/n)) foreach i within 0...(n / 2 - 1);
 
        z := complexMultiply(d, bottom);
    in
        x when n <= 1
    else
        complexAdd(top,z) ++ complexSubtract(top,z);