/*

      helicopter_engine_FFT.cpp

      by Tom Irvine
      

*/

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

//   262144 
//   131072 
//    65536 
//    32768
//    16384

#define MAX   262144
#define MAXP1 131072

//************ complex number handling ********************
typedef struct 
{
    double r,i;
	
}doublecomplex;

typedef struct 
{
    float r,i;
	
}floatcomplex;

#define  COMPLEXMAG(x,y,z) z=sqrt( pow( x,2) + pow(y,2) );
#define COMPLEXMAG2(x,y,z) z=    ( pow( x,2) + pow(y,2) );

#define COMPLEXRECIP(ar, ai, xr, xi)              \
{                                                 \
	   xr=  ar/( pow( ar,2) + pow(ai,2) );    \
	   xi= -ai/( pow( ar,2) + pow(ai,2) );    \
}                                                 \

#define COMPLEXMUL(ar, ai, br, bi,cr,ci)          \
{                                                 \
	   cr=  ar*br -ai*bi;                     \
	   ci=  ar*bi +ai*br;                     \
}                                                 \

#define COMPLEXDIV(nr, ni, dr, di,cr,ci)          \
{                                                 \
	   cr =  nr*dr +ni*di;                    \
	   cr/= ( pow( dr,2) + pow(di,2) );       \
	   ci =  ni*dr -nr*di;                    \
	   ci/= ( pow( dr,2) + pow(di,2) );       \
}                                                 \

//*********************************************************

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

floatcomplex y[MAX];
floatcomplex ww[MAXP1];
    
double ae, ave;
double df, dt, dur;
double grms;
double kurt; 
double pi;
double sr, std;
double time, time1, t, tp, tt;

unsigned long i, j, k, n, ns;
unsigned long ich;

unsigned long nstages;
unsigned long nsegment;

double fstart;

int iflag, imr, ihw;

//**** function prototypes ****

void files(void);
void read2(void);
void bitrev(void);
void fourier(void);
void meanremove(void);
void hanning(void);
void ctfourier(void);

