Fast Fourier Transform

fft.hpp

#include <artery-math.h>
using namespace artery;

// Fast Fourier Transform implementation by Viktor Chlumsky
template <typename T, typename U>
void fft(Complex<T> *dst, const U *src, int length, int stride = 1) {
    if (!(length >>= 1)) {
        *dst = *src;
        return;
    }
    fft(dst, src, length, stride<<1);
    fft(dst+length, src+stride, length, stride<<1);
    Complex<T> w = exp(-PI/length*Complex<T>::I);
    for (Complex<T> q = 1, *mid = dst+length; length--; q *= w) {
        Complex<T> k = q**mid;
        *mid++ = *dst-k;
        *dst++ += k;
    }
}