fft.cpp

快速傅立叶算法

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

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

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

#define pi 3.1415927
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CFFT::CFFT()
{
	pComplex = NULL;
}

CFFT::~CFFT()
{
	if ( pComplex )
		delete [] pComplex;
}

void CFFT::b2fft(complex x[], int n, int sign)
{
	/*---------------------------------------------------------------------
 Routine mcmpfft:to obtain the DFT of Complex Data x(n)
                               By Cooley-Tukey radix-2 DIT Algorithm .
 input parameters:
 x : complex array.input signal is stored in x(0) to x(n-1).
 n : the dimension of x and y.
 isign:if ISIGN=-1 For Forward Transform
       if ISIGN=+1 For Inverse Transform.
 output parameters:
 x : complex array. DFT result is stored in x(0) to x(n-1).
 Notes:
     n must be power of 2.
                                       in Chapter 5
--------------------------------------------------------------------*/
	complex t;
	double ths,cth,sth;
    double pisign;
    int mr,m,l,j,i,nn;
    for(i=1;i<=16;i++)
	{
		nn=(int)pow(2,i);
	    if(n==nn) break;
	}
	if(i>16)
	{
		printf(" N is not a power of 2 ! \n");
		return;
	}
	pisign=2*sign*pi;
    mr=0;
    for(m=1;m<n;m++)
	{
		l=n;
	    while( (mr+l)>=n )
		
			l=l/2;
            mr=mr%l+l;
		
        if(mr<=m)
            continue;
        t.real=x[m].real;t.imag=x[m].imag;
        x[m].real=x[mr].real;x[m].imag=x[mr].imag;
        x[mr].real=t.real;x[mr].imag=t.imag;
     }
	/*
	CString str;
	for ( m=0; m<8; m++)
	{
		str.Format ( "%f ",x[m].real );
		AfxMessageBox(str);
	}
	*/
/*-------------------------------------------------------------------*/
	l=1;
	while(1)
	{
	  if(l>=n)
	  {
	      if(sign==-1)
				return;
	       for(j=0;j<n;j++)
		   {
				x[j].real=x[j].real/n;
				x[j].imag=x[j].imag/n;
		    }
	       return;
	   }
	  for( m=0; m<l; m++ )
	  {
		for( i=m; i<n; i=i+2*l )
		{
		     ths=m*pisign/(2*l);
		     cth = cos (ths);
			 sth = sin (ths);			 
		     t.real=x[i+l].real*cth-x[i+l].imag*sth;
		     t.imag=x[i+l].real*sth+x[i+l].imag*cth;
		     x[i+l].real=x[i].real-t.real;
		     x[i+l].imag=x[i].imag-t.imag;
		     x[i].real=x[i].real+t.real;
		     x[i].imag=x[i].imag+t.imag;
		}
	  }
	  l=2*l;
	}
//do the forl change and version
}

complex CFFT::cexp(complex a)
{
	complex z;
	z.real=cos(a.imag);
	z.imag=sin(a.imag);
	return(z);
}

void CFFT::RealFFT(double realData[], int length, int sign)
{
	int i;
	if ( pComplex )
		delete [] pComplex;
	pComplex = new complex[length];
	for ( i=0; i<length; i++ )
	{
		pComplex[i].real = realData[i];
		pComplex[i].imag = 0.0;
	}
	this->b2fft( pComplex, length,sign );
	for( i=0; i<length; i++ )
		realData[i] = this->mabs ( pComplex[i] );
	//return the frequency in the realData[];
}

double CFFT::mabs(complex a)
{
	double m;
	m=a.real*a.real+a.imag*a.imag;
	m=sqrt(m);
	return(m);	
}