#include <complex>
#include <iostream>
#include "Function.h"
#include <algorithm>

const double pi = acos(-1);

//From numerical recipies
void four11( vector<double>& data, const int isign)
{
	//vector<complex<double> > ret(vals.size());
	//#include <cmath>
	//#include "nr.h"
	//using namespace std;
	//void NR: :four1(Vec_io_dp &data, const int isign)
	//	Replaces data [0. .2*nn-1] by its discrete Fourier transform, if isign is input as 1; or replaces
	//	data [0. .2*nn-1] by nn times its inverse discrete Fourier transform, if isign is input as -1. 
	//	data is a complex array of length nn stored as a real array of length 2*nn. nn MUST be an 
	//	integer power of 2 (this is not checked for!)
	{
		int n , mmax, m, j, istep, i;
		double  wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;

		int nn=data .size( )/2;
		n=nn << 1;
		j=1;
		for (i=1; i<n; i+=2)  
		{
			if  (j>i)  
			{
				swap (data [j-1], data [i-1]);
				swap (data [j], data [i]) ;
			}
			m=nn;
			while (m>=2 && j>m)
			{
				j-=m;
				m>>=1;
			}
			j+=m;
		}
		mmax=2;
		while (n>mmax) 
		{
			istep=mmax << 1;
			theta=isign*(6.28318530717959/mmax);
			wtemp=sin (0.5*theta);
			wpr = -2.0*wtemp*wtemp;
			wpi = sin(theta);
			wr=1.0;
			wi=0.0;
			for (m=1; m<mmax; m+=2) 
			{
				for (i=m; i<n; i+=istep) 
				{
					j=i+mmax;
					tempr=wr*data [j-1] - wi*data [j];
					tempi=wr*data [j] +wi*data [j-1];
					data [j-1]=data [i-1] - tempr;
					data [j] = data [i] - tempi;
					data [i-1] += tempr;
					data [i] += tempi;
				}
				wr=(wtemp=wr)*wpr-wi*wpi+wr;
				wi=wi*wpr+wtemp*wpi+wi;
			}
			mmax=istep;
		}
	}
}


// prints out the values of the samples and the spectrum 
void doFFT(Function& f, vector<double>& vals, double secs, double sams)
{
	double t;
	int i;

	vals.resize(secs*sams*2);

	cout << "# The sampled values are\n";
	cout << "# Time       \t       f(t)\n";

	for (t=0, i=0;t<secs ;t+=(1/sams))
	{
		cout << t << "\t" << (vals[i*2] = f(t)) << "\t" << (vals[i*2+1] = 0) << endl; 
		i++;
	}

	cout << "\n# Here is the spectrum of the sampled data\n";
	four11(vals, 1);
	for (i=0;i<vals.size()/2;i++)
		cout << i*2*pi/secs << "\t" << vals[i*2] << "\t" << vals[i*2+1]<< endl;
}


void main(int argc, char* argv[])
{
	double secs=1, sams;
	Function self;
	FunctionObject(cos);
	Function& f = cosF(self*(6*pi));
	vector<double> vals;

	//Problem 6.13
	do
	{
		cout.precision(14);
		cout.width(14);
		cout << "\n# Welcome to the Fast Fourier Transform Program " << endl;
		cout << "# We will first examine the function :" << endl;
		cout << "# f(t) = cos( 6 * pi * t )" << endl;
		cout << "# How many seconds would you like to sample this function for (0 to quit)? ";
		cin >> secs;
		cout << "# How many samples would you like to take per second? ";
		cin >> sams;

		doFFT(f, vals, secs, sams);
	} while (secs);

	// Problem 6.15
	secs = 1;
	Function& cos3t = cosF(self*3);
	do
	{
		cout << "\n# Please enter the number of seconds to sample cos3t  (or 0 to quit): ";
		cin >> secs;

		doFFT(cos3t, vals, secs, 1);
	} while (secs);
}