//
//  Fourier transform
//

#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>

#include <time.h>
#include <sys/types.h>
#include <sys/timeb.h>

#define MAX 200000

void elapsed(long *);
void hanning(const long &last);
void meanremove(const double &ave,const long &last);
void read_data(void);


void line_test(void);

typedef struct 
{
    double r,i;
	
}doublecomplex;

 time_t ltime;
 time_t ltimen;

double a[MAX];
double arg,ave;
double df,dt;
double fs,freq;

double sr;
double t,t1,tmax,ts;
double z;

double rms;

doublecomplex f;

double const  pi=   atan2(0.,-1);
double const tpi=2.*pi;

long i,iflag,ihw,k, last;
long n,nm1;

long irem,imr;


int kflag;

double tstart;
double tnow;


double peak_f=0.;

double peak_m=-1.;

FILE *pFile[10];
char filename[10][FILENAME_MAX];

long numBytes=300;

char string[300];



int mflag=0;


int main( int argc, char *argv[] )
{

	if(argc > 1 )
	{

		printf("\n\n  argc = %d ", argc );
		printf("\n argv[0] %s ", argv[0] );
		printf("\n argv[1] %s ", argv[1] );


		strcpy (filename[0], argv[1] );

		mflag=1;

	}


    printf("\n fourier.cpp,  ver 5.0, May 17, 2010\n");
	printf("\n by Tom Irvine ");
	printf("\n Email:  tomirvine@aol.com \n");

	printf("\n This program calculates the Fourier transform of a time history. \n");


	printf("\n It is suitable for small data sets < 20000 points. ");

	printf("\n Otherwise, an FFT should be calculated. \n");


    printf( "\n The input file must contain two columns: \n");
    printf( " time(sec)  amplitude(units)    \n");


	if(mflag==0)
	{
		printf( "\n Input filename \n");

		scanf("%s",filename[0]);

	}


    pFile[0] = fopen(filename[0], "rb");

	while(pFile[0] == NULL )
	{

		printf("\n Failed to open file: %s \n", filename[0]);

		printf("\n Please enter the input filename: \n");

		scanf("%s",filename[0]);

		pFile[0] = fopen(filename[0], "rb");

	}

	printf(" File: %s opened. \n", filename[0]);


	for(i=0; i< MAX; i++)
	{
		a[i]=0.;
	}

	read_data();


	printf("\n n= %ld ",n);
    
	printf("\n tmax= %10.4g sec \n",tmax);
	

//	printf("\n sr= %10.4g samples/sec ",fs);

//	printf("\n dt= %10.4g sec \n",1./fs);


	printf("\n df= %10.4g Hz \n",df);

    printf("\n          mean = %10.4g ",ave);
    printf("\n overall level = %10.4g rms \n",rms);

	

	printf("\n Mean removal? 1=yes 2=no \n");

	scanf("%d",&imr);



	printf("\n Window? 1=rectangular 2=Hanning \n");

	scanf("%d",&ihw);

	

	if(imr==1){ meanremove(ave,last);}

	if(ihw==2){ hanning(last);}

    
    
    strcpy(filename[1],"fourier.out");
	pFile[1]=fopen(filename[1],"w");

	strcpy(filename[2],"ft_half.out");
	pFile[2]=fopen(filename[2],"w");

	strcpy(filename[3],"ft_full.out");
	pFile[3]=fopen(filename[3],"w");

	strcpy(filename[4],"psd.out");
	pFile[4]=fopen(filename[4],"w");

	strcpy(filename[6],"FRMS.out");

	pFile[6]=fopen(filename[6],"w");


	strcpy(filename[7],"ft_mp.out");

	pFile[7]=fopen(filename[7],"w");


	time( &ltime );
	tstart=double(ltime);


    if(n> 200)printf( "\n\n Current time and date: \n   %s", ctime( &ltime ) );

    double ms=0.;  


    iflag =0;

    for( k=0; k < n; k++)
	{
          f.r=0.;
          f.i=0.;


		  if(  (iflag == 0 && k > 0)  )
		  {
			 elapsed(&k);
		  }
          
		  for( i=0; i < n; i++)
		  {
             arg=tpi*k*( double(i)/double(n) );

			 f.r += a[i]*cos(arg);
             f.i += a[i]*sin(arg);

		  }

          freq=k*df;
          
		  f.r/=double(n);
		  f.i/=-double(n);

		  z=sqrt(pow(f.r,2)+pow(f.i,2));

          fprintf(pFile[1],"%10.3e \t %10.3e \t %10.3e\n",freq,f.r,f.i);
		  fprintf(pFile[2],"%10.3e \t %10.3e \n",freq,z);

	      double phase=180*atan2(f.i,f.r)/pi;
		  

		  if(k<(n/2))  
		  {
				 if(k == 0)
				 {
					fprintf( pFile[3],"%15.6e \t %15.6e \n", freq, z);

					fprintf( pFile[7],"%15.6e \t %15.6e \t %15.6e\n", freq, z,phase);


					if(z>peak_m)
					{
						peak_m=z;
						peak_f=freq;

					}
				 }
				 else
				 {
					fprintf( pFile[3],"%15.6e \t %15.6e \n", freq, 2.*z);
					fprintf( pFile[4],"%15.6e \t %15.6e \n", freq, 0.5*pow((2.*z),2.)/df);
				 
					fprintf( pFile[7],"%15.6e \t %15.6e \t %15.6e\n", freq, 2.*z,phase);


				    z*=2.;


					if(z>peak_m)
					{

						peak_m=z;
						peak_f=freq;

					}

				 }

				
				 z*=0.707;

				 ms+=pow(z,2);

				 fprintf( pFile[6],"%15.6e \t %15.6e \n", freq, sqrt(ms));

		  }
    }
    fclose( pFile[1]);
	fclose( pFile[2]);
	fclose( pFile[3]);
	fclose( pFile[4]);
      
	fclose( pFile[6]);
	fclose( pFile[7]);


	printf("\n The peak magnitude occurs at: %9.5g Hz\n",peak_f);

    printf("\n\n **** Output Files **** ");


	printf("\n\n Complex output:  fourier.out ");
	printf("\n The format is: ");
	printf("  freq (Hz)  real amp(units)  imaginary amp(units)");


	printf("\n\n Double-sided, half-amplitude magnitude:  ft_half.out ");
	printf("\n The format is: ");
	printf("  freq (Hz)  magnitude(units) \n");


	printf("\n\n Single-sided, full-amplitude magnitude:  ft_full.out ");
	printf("\n The format is: ");
	printf("  freq (Hz)  magnitude(units) \n");


	printf("\n\n Single-sided, full-amplitude magnitude & phase:  ft_mp.out ");
	printf("\n The format is: ");
	printf("  freq (Hz)  magnitude(units)  phase(deg)");


    printf("\n\n Single-sided, power spectral density:  psd.out ");
	printf("\n The format is: ");
	printf("  freq (Hz)  magnitude(units^2/Hz) \n");


    printf("\n\n FRMS written to:  FRMS.out ");
	printf("\n The format is: ");
	printf("  freq (Hz)  Cumulative RMS \n");


	printf("\n\n Please call the output files into your own plotting program.");

	printf("\n Press any key to exit.");
	getch();

    exit(1);
}

