#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <complex>
#include<vector>
using namespace std;
#ifndef M_PI
#define M_PI           3.14159265358979323846
#endif
typedef struct {
	long N;
	double b, a, beta, alpha;
	double tmin, tmax,h,tl;
	int m;
    int sign;
} fft_param;

complex<double> expo(fft_param *fp,int k) {
    int sign = fp->sign;
    int N=fp->N;
	return exp(sign*2.0 * M_PI * 1i* ((double)(k/N)));
    
}

vector<complex<double>> fourierTrafo(fft_param *fp, vector<complex<double>> x) {
	long N = x.size();
    int sign=fp->sign;
	vector<complex<double>> gerPkt(N / 2);
	vector<complex<double>> ungerPkt(N / 2);
    
	if (N == 1) {
		
		return x;
	} else {
		for (int k = 0; k < N / 2; k++) {
			gerPkt[k] = x[2 * k];
			ungerPkt[k] = x[2 * k + 1];
		}
		
		vector<complex<double>> ftger = fourierTrafo(fp, gerPkt);   
		vector<complex<double>> ftunger = fourierTrafo(fp, ungerPkt); 
		vector<complex<double>> out(N);
        
		for (int k = 0; k < N / 2; k++) {
            if(sign==-1){
                out[k] = (ftger[k] + ftunger[k] * expo(fp,k));
			    out[k + N / 2] = (ftger[k] - ftunger[k] *expo(fp,k));
            }
            else{
                out[k] = (ftger[k] + ftunger[k] * expo(fp,k));
			    out[k + N / 2] = (ftger[k] - ftunger[k] *expo(fp,k));
            }
		}
		
		return out;
    }
}
vector<complex<double>> test(fft_param *fp){
    double h=fp->h;
    double tmin=fp->tmin;
    double tmax=fp->tmax;
    double t=tmin;
    int N=fp->N;
    
    vector<complex<double>> x(N);
    
    for(int k=0;k<N;k++){
        x[k]=sin(tmin+k*h);
    }
    vector<complex<double>> out = fourierTrafo(fp, x);
    fft_param n= *fp;
    n.sign=-1;
    vector<complex<double>> doppelout = fourierTrafo(&n, out);
    FILE *fp2 = fopen("testFFT.csv" ,"w");
    for(int k=0;k<N;k++){
        fprintf(fp2,"%f\t%f\t%f\n",(tmin+k*h),x[k].real(), doppelout[k].real()/(double)N);
    }
    
    return x;
}   
vector<complex<double>> makeS(fft_param *fp){
    double h=fp->h;
    double tmax=fp->tmax;
    double tmin=0;
    double t;
    double tl=fp->tl;
    double alpha=fp->alpha;
    int N=fp->N;
    vector<complex<double>> s(N);
    for(int k=0;k<N;k++){
        t=tmin+k*h;
        s[k]=sin(2*M_PI*alpha*(t-tl));
    }
    return s;
}
vector<complex<double>> makeE(fft_param *fp){
    double h=fp->h;
    double tmax=fp->tmax;
    double tmin=0;
    double t;
    double beta=fp->beta;
    double tl=fp->tl;
    double a=fp->a;
    double b=fp->b;
    int N=fp->N;
    vector<complex<double>> e(N);
    vector<complex<double>> s=makeS(fp);
    for(int k=0;k<N;k++){
        t=tmin+k*h;
       e[k]=b*sin(2*M_PI*beta*t)+ 2*a*(rand()-0.5)+s[k];
       //e[k]=b*sin(2*M_PI*beta*t)+s[k];
    }
    return e;
}
vector<complex<double>> aufgabe3 (fft_param *fp){
    double h=fp->h;
    double tmin=fp->tmin;
    double tmax=fp->tmax;
    double t=tmin;
    int N=fp->N;
    vector<complex<double>> e=makeE(fp);
    vector<complex<double>> s=makeS(fp);
    vector<complex<double>> Fe=fourierTrafo(fp,e);
    fft_param fps=*fp;
    fps.sign=-1;
    vector<complex<double>> Fs=fourierTrafo(&fps,s);
    vector<complex<double>> FFes(N);
    
    for(int k=0;k<N;k++){
        FFes[k]=Fe[k]*conj(Fs[k]);
        
    }
    vector<complex<double>> FFFes=fourierTrafo(fp, FFes);
    return FFFes;
    
}
int main(){
    fft_param fp;
    fp.N= pow(2,13);
    fp.tmin=0;
    fp.tmax = 100;
    fp.h=(fp.tmax-fp.tmin)/fp.N;
    fp.sign=1;
    test(&fp);
    
    
}