fft1.cpp

fft变换

#include <stdio.h>
#include <math.h>
#define N 32
#define PI 3.1415926

struct compx{
	double real;
	double imag;
};

struct compx addcom(struct compx v1,struct compx v2)
{
	struct compx v3;
	v3.real=v1.real+v2.real;
	v3.imag=v1.imag+v2.imag;
	return(v3);
}

struct compx minuscom(struct compx v1,struct compx v2)  //定义复数相减
{
	struct compx v3;
	v3.real=v1.real-v2.real;
	v3.imag=v1.imag-v2.imag;
	return(v3);
}

struct compx EE(struct compx v1,struct compx v2)
{
	struct compx v3;
	v3.real=v1.real*v2.real-v1.imag*v2.imag;
	v3.imag=v1.real*v2.imag+v1.imag*v2.real;
	return(v3);
}

struct compx FF(struct compx v,int n) //定义复数的n次方
{
	struct compx v1;
	int i;
	v1.real=1.0;v1.imag=0.0;
	for(i=1;i<=n;i++)
		v1=EE(v1,v);
	return(v1);
}


/*下面进行fft运算按时间抽取,先做旋转因子后做蝶形图*/
void FFT(struct compx xin[])
{
	int i,j,l,k,m,f;
	double p,ps;
	int le,lei,ip;
	struct compx w,t,tmp;
	
	f=N;
	for(m=1;(f=f/2)!=1;m++) {;}  //计算出共m级蝶形图,避免了对数运算
	
	for(l=0;l<m;l++)
	{
		le=(int)pow(2,m-l);
		lei=le/2;
		for(j=0;j<=lei-1;j++)
		{
			p=pow(2,l)*j;
			ps=2*PI/N*p;
			w.real=cos(ps);w.imag=-sin(ps);  //w即为旋转因子
			for(i=j;i<N;i=i+le)
			{
				ip=i+lei;
				t=xin[i];
				xin[i]=addcom(xin[i],xin[ip]);
				xin[ip]=minuscom(t,xin[ip]);
				xin[ip]=EE(xin[ip],w);
			}	
		}
	}


	j=N/2;
	for(i=1;i<N-1;i++)
	{
		if(i<j) {tmp=xin[i];xin[i]=xin[j];xin[j]=tmp;}
		k=N/2;
		while(j>=k) {j=j-k;k=k/2;}
		j=j+k;
	}
}







int main()
{
	int i,j;
	struct compx q,x[N];
	void FFT(struct compx x[]);
	q.real=0.9;q.imag=0.3;
	for(i=0;i<N;i++)
	{
		x[i]=FF(q,i);
		printf("x[%d] is %13.7f+j%13.7f\n",i,x[i].real,x[i].imag);
	}
	FFT(x);
	for(i=0;i<N;i++)
		printf("y[%d] is %13.7f+j%13.7f\n",i,x[i].real,x[i].imag);
	return 0;
}