void read_data()
{

	double dtmin=1.0e+90;

	double dtmax=-dtmin;

    double tb;

	ave=0.;

	rms=0.;


	i=0;


	while( fgets(string,numBytes,pFile[0] )>0)
    {
		kflag=0;

		line_test();

		if(kflag==0)
		{
			if(strrchr(string,',' ) != NULL )
			{
				sscanf(string,"%lf, %lf", &t,&a[i] );
			}
			else
			{
				sscanf(string,"%lf %lf", &t,&a[i] );
			}


			rms+=pow(a[i],2.);  

			ave+=a[i];


			if(i==0)t1=t;

			if(i== (MAX-1))
			{

				printf("\n Warning: maximum line number reached. \n");

				i++;

				break;

			}
			if(i>0)
			{
				if(t<tb)
				{

					printf("\n Time error.  Time values must ascend.");

					printf("\n Press any key to exit \n");

					getch();

					exit(1);

				}
				if(dtmin>(t-tb))
				{
					dtmin=(t-tb);
				}
				if(dtmax<(t-tb))
				{
					dtmax=(t-tb);
				}
			}

			tb=t;

			i++;

		}
	}
	fclose( pFile[0]);


	last = i;

	n = i;

	nm1=n-1;


    tmax=t-t1;

    ts=tmax/nm1;

    fs=1./ts;

    df=1/tmax;      

    rms/=n;

	rms=sqrt(rms);

	ave/=n;

	dt=ts;

	sr=1./dt;


	printf("\n Time span: %10.5g to %10.5g sec \n",t1,t);

	printf("\n dtmin= %11.6g sec	   ",dtmin);
	printf("\n dt   = %11.6g sec	   ",dt);
	printf("\n dtmax= %11.6g sec	 \n",dtmax);

	printf("\n srmin= %11.6g samples/sec	   ",1./dtmax);
	printf("\n sr   = %11.6g samples/sec	   ",sr);
	printf("\n srmax= %11.6g samples/sec	 \n",1./dtmin);


	if( dtmin < 0.99 *dtmax)
	{
	   int nsr;

	   printf("\n *** Time Step Warning ***   \n");

       printf("\n dtmax = %8.4g   dtmin = %8.4g \n",dtmax,dtmin);

	   printf("\n The time step must be constant. \n");

       printf("\n Enter new sample rate?  1=yes  2=no \n");

	   scanf("%d",&nsr);


	   if(nsr==1)
       {
			printf("\n\n Assume constant sample rate. Enter sample rate (samples/sec) \n");

			scanf("%lf",&sr);

			dt=1./sr;
	
			ts=dt;

			fs=sr;


			printf("\n\n sr   = %11.6g samples/sec	   ",sr);

			printf("\n Nyquist Freq = %lf Hz \n",sr/2.);

			printf("\n dt   = %11.6g sec	   \n",dt);

	   }
	   else
	   {

			printf("\n Press any key to exit \n");

			getch();

			exit(1);

	   }
	}

}
void elapsed(long *k)
{
   
	time( &ltimen );

	tnow=double(ltimen);

 //   printf("\n %ld %ld %ld ",*k, ltimen,tstart);

    irem=long( double(n-1)*(tnow-tstart)/double(*k)  );      

	if(  (tnow-tstart) > 20.)
	{
  	  long ihr= long( (irem                      ) / 3600.);
      long imin=long( (irem - ihr*3600           ) / 60.  );
      long isec=long( (irem - ihr*3600 - imin*60 ) / 1.   );

//	  printf("\n k= %ld\n",*k);
	  
      printf("\n\n Estimated remaining time (hr:min:sec)\n");
				
      printf("      = %ld:%ld:%ld \n", ihr,imin,isec );

	  iflag = 1;
	}
}
void hanning(const long &last)
{
	double ae=sqrt(8./3.);

	double dt=ts;

	strcpy(filename[5],"hanning.out");
	pFile[5]=fopen(filename[5],"w");

	printf("\n Hanning window ");

	for(i=0; i < last; i++)
    {
	    a[i]*=float( ae*pow(sin( (i*pi/last) ),2) );

		t=t1+ i*dt;

		fprintf(pFile[5],"%12.5e \t %10.4e \n",t,a[i]);
	}

	fclose( pFile[5]);
}

