fft.cpp

快速傅里叶变换

#include <stdio.h>
#include <iostream.h>
#include <math.h>
#define pi 3.1415926
#define N 4

struct Complex
	{
		double real;
		double image;
	};
//struct Complex y[N];

struct Complex comadd(struct Complex a,struct Complex b);
struct Complex comsub(struct Complex a,struct Complex b);
struct Complex commul(struct Complex a,struct Complex b);
void FFT(struct Complex extract[N],int M);

void main()
{
	int i,j;
	struct Complex x[N][N],X[N][N];	
	//暂且生成一个测试数组
	/*for(i=0;i<=N-1;i++)
		for(j=0;j<=N-1;j++)
		{
			x[i][j].real=i+1;
			x[i][j].image=j+1;
		}*/
	x[0][0].real=x[0][1].real=x[0][2].real=x[0][3].real=x[0][0].image=x[1][0].image=x[2][0].image=x[3][0].image=1;
	x[1][0].real=x[1][1].real=x[1][2].real=x[1][3].real=x[0][1].image=x[1][1].image=x[2][1].image=x[3][1].image=2;
	x[2][0].real=x[2][1].real=x[2][2].real=x[2][3].real=x[0][2].image=x[1][2].image=x[2][2].image=x[3][2].image=3;
	x[3][0].real=x[3][1].real=x[3][2].real=x[3][3].real=x[0][3].image=x[1][3].image=x[2][3].image=x[3][3].image=4;

	//对数组的维数进行计算，暂且要求其值为2的整数次幂，并将这个幂求出来。
	double test;
	test=log10(N)/log10(2);
	int M;
	M=int(test);
	if(test!=M)
	{
		cout<<"所设定的数组维数不是2的整数次幂！"<<endl;
		cout<<"请重新设定所要进行变换的数组:"<<endl;
		return;
	}
	cout<<"生成的数组为："<<endl;
	for(i=0;i<=N-1;i++)
		for(j=0;j<=N-1;j++)
		{
			if(x[i][j].image>=0)
			printf("x[%d][%d]=%f+%fi\n",i,j,x[i][j].real,x[i][j].image);
		else
			printf("x[%d][%d]=%f%fi\n",i,j,x[i][j].real,x[i][j].image);
	}

	//对二维数组的每一行进行提取，进行第一次傅里叶变换,变换完后的数组仍然放在x里面。
	struct Complex extractrow[N];
	for(i=0;i<=N-1;i++)
	{
		for(j=0;j<=N-1;j++)
		{
			extractrow[j]=x[i][j];
			//extractrow[j].real=x[i][j].real;
			//extractrow[j].image=x[i][j].image;
		}
		FFT(extractrow,M);
		for(j=0;j<=N-1;j++)
		{
			X[i][j]=extractrow[j];
		}
	}
	
	/*//行变换后的数组为：
	cout<<"行变换后的数组为："<<endl;
	for(i=0;i<=N-1;i++)
		for(j=0;j<=N-1;j++)
		{
			if(X[i][j].image>=0)
			printf("X[%d][%d]=%f+%fi\n",i,j,X[i][j].real,X[i][j].image);
		else
			printf("X[%d][%d]=%f%fi\n",i,j,X[i][j].real,X[i][j].image);
	}*/

//对第一次变换后的数组每一列进行提取，进行第二次傅里叶变换。
	struct Complex extractcol[N];
	for(i=0;i<=N-1;i++)
	{
		for(j=0;j<=N-1;j++)
		{
			extractcol[j]=X[j][i];
			//extractcol[j].real=x[j][i].real;
			//extractcol[j].image=x[j][i].image;
		}
		FFT(extractcol,M);
		for(j=0;j<=N-1;j++)
		{
			x[i][j]=extractcol[j];
		}
	}
	/*// 列变换后的数组为：
	cout<<"列变换后的数组为："<<endl;
	for(i=0;i<=N-1;i++)
		for(j=0;j<=N-1;j++)
		{
			if(x[i][j].image>=0)
				printf("x[%d][%d]=%f+%fi\n",i,j,x[i][j].real,x[i][j].image);
			else
				printf("x[%d][%d]=%f%fi\n",i,j,x[i][j].real,x[i][j].image);
		}*/
	//对变换后的数组进行转置,这样一来，结果就和matlab里面的fft2变换的结果完全相同了
	for(i=0;i<=N-1;i++)
		for(j=0;j<=N-1;j++)
		{
			X[i][j]=x[j][i];
		}
	cout<<"转置后的矩阵,也就是最后的变换结果为："<<endl;
		for(i=0;i<=N-1;i++)
			for(j=0;j<=N-1;j++)
			{
				if(X[i][j].image>=0)
					printf("X[%d][%d]=%f+%fi\n",i,j,X[i][j].real,X[i][j].image);
				else
					printf("X[%d][%d]=%f%fi\n",i,j,X[i][j].real,X[i][j].image);
			}
}
			
			