void main()
{
	iflag=0;
	pi=atan2(0.,-1.);
	tp=2.*pi;

    for(i=0; i< MAX; i++)
    {
		 y[i].r=float(0.);
	     y[i].i=float(0.); 
	}
	for(i=0; i< MAXP1; i++)
    {
		 ww[i].r=float(0.);
	     ww[i].i=float(0.); 
	}


//***** open files ****************
    
	files();
    
	
//*** read time history

	read2();
    
	if( iflag < 900 )
    {
//***  read mean removal and window parameters

/*
		printf("\n");

		printf("\n Mean removal? 1=yes 2=no \n");
		scanf("%d",&imr);

		printf("\n");

		printf("\n Window? 1=rectangular 2=Hanning \n");
		scanf("%d",&ihw);
*/
		imr=1;
		ihw=2;

//*** perform mean removal if selected.

	   if(imr==1){ meanremove();}

//*** perform hanning window if selected.

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

//****** Perform bit reversal *******

	   bitrev();

//******* Cooley-Tukey Method ******

	   for(i=0; i< MAXP1; i++)
       {
		  ww[i].r=float(0.);
	      ww[i].i=float(0.); 
	   }
       printf("\n call function ctfourier ");
       ctfourier();  
    }
	return;
}
void ctfourier(void)
{
	int i1;

	unsigned long ig;
	
	int iy=1;
	int jj;
	int jjj;
    int l;
	int lg;
	int lk;
	int lwd2;

	unsigned long m;

	int nworld;
	int wd2;
	int world=4;

	double arg;
	double freq;
	double z;

    float am;

	floatcomplex asum;
	floatcomplex a;
	floatcomplex b;
	floatcomplex c;

//      printf("\n  in ctfourier1 ");

    m=nsegment;

	printf("\n m= %ld  nstages= %ld \n\n",m, nstages );

    for(j=0; j<= m-2; j+=2)
    {
	     jj=j+1;
		 
		 a.r=  y[j].r;
	     a.i=  y[j].i;
		 b.r= y[jj].r;
		 b.i= y[jj].i;
		 
		 y[j].r = a.r + b.r;
		 y[j].i = a.i + b.i;

		 y[jj].r = a.r - b.r;
		 y[jj].i = a.i - b.i;
	}
	world=4;

//      printf("\n begin world loop ");
	
    for(i=2; i<=nstages; i++)
	{ 

        wd2=world/2;

	    nworld=int(m/world);    
	  
	    printf("\n i=%ld / %ld ",i,nstages );
//       printf("\n  (wd2-1)=%ld ", (wd2-1) );

        for(lk=0; lk<= (wd2-1); lk++)
		{
	       arg=tp*double(lk)/world;
		 
		    ww[lk].r =float(cos(arg));
		    ww[lk].i =float(-sin(arg));

	    }  
	    lwd2=0;

 //       printf("\n  nworld=%ld ",nworld); 
	    for( jjj=1; jjj<= nworld; jjj++)
		{  
		       i1=1+lwd2;
	      
			   for(lg=1; lg <= wd2; lg++)      
			   {       
					l=i1;
					lwd2=l+int( pow(2,iy) );
		     
					a.r=y[l-1].r;
					a.i=y[l-1].i;
					b.r=y[lwd2-1].r;
					b.i=y[lwd2-1].i;   
					 
					COMPLEXMUL(ww[lg-1].r,ww[lg-1].i,b.r,b.i,c.r,c.i)

					y[l-1].r  = a.r + c.r;
			        y[l-1].i  = a.i + c.i;

					y[lwd2-1].r  = a.r - c.r;
			        y[lwd2-1].i  = a.i - c.i;
					 
					 i1++;
				}
	    }
	    iy++;
	    world*=2;

//                printf("\n iy=%ld world=%ld", iy, world);
	   }

//*** am is the number of points in the segment
	
	   am=float(m);
    
//         printf("\n am=%f ",am);

       asum.r=float(0.);
	   asum.i=float(0.);

//         printf("\n m=%d ",m);

       double vms=0.;
	   double amag;
	   double vmag;
	   double vmag_rms;

       for(ig=0; ig< m; ig++)      
       {
	         y[ig].r/=am;
			 y[ig].i/=am;

			 freq=ig*df;

			 fprintf( pFile[2],"%15.6e %15.6e %15.6e \n", freq, y[ig].r, y[ig].i);
			 
			 
			 COMPLEXMAG(y[ig].r,y[ig].i,z)
			 
			 fprintf( pFile[3],"%15.6e %15.6e \n", freq, z);



	         if(ig<(m/2))  
			 {
				 asum.r += y[ig].r;
				 asum.i += y[ig].i; 

				 if(ig == 0)
				 {
					fprintf( pFile[4],"%15.6e %15.6e \n", freq, z);
				 }
				 else
				 {
					fprintf( pFile[4],"%15.6e %15.6e \n", freq, 2.*z);

					amag=2.*z;

					amag*=386.;

					vmag=amag/(2.*pi*freq);

					vmag_rms = vmag*0.707;

					if(freq >=fstart)
					{
						vms+=pow(vmag_rms,2.);
					}
				 }
			 }
       }

       double oav = sqrt(vms);

	   printf("\n\n Overall velocity = %8.4g in/sec RMS  above 200 Hz \n",oav);

		fclose( pFile[2] );
		fclose( pFile[3] );
		fclose( pFile[4] );


//     printf("\n\n sum real=%lf   sum imag=%lf",asum.r,asum.i);
       
/*		
	   printf("\n\n Complex output data written to:  fft.out \n");
	   printf("\n The format is: ");
	   printf("\n  freq (Hz)  real accel(G)  imaginary accel(G)\n");
	   
	   
	   printf("\n\n Double-sided, half-amplitude magnitude output data written to:  fft_half.out \n");
	   printf("\n The format is: ");
	   printf("\n  freq (Hz)  accel(G) \n");
*/

	   printf("\n\n Single-sided, full-amplitude magnitude output data written to:  fft_full.out \n");
	   printf("\n The format is: ");
	   printf("\n  freq (Hz)  accel(G) \n");


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

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

	   
	   
}
void read2(void)
{
    double sum, sum2, sum4, tn;

	printf("\n\n reading time history ");
	
	if( ( fscanf(pFile[1]," %lf %f ",&time1, &y[0].r) ) <= 0 )
    {
		printf("\n fscanf error ");
    }
	sum=     y[0].r;
    sum2=pow(y[0].r,2);
    sum4=pow(y[0].r,4);

    i=1;
	while ( fscanf(pFile[1]," %lf %f ",&tn, &y[i].r) > 0 )
    {
		sum+=     y[i].r;
		sum2+=pow(y[i].r,2);
		sum4+=pow(y[i].r,4);

		i++;
/*               if(i == nsegment)
		 { 
			 time=tn;
			 break;
		 }       */

		 if( i == MAX){break;} 
	}                               
	ns=i;

	for(i=1; i<=18; i++)
	{
		if( pow(2,i) > ns )
		{
			break;
		}
	}
	nsegment = unsigned long( pow(2, i-1 ) );
	nstages=i-1;

	rewind( pFile[1] );
	i=0;
	while ( fscanf(pFile[1]," %lf %f ",&time, &y[i].r) > 0 )
    {
		i++;
		if(i==nsegment){break;}
	}
	ns=i;
	i=ns;

	printf("\n\n  %ld samples will be used for FFT\n", nsegment);

	if(iflag != 905)
    {
	if( ns > 1)
    {
	    dt=(time-time1)/float(ns-1);
	    sr=1./dt;

		dur = (time-time1);
	    df=1./dur;
	
		ave=sum/i;
		grms=sqrt(sum2/i);
		std = sqrt( pow(grms,2) - pow(ave,2) );
		kurt=sum4/(i*pow(std,4));

		printf("\n samples= %ld ", i);
		printf("\n dt = %10.5g sec  sr= %12.6e samples/sec", dt, sr);
		printf("\n\n time span= %12.6g sec to %12.6g sec", time1, time);
   
		printf("\n\n  ave= %9.4g G    std dev = %9.4g G   rms= %9.4g G",
				    ave, std, grms);
		printf("\n kurtosis= %9.4g", kurt );
		printf("\n");

		printf("\n delta f= %14.4g ",df);

	}       
	else
    {
		 printf("\n  dt error ");
	}
	}
	fclose( pFile[1] );
}

