#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#define _USE_MATH_DEFINES
#include "math.h"
#include <float.h>
#include <vector>
#include <complex>
using namespace std;

double t1 = -M_PI;
double t2 = M_PI;
double T = t2-t1;// 2 * M_PI;

void fft(vector<complex<double>> & a) {
	int n = (int)a.size();
	if (n == 1)  return;

	vector<complex<double>> a0(n / 2), a1(n / 2);
	for (int i = 0, j = 0; i<n; i += 2, ++j) {
		a0[j] = a[i];
		a1[j] = a[i + 1];
	}
	fft(a0);
	fft(a1);

	double ang = 2 * M_PI / n * 1;//-1
	complex<double> w(1), wn(cos(ang), sin(ang));
	for (int i = 0; i<n / 2; ++i) {
		a[i] = a0[i] + w * a1[i];
		a[i + n / 2] = a0[i] - w * a1[i];

		w *= wn;
	}
}

vector<double> Tgrid(int n) {
	vector<double> t(n);
	for (int i = 0; i < n; ++i)
		t[i] = t1+1.0*i*T / n;
	return t;
}

vector<double> Wgrid(int n) {
	vector<double> w(n);
	for (int i = 0; i < n; ++i)
		w[i] = -1. * M_PI*n / T + 2. * M_PI*i / T;
	return w;
}

double func(double t, double w){
	return cos(w*t)*sin(w*t)* sin(w*t);
	//return sin(21.*t) + sin(23.*t);
}

vector<double> Fgrid(vector<double> t, double w, int n) {
	vector<double> f(n);
	for (int i = 0; i < n; ++i)
		f[i] = func(t[i], w);
	return f;
}

double hann_window(int k, int n) {
	return 1. / 2.*(1 - cos(2 * M_PI*k / (n - 1)));
}

double rect_window(int k, int n) {
	return 1.;
}

void write_to_file(string file, vector<double> w, vector<double> res, int n) {
	ofstream plot;
	plot.open(file);
	plot << "w" << " res" << endl;
	for (int i = 0; i < n; ++i){
		plot << w[i] << " " << res[i] << endl;
	}
	plot.close();
}

vector<double> ft(vector<double> t, vector<double> w, vector<double> f, int n, double win_f(int, int)) {
	vector<double> res(n);

	for (int j = 0; j < n; ++j) {
		double re = 0.;
		double im = 0.;
		for (int k = 0; k < n; ++k) {
			re += f[k] * win_f(k, n) * cos(w[j]*t[k]);
			im += f[k] * win_f(k, n) * sin(w[j]*t[k]);
		}
		res[j] = re*re + im*im;
	}
	return res;
}

	int main(){
		int N = 4000;
		double W = 100.;

		vector<double> f(N);
		vector<double> t(N);
		vector<double> w(N);
		vector<double> res(N);

		t = Tgrid(N);
		w = Wgrid(N);
		f = Fgrid(t, W, N);
		ofstream timefile;
		timefile.open("timve.txt");

		int seconds = time(NULL);

		res = ft(t, w, f, N, &rect_window);
		seconds = time(NULL) - seconds;
		timefile << seconds << endl;
		write_to_file("plot_rect.txt", w, res, N);

		//res = ft(t, w, f, N, &hann_window);
		//write_to_file("plot_hann.txt", w, res, N);

		//N = 1024;
		//t = Tgrid(N);
		//w = Wgrid(N);
		//f = Fgrid(t, W, N);
		vector<complex<double>> fc(N);
		for (int i = 0; i < N; ++i)
			fc[i] = complex<double> (f[i], 0);
		vector<double> fcm(N);
		seconds = time(NULL);
		fft(fc);
		seconds = time(NULL) - seconds;
		timefile << seconds << endl;
		timefile.close();
		for (int i = 0; i < N; ++i)
			fcm[i] = fc[i].real()*fc[i].real() + fc[i].imag()*fc[i].imag();
		write_to_file("plot_ftt.txt", w, fcm, N);

		return 0;
	}