/***************************************************************
 *                                                             *
 * Take the FFT of a dataset and export as a power spectrum    *
 *                                                             *
 * to compile:                                                 *
 *    g++ pspectra_Vx.cpp -lhdf5_hl -lhdf5 -Wall -pedantic \   *
 *        -lfftw3 readinput.cpp *                              *
 *							       * 
 *  v5.2  Changed formula for binning        		       *
 *  v5.1  Added labview v3 fileformat support		       *
 *  v5.0  Attempt to add welch windowign                       *
 *        added flat top window                                * 
 *  v4.10 Removed longdouble				       *
 *  v4.9 Introduced longdouble                                 *
 *  v4.8 Added max checking and output                         *
 *  v4.8 Added max checking and output                         *
 *  v4.7 Clariefied help section a bit                         *
 *  v4.6 Fix bug with freq calculation, add log binnig         *
 *  v4.5 Attempt at adding hdf5 output, added a couple of      *
 *       window functions and set Hamming default              *
 *       Run through valgrind                                  *
 *  v4.4 Added sine window      			       *
 *  V4.3 Modified code to accept a double time,added no window *
 *  V4.2 Added colour and shift of data,default no normalize   *
 *  V4.1 Changed default to not normalize or bin               *
 *  V4 Checked working with variance = sum(S(f)*df)            *
 *     matlab var(test_raw_windowed),sum(test_psd(2:end,:)*0.1)*
 *  V3 Small corrections and choice of windowing function      *
 *  V2 History                                                 *
 *  V1 Ancient History                                         *
 *                                                             *
 *                                                             *
 *                                                             *
 *                                                             *
 *                                                             *
 ***************************************************************/

#include <iostream>
#include <cmath>
#include <complex>
#include "fftw3.h"
#include <hdf5.h>
#include <H5LT.h>
#include "readinput.h"

//Add colour output to make things more readable
//Check if in cygwin, and alter coloura approriately.
#if defined(__GNUC__) && defined(__CYGWIN__)
#define redc "\033[22;31m" 
#define normc "\033[01;37m" 
#define bluec "\033[22;34m"
#define yellowc "\033[01;33m"
#else
#define redc "\033[22;31m" 
#define normc "\033[01;37m" 
#define bluec "\033[22;34m"
#define yellowc "\033[01;33m"

#endif
//Max voltage before warning
#define maxV 9.98

//Allocate a 2d array compatible with hdf5 library
using namespace std;
    template < typename T >
T **Allocate2DArray( int nRows, int nCols)
{
    T **ppi;
    T *pool;
    T *curPtr;
    //(step 1) allocate memory for array of elements of column

    ppi = new T*[nRows];

    //(step 2) allocate memory for array of elements of each row
    pool = new T [nRows * nCols];

    // Now point the pointers in the right place
    curPtr = pool;
    for( int i = 0; i < nRows; i++)
    {
	*(ppi + i) = curPtr;
	curPtr += nCols;
    }
    return ppi;
}

//Deallocate the 2d array
template < typename T >
void Free2DArray(T** Array)
{
    delete [] *Array;
    delete [] Array;
}

//Initialize function definitions
long int power(int x, int n); 
long count_lines (string filename);
double average(double * array,int i,int f);
inline double hann(int n,int N);
void version(void);
void help(string filename);
inline double trianglew(int n,int N);
inline double sinew(int n, int N);
inline double bartletthann(int n,int N);
inline double fourtBlackHarr(int n,int N);
inline double fourtBlackHarr90db(int n,int N);
inline double Hamming(int n,int N);
inline double flattop(int n, int N);
double window(int n,int N, int argc, char* argv[]);
inline void coutr(string red);