void meanremove()
{
	printf("\n mean removal");
	for(i=0; i < nsegment; i++)
    {
	    y[i].r-=float(ave);
	}
}

void hanning()
{
	ae=sqrt(8./3.);

	printf("\n Hanning window ");
	for(i=0; i < nsegment; i++)
    {
	    y[i].r*=float( ae*pow(sin( (i*pi/ns) ),2) );
	}
}

void files(void)
{
	
//      strcpy(filename[0], "fft.in");
//   pFile[0]=fopen(filename[0], "rb");

   printf("\n  helicopter_engine_FFT - version 1.1, May 10, 2005 ");
   printf("\n  by Tom Irvine ");
   printf("\n  Email:  tomirvine@aol.com \n");

   printf("\n This program calculates the Fast Fourier transform ");
   printf("\n of an acceleration time history. \n");
   printf("\n It then calculated the overall velocity.\n");

   printf("\n The input file must be time(sec) and accel(G) ");
   printf("\n The format is free, but no header lines allowed.\n");
   printf("\n Please enter the input filename. \n");

   scanf("%s", filename[1]);
   pFile[1]=fopen(filename[1], "rb");

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

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

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


   for(i=1; i<=4; i++)
   {
      if(pFile[i] == NULL )
      {
			printf(" Failed to open file: %s \n", filename[i]);
			fclose(pFile[i]);
      }
      else
      {
			printf(" File: %s opened. \n", filename[i]);
      }
   }
   printf("\n Enter the highpass frequency (Hz) \n");
   scanf("%lf",&fstart);

} 

void bitrev(void)
{
	   float dd;
 
	   unsigned long ia, ib, ik, ilk;
       unsigned long kk;
	   unsigned long nm1;
	   unsigned long m, md2, mm1;
	  
	   m=nsegment;
       mm1=m-1;
       dd=float( nstages );
       nm1=unsigned long(dd);
       md2=m/2;

       printf("\n nm1= %ld  md2=%ld  ns=%ld ", nm1, md2, ns);
	   
       for( i=1; i<= nm1; i++)   
       { 

		  for(ik=0; ik<m; ik++)
		  {
		    ww[ik].r=y[ik].r;
		  }
	  
	  kk=0;
	  ia=0;

	  ilk=unsigned long( pow(2,(i-1)) );

	  for(j=1; j<=ilk; j++)
		  {
			  ib=ia+md2-1;
	     
	     for(k=ia; k<=ib; k++)
			 {
		if(k==ia)
				{
					kk=k;
				}
		
		y[k].r = ww[kk].r;
		y[k].i = float(0.);
		
				y[k+md2].r = ww[kk+1].r;
				y[k+md2].i =float(0.);
		
				kk=kk+2;

	     }
	     ia=ia + md2*2;
	  
	  }
	  md2=(md2/2);   
		  
	}
 }