fft.cpp

FFT快速傅立叶算法程序

# include <math.h>
# include <malloc.h>
# include <memory.h>
# include <iostream.h>
//////////////////////////////////////////////////////////
// internal definitions

#define PI (double)3.14159265359

/*complex number*/
typedef struct
{
	double re;
	double im;
}COMPLEX;

/*complex add*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}

/*complex substract*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}

/*complex multiple*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}
//////////////////////////////////////////////////////////

/*
void FFT(COMPLEX * TD, COMPLEX * FD, int power);
void IFFT(COMPLEX * FD, COMPLEX * TD, int power);
void DCT(double *f, double *F, int power);
void IDCT(double *F, double *f, int power);
void WALh(double *f, double *W, int power);
void IWALh(double *W, double *f, int power);
*/

/****************************************************
	FFT()

	参数：

		TD为时域值
		FD为频域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数实现快速傅立叶变换
****************************************************/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
	int count;
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W,*X1,*X2,*X;

	/*计算傅立叶变换点数*/
	count=1<<power;
	
	/*分配运算所需存储器*/
	W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	
	/*计算加权系数*/
	for(i=0;i<count/2;i++)
	{
		angle=-i*PI*2/count;
		W[i].re=cos(angle);
		W[i].im=sin(angle);
	}
	
	/*将时域点写入存储器*/
	memcpy(X1,TD,sizeof(COMPLEX)*count);
	
	/*蝶形运算*/
	for(k=0;k<power;k++)
	{
		for(j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-k);
			for(i=0;i<bfsize/2;i++)
			{
				p=j*bfsize;
				X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
				X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
			}
		}
		X=X1;
		X1=X2;
		X2=X;
	}
	
	/*重新排序*/
	for(j=0;j<count;j++)
	{
		p=0;
		for(i=0;i<power;i++)
		{
			if (j&(1<<i)) p+=1<<(power-i-1);
		}
		FD[j]=X1[p];
	}
	
	/*释放存储器*/
	free(W);
	free(X1);
	free(X2);
}


void main(void)
{
	int power=6;
	COMPLEX TD[64],FD[64];
	for(int i=0;i<64;i++)
	{
		TD[i].re=exp(-0.1*(i+0.5));
		TD[i].im=0.0;
	}

	cout.setf(ios::fixed);
	cout.precision(8);

	cout << "序列TD[i]为：  "<<endl;
    for(int i1=0;i1<16;i1++)
	{
		for(int j=0;j<4;j++)
		{
			cout<<TD[i1*4+j].re<<"\t"<<TD[i1*4+j].im<<"\n";
		}
		
	}

	//调用fft函数；
	FFT  (TD,   FD,  power);

	cout<<"/////////////////////////////////////////////////"<<endl;
	for(int i2=0;i2<64;i2++)
	{
		cout<<TD[i2].re<<"\t"<<TD[i2].im<<endl;
	}

	cout<<"/////////////////////////////////////////////////"<<endl;
	for(int i3=0;i3<64;i3++)
	{
		cout<<FD[i3].re<<"\t"<<FD[i3].im<<endl;
	}


	cout<<"某////////////////////////////////////////////////"<<endl;
	for(int i4=0;i4<64;i4++)
	{
		cout<<pow(pow(FD[i4].re,2)+pow(FD[i4].im,2),0.5)<<endl;
	}

	cout<<"幅角//////////////////////////////////////////////"<<endl;
	for(int i5=0;i5<64;i5++)
	{
		cout<<atan2(FD[i5].im,FD[i5].re)*180/PI<<endl;
	}

//	cout<<atan(1)<<endl;
}