int main(int argc, char* argv[]) {
    cout << endl<<endl;
    // Check input for version call
    if (readinputb(argc,argv,"V")){
	version();
    }
    bool showHelp=false;
    //Check input for help call
    if (readinputb(argc,argv,"H") ||  readinputb(argc,argv,"h") )
	showHelp=true;
    //Initialize all variables
    long warn=0;
    int bins=0;
    bool Force=false;
    int nan=0;
    string filename,outfilename, infofilename, xfilename, yfilename, zfilename;
    bool shift = !readinputb(argc,argv,"S");
    bool welch = true;
    int segments = 1;
    int overlap= 0;
    if(welch){
	segments = readinputi(argc,argv,"s",1);
	overlap=readinputi(argc,argv,"l",0);
    }
    double shiftx=0,shifty=0;
    //Make filename variables and set up program depending on number of inputs
    if(2 == argc ){
	filename=argv[1];
	infofilename=filename + ".txt";
	if ( !file_exists(infofilename))
	    help(argv[0]);
	outfilename=filename;
	xfilename=filename + "_x.hws";
	yfilename=filename + "_y.hws";
	zfilename=filename + "_z.hws";
	bins=0;
    }
    else if (5>argc)
	showHelp=true;
    else{
	bins=readinputi(argc,argv,"b",0);
	Force = readinputb(argc,argv,"F");
	filename=readinputs(argc,argv,"i","infile");
	outfilename=readinputs(argc,argv,"o","outfile");
	infofilename = filename + ".txt";
	xfilename =filename +"_x.hws";
	yfilename =filename +"_y.hws";
	zfilename =filename +"_z.hws";
    }
    string outifilename=outfilename+"_psd.txt";
    if(showHelp)
	help(argv[0]);
    //Set variables from stdin
    bool normalize = !readinputb(argc,argv,"n");
    bool nocomment = readinputb(argc,argv,"D");
    int find_max  = readinputb(argc,argv,"M");
    // Check all files exist
    if( !file_exists(infofilename)){
	cout <<"Please specify a valid input information file"<<endl<<
	    "e.g. E01.txt " <<endl;
	exit(EXIT_FAILURE);
    }
    if( !file_exists(xfilename) ||  !file_exists(xfilename)  ||  !file_exists(xfilename) ){
	cout <<"Please verify valid data files"<<endl<<
	    "e.g. E01_x.txt, E01_y.txt and E01_z.txt " <<endl;
	exit(EXIT_FAILURE);
    }
    FILE * InfoFile;
    InfoFile = fopen(infofilename.c_str(),"r");
    cout << "Reading infofile:"<<infofilename<<endl;
    // Create a variable to hold experimental info
    char test[200];
    //Read last line (bit silly but works)
    for(int i =0;i<99;i++){
	fgets(test,99,InfoFile);
    }
    // initialize variables used in data analysis
    int samples,version=-1;
    double  tempavg,tempi,tempf,bandwidth,time,particlesize,drivefreq,drivedist,power,samples_d;
    char version_ascii=0;
    //Determine length of line (1st or second version
    for(int j=0;version<1;j++){
	version_ascii=test[strlen(test)-1-j];
	if(version_ascii<'1' || version_ascii>'9')
	    version_ascii='0';
	version=version_ascii-48;
	if (j>5 || version>3){
	    cout << "Wrong specifications file, please correct"<<endl;
	    cout << "Final line should be:"<<endl;
	    cout << "\"Sampling Freq\"   \"Time\"  \"Avg Temp\" \"I temp\" \"F temp\" \"Version\""<<endl;
	    exit(EXIT_FAILURE);
	}
    }
    //Read input file, depending on version
    if(version==3)
	sscanf(test,"%le %le %le %le %le %le %le %le %le %d",&samples_d,&time,&tempavg,&tempi,&tempf,&particlesize,&drivefreq,&drivedist,&power,&version);
    else if(version==2)
	sscanf(test,"%le %le %le %le %le %le %le %le %d",&samples_d,&time,&tempavg,&tempi,&tempf,&particlesize,&drivefreq,&drivedist,&version);
    else if (version==1)
	sscanf(test,"%le %le %le %le %le %d",&samples_d,&time,&tempavg,&tempi,&tempf,&version);
    samples=(int)samples_d;
    //read from the test variable and read info
    fclose(InfoFile);
    //If we are splitting up the datafiles, find size of segments
    long orglines=(long)(samples*time);
    long lines=orglines/segments;
    long freq=samples;
    long overlaplines=orglines/segments*overlap/100;
    long discardedpoints=orglines+overlaplines*segments-lines*segments;
    while(discardedpoints>lines){
	segments++;
	discardedpoints=orglines+overlaplines*segments-lines*segments;
	cout << "Added segment, due to overlap"<<endl;
    }

    //long lines=count_lines(filename);
    cout << "\n\nFFT of Optical tweezer data\n"<<
	"Note that the number of bins is taken as an approximate\n"<<
	"Sampling frequency : " << freq << endl<<
	"In Filename:"<<filename<<"\nOutfilename :"<<outifilename<<endl<<
	"Sampling rate:"<<samples<<" samples/s"<<endl<<
	"Sampling time:"<<time<<" s"<<endl<<
	"File version:"<<version<<endl<<endl;
    if(!normalize) 
	cout << bluec<<"You have choosen not to normaize the data"<<normc<<endl<<endl;
    double * freq_list= new double[lines/2];
    cout << "Generating frequency data";
    for(int i=0;i<lines/2;i++){
	freq_list[i]=samples*(double)i/lines;
    }
    cout <<", done"<<endl<<endl;

    hid_t       file_id; 
    herr_t      status;
    bandwidth=(double)lines*(double)lines/(double)time/2.0;
    /******************************************
     *                                        *
     *   Import Z Data and take FFT           *
     *   Conjugate and save to standard array *
     *                                        *
     ******************************************/
    double * zdata = new double[lines/2];
    double * raw_z_data = new double[orglines];
    cout << "Loading Zdata into memory: " << zfilename;
    // Open HDF File
    file_id = H5Fopen (zfilename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); 
    status = H5LTread_dataset_double(file_id,"/wfm_group0/axes/axis1/data_vector/data",raw_z_data);
    status = H5Fclose (file_id);
    std::complex<double>*  inz=new std::complex<double> [lines];
    std::complex<double>* outz=new std::complex<double> [lines];
    cout << ", Starting FFT\n";
    //Segment and analyze
    for(int segs=0;segs<segments;segs++){
	for(int i=0;i<lines;i++){
	    inz[i].real()=raw_z_data[i]*window(i,lines,argc,argv);
	    if(raw_z_data[i]>9.98)
		warn++;
	    inz[i].imag()=0.0;
	}
	//Perform FFT
	fftw_plan pz = fftw_plan_dft_1d(lines,reinterpret_cast<fftw_complex*>(inz), 
		reinterpret_cast<fftw_complex*>(outz), 
		FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(pz); 
	for(int i=0;i<lines/2;i++)
	    zdata[i]=real(outz[i]*conj(outz[i]))/bandwidth;
	fftw_destroy_plan(pz);
    }
    delete [] inz;
    delete [] outz;
    cout << "FFT'd and conjugated Z data"<<endl<<endl;
    /******************************************
     *                                        *
     *   Import X Data and take FFT           *
     *   Conjugate and save to standard array *
     *                                        *
     ******************************************/
    std::complex<double>*  inx=new std::complex<double> [lines];
    std::complex<double>* outx=new std::complex<double> [lines];
    double * xdata = new double[lines/2];
    for(int i=0;i<lines/2;i++)
	xdata[i]=0;
    double * tempdata = new double[lines];
    double * rawxdata = new double[orglines];
    cout << "Loading Xdata into memory: " << xfilename;
    file_id = H5Fopen (xfilename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); 
    status = H5LTread_dataset_double(file_id,"/wfm_group0/axes/axis1/data_vector/data",rawxdata);
    status = H5Fclose (file_id);
    cout << ", Starting FFT\n";
    for(int segs=0;segs<segments;segs++){
	int segstart=segs*lines-segs*overlaplines;
	for(int p=0;p<lines;p++)
	    tempdata[p]=rawxdata[p+segstart];
	if(shift)
	    shiftx=average(tempdata,0,lines-1);
	if(normalize)
	    for(int i=0;i<lines;i++){
		if(raw_z_data[i]==0){
		    cout << "Please rerun without normalization, divide by zero"<<endl;
		    exit(EXIT_FAILURE);
		}
		inx[i].real()=(tempdata[i]-shiftx)/raw_z_data[i+segstart]*window(i,lines,argc,argv);
		inx[i].imag()=0.0;
	    }
	else 
	    for(int i=0;i<lines;i++){
		inx[i].real()=(tempdata[i]-shiftx)*window(i,lines,argc,argv);
		inx[i].imag()=0.0;
	    }

	fftw_plan px = fftw_plan_dft_1d(lines,reinterpret_cast<fftw_complex*>(inx), 
		reinterpret_cast<fftw_complex*>(outx), 
		FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(px); 
	for(int i=0;i<lines/2;i++){
	    xdata[i]+=real(outx[i]*conj(outx[i]))/bandwidth;
	}
	fftw_destroy_plan(px);
    }
    delete [] rawxdata;
    delete [] inx;
    delete [] outx;
    cout << "FFT'd and conjugated X data"<<endl<<endl;
    /******************************************
     *                                        *
     *   Import Y Data and take FFT           *
     *   Conjugate and save to standard array *
     *                                        *
     ******************************************/
    std::complex<double>*  iny=new std::complex<double> [lines];
    std::complex<double>* outy=new std::complex<double> [lines];
    double * ydata = new double[lines/2];
    for(int i=0;i<lines/2;i++)
	ydata[i]=0;
    double * rawydata = new double[orglines];
    cout << "Loading Ydata into memory: " << yfilename;
    file_id = H5Fopen (yfilename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); 
    status = H5LTread_dataset_double(file_id,"/wfm_group0/axes/axis1/data_vector/data",rawydata);
    status = H5Fclose (file_id);
    tempdata[0]=5;
    cout << "Starting FFT\n";
    for(int segs=0;segs<segments;segs++){
	int segstart=segs*lines-segs*overlaplines;
	for(int p=0;p<lines;p++)
	    tempdata[p]=rawydata[p+segstart];
	if(shift)
	    shifty = average(tempdata,0,lines-1);
	if(normalize)
	    for(int i=0;i<lines;i++){
		if(raw_z_data[i]==0){
		    cout << "Please rerun without normalization, divide by zero"<<endl;
		    exit(EXIT_FAILURE);
		}
		iny[i].real()=(tempdata[i]-shifty)/raw_z_data[i+segstart]*window(i,lines,argc,argv);
		iny[i].imag()=0.0;
	    }
	else
	    for(int i=0;i<lines;i++){
		iny[i].real()=(tempdata[i]-shifty)*window(i,lines,argc,argv);
		iny[i].imag()=0.0;
	    }

	fftw_plan py = fftw_plan_dft_1d(lines,reinterpret_cast<fftw_complex*>(iny), 
		reinterpret_cast<fftw_complex*>(outy), 
		FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(py); 
	for(int i=0;i<lines/2;i++)
	    ydata[i]+=real(outy[i]*conj(outy[i]))/bandwidth;
	fftw_destroy_plan(py);
    }
    delete [] tempdata;
    delete [] rawydata;
    delete [] iny;
    delete [] outy;
    cout << "FFT'd and conjugated Y data"<<endl<<endl;
    for(int i =0;i<lines/2;i++){
	ydata[i]=ydata[i]/segments;
	xdata[i]=xdata[i]/segments;
    }

    FILE * outfile;
    outfile = fopen(outifilename.c_str(),"w");
    double xmax_freq=0,xmax_value=0;
    double ymax_freq=0,ymax_value=0;

    /*****************************************
     * Find Max
     * **************************************/

    if(find_max){
	for(int i=0;i<lines/2;i++){
	    if(freq_list[i]>12 && freq_list[i]<42){
		if(xdata[i]>xmax_value){
		    xmax_freq=freq_list[i];
		    xmax_value=xdata[i];
		}
	    }
	}
	for(int i=0;i<lines/2;i++){
	    if(freq_list[i]>12 && freq_list[i]<42){
		if(ydata[i]>ymax_value){
		    ymax_freq=freq_list[i];
		    ymax_value=ydata[i];
		}
	    }
	}
	printf("XMax at %g Hz with %g\n",xmax_freq,xmax_value);
	printf("YMax at %g Hz with %g\n",ymax_freq,ymax_value);
    }

    /*****************************************
     * 
     * Write comments to file 
     *
     * **************************************/

    if(!nocomment){
	if((warn>0) && normalize)
	    fprintf(outfile,"%%WARNING NORMALIZATION MAY NOT BE CORRECT, VOLTAGE Exceeded %gV %ld times (%ld%%)\n",maxV,warn,warn*100/lines);
	fprintf(outfile,"%%%s %s\n","Source file name: ",__FILE__);
	fprintf(outfile,"%%%s %s %s %s %s","Compiled at: ",__DATE__," ",__TIME__," with GNU g++ ");
	fprintf(outfile,"%%%d.%d.%d\n", __GNUC__ ,__GNUC_MINOR__ ,__GNUC_PATCHLEVEL__);
	fprintf(outfile,"%%%s %s\n","FFTW Version: ",fftw_version);
	fprintf(outfile,"%%%s %d.%d.%d.%s\n","HDF5 Version: ",H5_VERS_MAJOR,H5_VERS_MINOR,H5_VERS_RELEASE,H5_VERS_SUBRELEASE);
	if (Force)
	    fprintf(outfile,"%%Input version not correct, user override\n");
	if(normalize)
	    fprintf(outfile,"%%Data was normalized using Z\n");
	if(find_max)
	    fprintf(outfile,"%%Max at %g Hz with %g\n",xmax_freq,xmax_value);
	switch (readinputc(argc,argv,"W",'Q')){
	    case 'H':
		fprintf(outfile,"%%Windowing Function: Hann\n");
		break;
	    case 'T':
		fprintf(outfile,"%%Windowing Function: Triangle\n");
		break;
	    case '9':
		fprintf(outfile,"%%Windowing Function: 4term BlackmannHarris(90db)\n");
		break;
	    case '4':
		fprintf(outfile,"%%Windowing Function: 4term BlackmannHarris(72db)\n");
		break;
	    case 'B':
		fprintf(outfile,"%%Windowing Function: Barlett-Hann\n");
		break;
	    case '1':
		fprintf(outfile,"%%Windowing Function: None/rectangular\n");
		break;
	    case 'S':
		fprintf(outfile,"%%Windowing Function: Sine\n");
		break;
	    case 'A':
		fprintf(outfile,"%%Windowing Function: Hamming\n");
		break;
	    case 'F':
		fprintf(outfile,"%%Windowing Function: Flat top\n");
		break;
	    default:
		fprintf(outfile,"%%Windowing Function: Hamming(Default)\n");
	}
	if(nan >0 )
	    fprintf(outfile,"\n\nWARNING ENCOUNTERED NaN's, %d times\n\n",nan);
    }

    /**************************************************
     *
     * Put data in correct array for writing to file.
     *
     * ***********************************************/
    double * fxdata = new double[lines/2];
    double * fydata = new double[lines/2];
    double * fzdata = new double[lines/2];
    double * ffdata = new double[lines/2];
    int numlines=0;
    for(int i=0;i<lines/2;i++)
	fxdata[i]=fydata[i]=fzdata[i]=ffdata[i]=0.0;
    if(bins ==0){
	if(!nocomment)
	    fprintf(outfile,"%%Binning was not performed\n\n");
	numlines=lines/2;
	for(int i=0;i<lines/2;i++){
	    ffdata[i]=freq_list[i];
	    fxdata[i]=xdata[i];
	    fydata[i]=ydata[i];
	    fzdata[i]=zdata[i];
	}
    }
    else {
	double min=freq_list[1];
	double max=lines/2;
	double delta=log(max/min)/bins;
	int * count=new int[lines/2];
	for(int i=0;i<lines/2;i++){
	    count[i]=0;
	}
	int bin=0;
	ffdata[0]=freq_list[0];
	fxdata[0]=xdata[0];
	fydata[0]=ydata[0];
	fzdata[0]=zdata[0];
	count[0]=1;
	int i;
	if (drivefreq <= 0 )
	    drivefreq=20;
	    
	for(i=1;ffdata[i-1]<drivefreq+4;i++){
	    ffdata[i]=freq_list[i];
	    fxdata[i]=xdata[i];
	    fydata[i]=ydata[i];
	    fzdata[i]=zdata[i];
	    count[i]++;
	}
	for(;i<lines/2;i++){
	    bin=(int)(log((double)i/min)/delta);
	    ffdata[bin]+=freq_list[i];
	    fxdata[bin]+=xdata[i];
	    fydata[bin]+=ydata[i];
	    fzdata[bin]+=zdata[i];
	    count[bin]++;
	}
	int j=0;
	for(int i=0;i<lines/2;i++){
	    ffdata[j]=ffdata[i]/count[i];
	    if(isnan(ffdata[j]))
		j--;
	    else{
		fxdata[j]=fxdata[i]/count[i];
		fydata[j]=fydata[i]/count[i];
		fzdata[j]=fzdata[i]/count[i];
	    }
	    j++;
	}
	numlines=j;
	//if(!nocomment)
	//    fprintf(outfile,"%%Binning was performed with %d specified restulting in a binning coeffient of %g\n\n", bins, bincoef);
    }
    delete [] xdata;
    delete [] ydata;
    delete [] zdata;
    delete [] raw_z_data;
    delete [] freq_list;

    /***************************************************
     * 
     * Write to file 
     *
     * ************************************************/
    cout << "Saving data..."<<endl;
    if(readinputb(argc,argv,"t")){
	for(int i=0;i<numlines;i++){
	    if(ffdata[i]-0.00001>ffdata[i-1])
		fprintf(outfile,"%20.9e,%20.9e,%20.9e,%20.9e\n", ffdata[i],fxdata[i],fydata[i],fzdata[i]);
	}
    }
    else{	
	hid_t   file_oid,dataspaceo,dataset_oid;
	hsize_t dimso[2];
	herr_t  statuso;
	double ** datx=Allocate2DArray<double>(4,numlines);
	for(int i=0;i<numlines;i++){
	    datx[0][i]=ffdata[i];
	    datx[1][i]=fxdata[i];
	    datx[2][i]=fydata[i];
	    datx[3][i]=fzdata[i];//freq_list[i];
	}
	string h5file=outfilename+"_psd.hdf5";
	file_oid=H5Fcreate(h5file.c_str(),H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
	dimso[0]=4;
	dimso[1]=numlines;
	dataspaceo=H5Screate_simple(2,dimso,NULL);
	dataset_oid=H5Dcreate(file_oid,"dset",H5T_IEEE_F64BE,dataspaceo,H5P_DEFAULT);
	statuso=H5Dwrite(dataset_oid,H5T_NATIVE_DOUBLE,H5S_ALL,H5S_ALL,H5P_DEFAULT,&datx[0][0]);
	statuso=H5Dclose(dataset_oid);
	statuso=H5Sclose(dataspaceo);
	statuso=H5Fclose(file_oid);
	Free2DArray(datx);
    }
    /********************************************
      Delete stuff
     ****************************************/

    delete [] ffdata;
    delete [] fxdata;
    delete [] fydata;
    delete [] fzdata;
    fclose(outfile);
    if(warn && normalize)
	cout << redc << "\n\nWARNING NORMALIZATION MAY NOT BE CORRECT, VOLTAGE EXceeded "<<maxV<< "V "<<warn<<" times ("<<warn*100/lines << "%%)\n\n"<<normc;
    if(nan >0 )
	cout << redc << "\n\nWARNING ENCOUNTERED NaN's, "<< nan<< " times\n\n"<<normc;
    return 0;
}

//Custom power function which takes two integers
long int power(int x, int n) {
    long int sol=x;
    for(int i=1;i<n;i++){
	sol*=x;
    }
    return sol;
}

//Find how many files a file contains
long count_lines (string filename)
{
    FILE *wcdata;
    wcdata = fopen(filename.c_str(),"r");
    int c;
    unsigned long num_lines = 0;

    while ((c = getc (wcdata)) != EOF)
	if (c == '\n')
	    num_lines++;
    fclose(wcdata);
    return num_lines;
}

//Average an array
double average(double * array,int i,int f){
    double sum=0;
    for(int x=i;x<=f;x++){
	sum+=array[x];
    }
    sum/=(f-i+1);
    return sum;
}

//Return windowed function
double window(int n,int N, int argc, char* argv[]){
    /* window with the appropriate windowing function */
    switch (readinputc(argc,argv,"W",'Q')){ //No default return
	case 'H':
	    return hann(n,N);
	    break;
	case 'T':
	    return trianglew(n,N);
	    break;
	case 'B':
	    return bartletthann(n,N);
	    break;
	case 'R':
	    return 1;
	    break;
	case 'S':
	    return sinew(n,N);
	    break;
	case '4':
	    return fourtBlackHarr(n,N);
	    break;
	case '9':
	    return fourtBlackHarr90db(n,N);
	    break;
	case 'A':
	    return Hamming(n,N);
	    break;
	case 'F':
	    return flattop(n,N);
	    break;
	default:
	    return Hamming(n,N);
    }
}

/**********************************************************
 *                                                        *
 * Windowing functions                                    *
 *                                                        *
 *********************************************************/

inline double Hamming(int n,int N){
    return 0.54 -0.46*cos(2*M_PI*n/N);
}
inline double fourtBlackHarr(int n,int N){
    return 0.40217 -0.49703*cos(2*M_PI*n/N)+0.09392*cos(4*M_PI*n/N)-0.00183*cos(6*M_PI*n/N);
}

inline double fourtBlackHarr90db(int n,int N){
    return 0.35875 -0.48829*cos(2*M_PI*n/N)+0.14128*cos(4*M_PI*n/N)-0.01168*cos(6*M_PI*n/N);
}
inline double hann(int n,int N){
    /* calculates the hann coeffients 
     * w(n)=0.5*(1-cos(2*pi*n/N)) */
    return 0.5*(1-cos(2*M_PI*(double)n/((double)N)));

}

inline double sinew(int n, int N){
    /* Calculate the sine window 
     * w(n) =  sin(pi*n/(N-1)) */
    return sin(M_PI*n/(N-1));
}
inline double trianglew(int n,int N){
    /* Calculate the triangle window coeffiecntes
     * w(n)=2/N*(N/2-abs(n-(N-1)/2)) */
    return 2.0/N*(N/2.0-fabs(n-(N-1.0)/2.0));

}

inline double bartletthann(int n,int N){
    /* calculate the barlett-hann window 
     * w(n)= 0.62--0.48*abs(n/(N-1)-0.5)-0.38*cos(2*pi*n/(N-1)) */
    return 0.62-0.48*fabs(n/(N-1)-0.5)-0.38*cos(2*M_PI*n/(N-1));
}

inline double flattop(int n, int N){
    /* Calculate the flattop window
     * w(n) = a0 - a1 cos(2 pi n/(N-1)) +a2 cos(4 pi n/(N-1))-a3cos(6 pi n)/(N-1))+a4 cos(8 pi n/(N-1)) */
    return 1 -1.93*cos(2*M_PI*n/(N-1))+1.29*cos(4*M_PI*n/(N-1))+0.388*cos(6*M_PI*n/(N-1))+0.032*cos(8*M_PI*n/(N-1));
}


//Return version information
void version(){
    cout << 
	"Source file name: "<< __FILE__ <<endl<<
	"Compiled at: "<< __DATE__ <<" "<< __TIME__<< " with GNU g++ " << __GNUC__ << "." << __GNUC_MINOR__ << "." << __GNUC_PATCHLEVEL__ << endl << endl<<
	"FFTW Version: "<<fftw_version<<endl<<
	"HDF5 Version: "<<H5_VERS_MAJOR<<"."<<H5_VERS_MINOR<<"."<<H5_VERS_RELEASE<<"."H5_VERS_SUBRELEASE<<endl<<endl;
    exit(EXIT_SUCCESS);
}

//Return help information
void help(string filename){
    cout << "Minimum call: "<< filename<<  "  -i infile -o outfile.txt" << endl << endl <<
	"General:\n"<<
	"    -i    input file name, without extension or channel\n"<<
	"    -o    output file name, where the psd will be placed it will be appended by _psd\n"<<
	"    -t    do not write hdf5 file, write to txt file\n"<<
	"    -F    Force using data even if version can't be confirmed\n"<<
	"    -H    Show help\n"<<
	"    -V    Version info"<<
	endl<<
	"Options Raw data:\n" <<
	"    -n    do not normalize the output with z voltage\n"<<
	"    -S    do not center around 0\n"<<
	"    -W    Window functions:"<<
	"          H    hann window\n"<<
	"          B    bartlett-Hann window\n"<<
	"          T    triangle window \n"<<
	"          1    no window function\n"<<
	"          S    sine window\n"<<
	"          4    4 term blackmannHarris (72db)\n"<<
	"          9    4 term blackmannHarris (90db)\n"<<
	"          A    Hamming\n"<<
	"	   F    Flat-top\n"<<
	"   -S	#  Segmented windowing\n"
	"   -l	#  with # segments and #%% overlap\n"<<
	endl<<
	"Options PSD Data\n"<<
	"    -b    Number of bins to use, 0 mean no binning, this is used as an appoximate\n"<<
	"    -D    make output datafile comment free\n"<<
	"    -M    report the maximum frequency and value \n"<<
	endl;
    exit(EXIT_SUCCESS);
}


/* ************************************************
 * String colors :
 * 033[22;30m - black
 \033[22;31m - red
 \033[22;32m - green
 \033[22;33m - brown
 \033[22;34m - blue
 \033[22;35m - magenta
 \033[22;36m - cyan
 \033[22;37m - gray
 \033[01;30m - dark gray
 \033[01;31m - light red
 \033[01;32m - light green
 \033[01;33m - yellow
 \033[01;34m - light blue
 \033[01;35m - light magenta
 \033[01;36m - light cyan
 \033[01;37m - white 
 *************************************************/

inline void coutr(string red){
    cout << redc <<red <<normc;
}


