C++ Code for Fast Fourier Transform (FFT)

            
#include <iostream>
#include <complex>
#include <valarray>
#include <cmath>

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley-Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x) {
    const size_t N = x.size();
    if (N <= 1) return;

    // divide
    CArray even = x[std::slice(0, N/2, 2)];
    CArray odd = x[std::slice(1, N/2, 2)];

    // conquer
    fft(even);
    fft(odd);

    // combine
    for (size_t k = 0; k < N/2; ++k) {
        Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
        x[k] = even[k] + t;
        x[k + N/2] = even[k] - t;
    }
}

// inverse fft (in-place)
void ifft(CArray& x) {
    // conjugate the complex numbers
    x = x.apply(std::conj);

    // forward fft
    fft( x );

    // conjugate the complex numbers again
    x = x.apply(std::conj);

    // scale the numbers
    x /= x.size();
}

int main() {
    const Complex test[] = {1, 2, 3, 4, 4, 3, 2, 1};
    CArray data(test, 8);

    // forward fft
    fft(data);

    std::cout << "Forward FFT:" << std::endl;
    for (int i = 0; i < 8; ++i)
        std::cout << data[i] << std::endl;

    // inverse fft
    ifft(data);

    std::cout << std::endl << "Inverse FFT:" << std::endl;
    for (int i = 0; i < 8; ++i)
        std::cout << data[i] << std::endl;

    return 0;
}