// FFT.cpp: implementation of the CFFT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "mon.h"
#include "FFT.h"
#include "math.h"

#define M_PI 3.1415926

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CFFT::CFFT()
{

}

CFFT::~CFFT()
{

}

float Hanning(float n,float j)
{
   float fret;

   fret = 0.5 * (1.0 - cos(2.0 * M_PI * j / (n-1.0)));
   return(fret);
}

void CFFT::WindowData(float *x,int numdat,int window)
{
   int i;
   float multiplier;

   for ( i = 0; i <= numdat - 1; ++i ) {
      switch ( window ) {
/*	 case 0:
	    multiplier = 1.0;
	    break;
	 case 1:
	    multiplier = Parzen(numdat,i);
	    break;
*/	 case 2:
	    multiplier = Hanning(numdat,i);
	    break;
/*	 case 3:
	    multiplier = Welch(numdat,i);
	    break;
	 case 4:
	    multiplier = Hamming(numdat,i);
	    break;
	 case 5:
	    multiplier = ExactBlackman(numdat,i);
	    break;
*/	 default:
	    multiplier = 1.0;
	    break;
      }
      x[i] = multiplier * x[i];
   }
}

/*void WindowFFTData(float *xdata, float *ydata, int numdat, int window)
{
   WindowData(xdata,numdat, window);
   WindowData(ydata,numdat, window );
}
*/

void FFT_swap(float *s1, float *s2)
{
   float temp;

   temp = (*s1);
   (*s1) = (*s2);
   (*s2) = temp;
}

void CFFT::fft(float *xreal, float *yimag, int numdat, int flag)
{
   int maxpower, arg,  cntr, pnt0, pnt1, i;
   int j, a, b, k;
   float sign, prodreal, prodimag, harm;
	float *cosary;
	float *sinary;

	cosary =  new float[numdat];

	sinary =  new float[numdat];

	j = 0;
	if ( flag != 0 ) 
	{
      sign = 1.0;
      for ( i = 0; i <= numdat - 1; ++i ) 
	  {
		xreal[i] = xreal[i] / numdat;
		yimag[i] = yimag[i] / numdat;
      }
	}
   else {
      sign = -1.0;
   }
   for ( i = 0; i <= numdat - 2; ++i ) 
   {
      if ( i < j ) 
	  {
		 FFT_swap(&xreal[i],&xreal[j]);
		 FFT_swap(&yimag[i],&yimag[j]);
      }
      k = numdat / 2;
      while ( k <= j ) 
	  {
		  j = j - k;
		  k = k / 2;
      }
      j = j + k;
   }
   maxpower = 0;
   i = numdat;
   while ( i != 1 ) 
   {
      maxpower = maxpower + 1;
      i = i / 2;
   }
   harm = 2 * M_PI / numdat;
   for ( i = 0; i <= numdat - 1; ++i ) {
      sinary[i] = sign * sin(harm * i);
      cosary[i] = cos(harm * i);
   }
   a = 2;
   b = 1;
   for ( cntr = 1; cntr <= maxpower; ++cntr ) {
      pnt0 = numdat / a;
      pnt1 = 0;
      for ( k = 0; k <= b - 1; ++k ) {
	 i = k;
	 while ( i < numdat ) {
	    arg = i + b;
	    if ( k == 0 ) {
	       prodreal = xreal[arg];
	       prodimag = yimag[arg];
	    }
	    else {
	       prodreal = xreal[arg] * cosary[pnt1] - yimag[arg] * sinary[pnt1];
	       prodimag = xreal[arg] * sinary[pnt1] + yimag[arg] * cosary[pnt1];
	    }
	    xreal[arg] = xreal[i] - prodreal;
	    yimag[arg] = yimag[i] - prodimag;
	    xreal[i] = xreal[i] + prodreal;
	    yimag[i] = yimag[i] + prodimag;
	    i = i + a;
	 }
	 pnt1 = pnt1 + pnt0;
      }
      a = 2 * a;
      b = b * 2;
   }
   delete cosary;
   delete sinary;
}

void CFFT::FFTCalc(float *xreal,float *yimag,int numdat)
{

   fft(xreal,yimag,numdat,0);
}