#include <iostream>
#include <string>
#include <vector>
#include <complex>
#include <cmath>

const double PI = std::acos(-1);
const std::complex<double> i(0, 1);

using namespace std::complex_literals;

// (Slow) DFT
std::vector<std::complex<double>> dft(const std::vector<std::complex<double>> &v) {
	std::vector<std::complex<double>> result(v.size(), 0);

	const double N = v.size();

	for (int k = 0; k < v.size(); k++) {
		for (int j = 0; j < result.size(); j++) {
			result[k] += v[j] * std::exp(-i * 2. * PI * double(k * j) / N);
		}
	}

	return result;
}

// Cooley-Tukey
std::vector<std::complex<double>> fft1(const std::vector<std::complex<double>> &v) {
	std::vector<std::complex<double>> result(v.size(), 0);
	
	return result;
}

std::vector<double> freqs{ 2, 5, 11, 17, 29 }; // known freqs for testing

void generate(std::vector<std::complex<double>>& v) {
    for (int i = 0; i < v.size(); i++) {
        v[i] = 0.;
        // sum several known sinusoids into v[]
        for(int j = 0; j < freqs.size(); j++)
            v[i] += sin(2 * M_PI * freqs[j] * i / v.size() );
    }
}


// separate even/odd elements to lower/upper halves of array respectively.
// Due to Butterfly combinations, this turns out to be the simplest way 
// to get the job done without clobbering the wrong elements.
void separate (complex<double>* a, int n) {
    complex<double>* b = new complex<double>[n/2];  // get temp heap storage
    for(int i=0; i<n/2; i++)    // copy all odd elements to heap storage
        b[i] = a[i*2+1];
    for(int i=0; i<n/2; i++)    // copy all even elements to lower-half of a[]
        a[i] = a[i*2];
    for(int i=0; i<n/2; i++)    // copy all odd (from heap) to upper-half of a[]
        a[i+n/2] = b[i];
    delete[] b;                 // delete heap storage
}

// N must be a power-of-2, or bad things will happen.
// Currently no check for this condition.
//
// N input samples in X[] are FFT'd and results left in X[].
// Because of Nyquist theorem, N samples means 
// only first N/2 FFT results in X[] are the answer.
// (upper half of X[] is a reflection with no new information).
void fft2 (complex<double>* X, int N) {
    if(N < 2) {
        // bottom of recursion.
        // Do nothing here, because already X[0] = x[0]
    } else {
        separate(X,N);      // all evens to lower half, all odds to upper half
        fft2(X,     N/2);   // recurse even items
        fft2(X+N/2, N/2);   // recurse odd  items
        // combine results of two half recursions
        for(int k=0; k<N/2; k++) {
            complex<double> e = X[k    ];   // even
            complex<double> o = X[k+N/2];   // odd
                         // w is the "twiddle-factor"
            complex<double> w = exp( complex<double>(0,-2.*M_PI*k/N) );
            X[k    ] = e + w * o;
            X[k+N/2] = e - w * o;
        }
    }
}

int main(int argc, char* argv[]) {
	std::vector<std::complex<double>> v(1024, 0);

	generate(v);

	std::vector<std::complex<double>> v_dft = dft(v);

	std::vector<std::complex<double>> v_fft = fft1(v);

	if (v_dft != v_fft) {
		std::cout << "Error: Results differ!\n";
		return 1;
	}
	
	return 0;
}