void FFT(struct Complex extract[N],int M)
{
	int i,j;		
	int *temp0=new int[M];
	int temp,temp1;
	struct Complex X[N];
	for(i=0;i<=N-1;i++) 
	{
		temp=0;
		temp1=i;
		for(j=0;j<=M-1;j++)
		{
			temp0[j]=temp1%2;
			temp1=int(temp1/2);			
			temp=temp+int(temp0[j]*pow(2,(M-j-1)));			
		}
		X[i]=extract[temp];
		//printf("X[%d]=x[%d]  ",temp,i);				
		cout<<endl;
	}
	delete []temp0;

	/*//对倒位序后的数组进行输出
	for(i=0;i<=N-1;i++)
	{
		if(x[i].image>=0)
			printf("X[%d]=%f+%fi\n",i,X[i].real,X[i].image);
		else
			printf("X[%d]=%f%fi\n",i,X[i].real,X[i].image);
	}*/

	//开始进行傅里叶变换
	int k,m,B,J,P;
	Complex W,tempX;//tmepXB;
	for(m=1;m<=M;m++)
	{
		B=int(pow(2,(m-1)));

		for(J=0;J<=B-1;J++)
		{
			P=int(J*pow(2,(M-m)));
			W.real=cos(2*pi*P/N);  //计算旋转因子的实部
			W.image=-sin(2*pi*P/N); //计算旋转因子的虚部

			for(k=J;k<=N-1;k=int(k+pow(2,m)))
			{							

	
				tempX=X[k];
				//tmepXB=X[k+B];
				X[k]=comadd(tempX,commul(X[k+B],W));
				X[k+B]=comsub(tempX,commul(X[k+B],W));
				//X[k].real=tempX.real+tmepXB.real*W.real-tmepXB.image*W.image;
				//X[k].image=tempX.image+tmepXB.real*W.image+tmepXB.image*W.real;

				//X[k+B].real=tempX.real-(tmepXB.real*W.real-tmepXB.image*W.image);
				//X[k+B].image=tempX.image-(tmepXB.real*W.image+tmepXB.image*W.real);
			}
		}
	}
	for(i=0;i<=N-1;i++)
	{
		extract[i].real=X[i].real;
		extract[i].image=X[i].image;
	}

	//变换后的数组为：
/*	cout<<"变换后的数组为："<<endl;
	for(i=0;i<=N-1;i++)
	{
		if(X[i].image>=0)
			printf("X[%d]=%f+%fi\n",i,X[i].real,X[i].image);
		else
			printf("X[%d]=%f%fi\n",i,X[i].real,X[i].image);
	}*/
	
}

struct Complex comadd(struct Complex a,struct Complex b)
{
	struct Complex c;
	c.real=a.real+b.real;
	c.image=a.image+b.image;
	return c;
}

struct Complex comsub(struct Complex a,struct Complex b)
{
	struct Complex c;
	c.real=a.real-b.real;
	c.image=a.image-b.image;
	return c;
}
struct Complex commul(struct Complex a,struct Complex b)
{
	struct Complex c;
	c.real=a.real*b.real-a.image*b.image;
	c.image=a.real*b.image+a.image*b.real;
	return c;
}