void meanremove(const double &ave,const long &last)
{

 //   printf("\n ave = %12.4g \n",ave);

	printf("\n mean removal");

	for(i=0; i < last; i++)
    {
	    a[i]-=float(ave);
	}

}
void line_test()
{


//		printf("j=%ld.  %s\n",j,string);



		if( strrchr(string,'A' ) != NULL ){kflag=1;}
		if( strrchr(string,'B' ) != NULL ){kflag=2;}
		if( strrchr(string,'C' ) != NULL ){kflag=3;}
		if( strrchr(string,'D' ) != NULL ){kflag=4;}
		if( strrchr(string,'E' ) != NULL && strstr(string,"E+")==NULL && strstr(string,"E-")==NULL ){kflag=5;}
		if( strrchr(string,'F' ) != NULL ){kflag=6;}
		if( strrchr(string,'G' ) != NULL ){kflag=7;}
		if( strrchr(string,'H' ) != NULL ){kflag=8;}
		if( strrchr(string,'I' ) != NULL ){kflag=9;}
		if( strrchr(string,'J' ) != NULL ){kflag=10;}
		if( strrchr(string,'K' ) != NULL ){kflag=11;}
		if( strrchr(string,'L' ) != NULL ){kflag=12;}
		if( strrchr(string,'M' ) != NULL ){kflag=13;}
		if( strrchr(string,'N' ) != NULL ){kflag=14;}
		if( strrchr(string,'O' ) != NULL ){kflag=15;}
		if( strrchr(string,'P' ) != NULL ){kflag=16;}
		if( strrchr(string,'Q' ) != NULL ){kflag=17;}
		if( strrchr(string,'R' ) != NULL ){kflag=18;}
		if( strrchr(string,'S' ) != NULL ){kflag=19;}
		if( strrchr(string,'T' ) != NULL ){kflag=20;}
		if( strrchr(string,'U' ) != NULL ){kflag=21;}
		if( strrchr(string,'V' ) != NULL ){kflag=22;}
		if( strrchr(string,'W' ) != NULL ){kflag=23;}
		if( strrchr(string,'X' ) != NULL ){kflag=24;}
		if( strrchr(string,'Y' ) != NULL ){kflag=25;}
		if( strrchr(string,'Z' ) != NULL ){kflag=26;}

//

		if( strrchr(string,'a' ) != NULL ){kflag=27;}
		if( strrchr(string,'b' ) != NULL ){kflag=28;}
		if( strrchr(string,'c' ) != NULL ){kflag=29;}
		if( strrchr(string,'d' ) != NULL ){kflag=30;}
		if( strrchr(string,'e' ) != NULL  && strstr(string,"e+")==NULL && strstr(string,"e-")==NULL ){kflag=31;}
		if( strrchr(string,'f' ) != NULL ){kflag=32;}
		if( strrchr(string,'g' ) != NULL ){kflag=33;}
		if( strrchr(string,'h' ) != NULL ){kflag=34;}
		if( strrchr(string,'i' ) != NULL ){kflag=35;}
		if( strrchr(string,'j' ) != NULL ){kflag=36;}
		if( strrchr(string,'k' ) != NULL ){kflag=37;}
		if( strrchr(string,'l' ) != NULL ){kflag=38;}
		if( strrchr(string,'m' ) != NULL ){kflag=39;}
		if( strrchr(string,'n' ) != NULL ){kflag=40;}
		if( strrchr(string,'o' ) != NULL ){kflag=41;}
		if( strrchr(string,'p' ) != NULL ){kflag=42;}
		if( strrchr(string,'q' ) != NULL ){kflag=43;}
		if( strrchr(string,'r' ) != NULL ){kflag=44;}
		if( strrchr(string,'s' ) != NULL ){kflag=45;}
		if( strrchr(string,'t' ) != NULL ){kflag=46;}
		if( strrchr(string,'u' ) != NULL ){kflag=47;}
		if( strrchr(string,'v' ) != NULL ){kflag=48;}
		if( strrchr(string,'w' ) != NULL ){kflag=49;}
		if( strrchr(string,'x' ) != NULL ){kflag=50;}
		if( strrchr(string,'y' ) != NULL ){kflag=51;}
		if( strrchr(string,'z' ) != NULL ){kflag=52;}

//
		if( strrchr(string,'!' ) != NULL ){kflag=53;}
		if( strrchr(string,'@' ) != NULL ){kflag=54;}
		if( strrchr(string,'#' ) != NULL ){kflag=55;}
		if( strrchr(string,'$' ) != NULL ){kflag=56;}
		if( strrchr(string,'%' ) != NULL ){kflag=57;}
		if( strrchr(string,'^' ) != NULL ){kflag=58;}
		if( strrchr(string,'&' ) != NULL ){kflag=59;}
		if( strrchr(string,'*' ) != NULL ){kflag=60;}
		if( strrchr(string,'/' ) != NULL ){kflag=61;}
		if( strrchr(string,'>' ) != NULL ){kflag=62;}
		if( strrchr(string,'<' ) != NULL ){kflag=63;}
		if( strrchr(string,'?' ) != NULL ){kflag=64;}
		if( strrchr(string,'=' ) != NULL ){kflag=65;}
		if( strrchr(string,';' ) != NULL ){kflag=66;}

		if( strrchr(string,'1')==NULL && strrchr(string,'2')==NULL && strrchr(string,'3')==NULL && strrchr(string,'4')==NULL  && strrchr(string,'5')==NULL && strrchr(string,'6')==NULL && strrchr(string,'7')==NULL && strrchr(string,'8')==NULL  && strrchr(string,'9')==NULL  && strrchr(string,'0')==NULL){ kflag=67;}

}