fft.cpp

FFT 变换 该程序可以实现运行时输入时域的采样数据

# include <math.h>
# include <malloc.h>
# include <memory.h>
# include <iostream.h>
# include <stdio.h>
# include <stdlib.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;//代表2的power幂次方，即2^power
	cout<<"input the power"<<endl;
	cin>>power;
	int m;
	m=pow(2,power);

	COMPLEX *T=new COMPLEX [m];
	COMPLEX *F=new COMPLEX [m];
    COMPLEX *TD=T;
	COMPLEX *FD=F;
//	cout<<(int *)T<<endl;
	
	for(int i=0;i<m;i++)//这个地方设置时域的实部和虚部!!
	{
		cout<<"输入第"<<(i+1)<<"个时部："<<endl;//可以在这里输入时域
		cin>>TD[i].re;
	//	TD[i].re=exp(-0.1*(i+0.5));//可以在这里自己设定步长
		TD[i].im=0.0;
	}

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

	

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

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

	cout<<"/////////////////////////////////////////////////"<<endl;
	cout<<"频域："<<endl;
	for(int i3=0;i3<m;i3++)
	{
		cout<<FD[i3].re<<"\t"<<FD[i3].im<<endl;
	}
//cout<<(int *)T<<endl;
	FILE *fp;
	if((fp=fopen("a.txt","w"))==NULL)
	{
		printf("can't open file\n");
		exit(1);
	}
    char ch1[10]="时域数据\n";
	char ch2[10]="频域数据\n";
	fprintf(fp,"%s",ch1);
	for( i=0;i<m;i++)
	{
		fprintf (fp, "%lf",T[i].re);
        fprintf (fp, " %lf \n",T[i].im);
	   
	}
	 fprintf(fp,"%s",ch2);
	for( i=0;i<m;i++)
	{
		fprintf (fp, "%lf",F[i].re);
        fprintf (fp, " %lf \n",F[i].im);
	   
	}

//	fwrite(T,sizeof(double),2*m,fp);
	fclose(fp);
	
	//cout<<sizeof(COMPLEX)<<endl;
  delete []T;
  delete []F;
}