fft.cpp

用visual c++6.0编写的 FFT程序

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

#include "stdafx.h"
#include "zhineng.h"
#include "FFT.h"
#include "iostream.h"
#include "math.h"


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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
#define  pi  3.141592653
MultiNumber FFT::one = {-1,0};
FFT *FFT::myFFT=NULL;//2008.9.10李琳

FFT::FFT()
{
	myFFT=this;//2008.9.10李琳
	rest =0;////////问题
// 	k =0;////////////问题
    for (int i=0;i<MAXCHAR;i++)
    {
		ashow[i]=0;
    }
}

FFT::~FFT()
{

}


MultiNumber FFT::ComMultiply(MultiNumber x , MultiNumber y)
{
	MultiNumber a;

    a.real = x.real*y.real - x.imag*y.imag;

    a.imag = x.real*y.imag + x.imag*y.real;

	return a;
}

MultiNumber FFT::ComAdd(MultiNumber x , MultiNumber y)
{
	MultiNumber a;

	a.real = x.real + y.real;

	a.imag = x.imag + y.imag;

	return a;
}

MultiNumber FFT::W(int m , int k)
{
	MultiNumber a;

    double b = pow(2,double(m-1));

	double c = double(k%(int(b)));

	b = pow(2,m);

	a.real = cos(2*pi*c/b);

	a.imag = (-1)*sin(2*pi*c/b);

	return a;
}


int FFT::GetRest(int v , int k , int m)

{
	int rest =0;

	int t1 = 0;

	int t2 = 0;

	int t3 = 0;

	int k1 = k;

	for(int d = 1 ; d<=m-1 ; d++)

	{
		t1 = int(pow(2 , double(v-d)));

		t2 = int(pow(2 , double(d-1)));

        t3 = k/t1;
 
        rest = t3*t2+rest;

		k = k - t1*(t3);


	}

	return rest;

}


MultiNumber *FFT::FFT1(int v , MultiNumber *h)

{

int static N = int(pow(2,double(v)));

int step = 1;

MultiNumber a ,b;

for(int m = 1; m<=v ; m++)

{


step = int(pow(2 , double(v-m)));

for(  int i = 0 ; i<N-step  ; i = i + step*2)

{   

	for( int k=i; k<=i+step-1 ; k++ )

{

	a = ComAdd(h[k+0] , ComMultiply(h[k+step] , W(m,GetRest(v,k,m))));

	b = ComAdd(h[k+0] , ComMultiply(h[k+step] ,ComMultiply(one, W(m,GetRest(v,k,m)))));

	h[k+0] = a;

	h[k+step] = b;

    
}

}

}

return h;



}

int FFT::InverseBit(int k , int v)

{

	int t = 1;

	int t1 =1;

	int k1 = k;

	int result = 0;

for(int i=v-1 ; i>=0 ; i--)

{

  t = int(pow(2,i));

  t1 = int(pow(2,v-i-1));

  result = result + (k1/t)*t1;

  k1 = k1 - (k1/t)*t;

}

return result;

}


MultiNumber *FFT::GetLastRe(MultiNumber *h , int v)

{
    int y = 0;

	int N = int(pow(2,v));

	MultiNumber *h1 = new MultiNumber[N];

	for(int i=0 ; i<N ; i++)

	{

        y = InverseBit(i,v);

        h1[i] = h[y];

	}

	return h1;

}



int FFT::mainFFT()
{   


	int m = 9;
    
	int N = int(pow(2, m));
    MultiNumber *a= new MultiNumber[N];//////////////////问题
///////////////////////////2008.12.22李琳
	HANDLE hSearch; 
	BOOL fFinished = FALSE; 
	char szFilePath[MAXCHAR];
	int pathlen;
	pathlen=::GetModuleFileName(NULL,szFilePath,MAX_PATH); 
	while(TRUE)   
	{   
		if   (szFilePath[pathlen]!='\\')
		{szFilePath[pathlen]=NULL;
		pathlen--;}
		else
			break;   
	}   	
	strcat(szFilePath,"code\\006modify.txt"); //2008.12.11文件名李琳

	FILE* stlfp;
	CString fileName;
	CString tempstring;
	int line=0;
	char readline[100];
	double x=0;
	fileName=szFilePath;
	stlfp = fopen(fileName,"r");
	if(!stlfp)
		return false;
	else
	{
		do
		{
			fgets(readline,100,stlfp);
			tempstring.Format("%s", readline);  
			a[line].real = atof(tempstring);
			a[line].imag = 0;
// 			ashow[line]=a[line].real;
			line++;
		}while(!feof(stlfp)&&line<400);//////////数目不能超过479，否则溢出，李琳2008.12.29解决问题
		fclose(stlfp);
	}

//     for (int acount=0;acount<line;acount++)
//     {
// 		ashow[acount]=a[acount].real;
//     }

	MultiNumber *b = FFT1(m,a);

	b = GetLastRe(b,m);

	for( int i =0 ; i<N ; i++)

	{
    	 cout<<b[i].real<<"  "<<b[i].imag<<endl; 
	}
//     for (int bcount=0;bcount<N;bcount++)
//     {
// 		bshow[bcount]=b[bcount].real;/////////////问题
//     }
	return 0;
}