dft.cpp

实现FFT功能的强大源码

#include<math.h>
#include"complex.h"
#include<iostream>
#include<fstream>
#include<ctime>
using namespace std;
void buffer(complex *X,int *k,int m,int p,int q,complex W);
void newbuffer(complex *X,int* j,int m,int p,int q,complex W,int zhishu);//蝶形计算
void for_dit_fft(int zhishu,complex* X);
void for_dif_fft(int zhishu,complex* X);
complex* gui_dit_fft(int zhishu,complex* X);
complex* gui_dif_fft(int zhishu,complex* X);
void DFT(int zhishu,complex*X);
void main()
{///////////////////////////////2的8次方
	int zhishu;
	for(zhishu=8;zhishu<13;zhishu++)
	{
	int L=int(pow(double(2),double(zhishu)));
	complex* X=new complex[L];
	for(int i=0;i<L;i++)
	{
		X[i].x=1;
		X[i].y=0;
	}
	long beginTime=clock();//获得开始时间，单位为毫秒
	for_dit_fft(zhishu,X);//结果保存于文件中
	long endTime=clock();//获得结束时间
	cout<<"循环方式的时间抽取的FFT计算2的"<<zhishu<<"次方的点所用的时间为"<<endTime-beginTime<<"ms"<<endl;
	for(int i=0;i<L;i++)
	{
		X[i].x=1;
		X[i].y=0;
	}
	beginTime=clock();
	for_dif_fft(zhishu,X);
	endTime=clock();
	cout<<"循环方式的频率抽取的FFT计算2的"<<zhishu<<"次方的点所用的时间为"<<endTime-beginTime<<"ms"<<endl;
	for(int i=0;i<L;i++)
	{
		X[i].x=1;
		X[i].y=0;
	}
	beginTime=clock();
	X=gui_dif_fft(zhishu,X);
	ofstream gui_dif;
	gui_dif.open("gui_dif.txt");
	for(int i=0;i<L;i++)
	{
		gui_dif<<X[i].x<<"+"<<X[i].y<<"i"<<endl;
	}
	endTime=clock();
	cout<<"递归方式的频率抽取的FFT计算2的"<<zhishu<<"次方的点所用的时间为"<<endTime-beginTime<<"ms"<<endl;
	for(int i=0;i<L;i++)
	{
		X[i].x=1;
		X[i].y=0;
	}
	beginTime=clock();
	X=gui_dit_fft(zhishu,X);
	ofstream gui_dit;
	gui_dit.open("gui_dit.txt");
	for(int i=0;i<L;i++)
	{
		gui_dit<<X[i].x<<"+"<<X[i].y<<"i"<<endl;
	}
	endTime=clock();
	cout<<"递归方式的时间抽取的FFT计算2的"<<zhishu<<"次方的点所用的时间为"<<endTime-beginTime<<"ms"<<endl;
	for(int i=0;i<L;i++)
	{
		X[i].x=1;
		X[i].y=0;
	}
	beginTime=clock();
	DFT(zhishu,X);
	endTime=clock();
	cout<<"用DFT计算2的"<<zhishu<<"次方的点所用的时间为"<<endTime-beginTime<<"ms"<<endl;
	cout<<"*************************************************************"<<endl;
	}
}
void for_dit_fft(int zhishu,complex* X)
{
	int j,m,diexing;
	complex W;
	double pi=3.141592654;
	int *k,L=int(pow(double(2),double(zhishu)));
	k=new int[L];
	for(j=0;j<L;j++)
		k[j]=0;
	for(m=0;m<zhishu;m++)//进行zhishu次循环δβχχχχη☆◎※§
	{
		for(diexing=0;diexing<int(pow(double (2),double(zhishu)));)//求解第m次循环所需进行的蝶形计算的次数
		{
			for(j=diexing;j<diexing+int(pow(double(2),double(zhishu-1-m)));j++)//对每个蝶形分别计算
			{
				W.x=cos(2*pi*double((k[j+int(pow(double(2),double(zhishu-1-m)))]%int(pow(double(2),m))))/pow(double(2),m+1));
				W.y=-sin(2*pi*double((k[j+int(pow(double(2),double(zhishu-1-m)))]%int(pow(double(2),m))))/pow(double(2),m+1));
				buffer(X,k,m,j,j+int(pow(double(2),double(zhishu-1-m))),W);//只需再计算出旋转因子即可
			}
		    diexing=diexing+int(pow(double(2),double(zhishu-m)));//注意蝶形要加上偏移量
		}
	}
	ofstream  for_dit;
	for_dit.open("for_dit.txt");
	for( j=0;j<8;j++)
	{
		for_dit<<X[j].x<<"+"<<X[j].y<<"i"<<endl;
	}

}
void buffer(complex *X,int* k,int m,int p,int q,complex W)//蝶形计算
{
	complex A;
	X[q]=cmul(X[q],W);
	A=csub(X[p],X[q]);k[q]=k[q]+int(pow(double(2),double(m)));
	X[p]=cadd(X[p],X[q]);k[p]=k[p]+0;
	X[q]=A;
}

void for_dif_fft(int zhishu,complex* X)
{
	int j,m,diexing;
	complex W;
	double pi=3.141592654;
	int *k,L=int(pow(double(2),double(zhishu)));
	k=new int[L];
	for(j=0;j<L;j++)
		k[j]=j;
	for(m=0;m<zhishu;m++)//进行zhishu次循环δβχχχχη☆◎※§
	{
		for(diexing=0;diexing<int(pow(double (2),double(zhishu)));)//求解第m次循环所需进行的蝶形计算的次数
		{
			for(j=diexing;j<diexing+int(pow(double(2),double(zhishu-1-m)));j++)//对每个蝶形分别计算
			{
				W.x=cos(2*pi*double((k[j+int(pow(double(2),double(zhishu-1-m)))]%int(pow(double(2),zhishu-m-1))))/pow(double(2),zhishu-m));
				W.y=-sin(2*pi*double((k[j+int(pow(double(2),double(zhishu-1-m)))]%int(pow(double(2),zhishu-m-1))))/pow(double(2),zhishu-m));
				newbuffer(X,k,m,j,j+int(pow(double(2),double(zhishu-1-m))),W,zhishu);//只需再计算出旋转因子即可
			}
		    diexing=diexing+int(pow(double(2),double(zhishu-m)));//注意蝶形要加上偏移量
		}
	}
	ofstream  for_dif;
	for_dif.open("for_dif.txt");
	for( j=0;j<8;j++)
	{
		for_dif<<X[j].x<<"+"<<X[j].y<<"i"<<endl;
	}
}
void newbuffer(complex *X,int* j,int m,int p,int q,complex W,int zhishu)//蝶形计算
{
	complex A;
	A=cadd(X[p],X[q]);
	X[q]=csub(X[p],X[q]);j[q]=j[q]-int(pow(double(2),double(zhishu-m-1)));
	X[q]=cmul(X[q],W);
	X[p]=A;
}
complex* gui_dit_fft(int zhishu,complex* X)
{
	double pi=3.141592654;
	int L=int(pow(double(2),double(zhishu)));
	if(L==2)
	{
		complex a=X[0];
		complex b=X[1];
		X[0]=cadd(a,b);
		X[1]=csub(a,b);
		return X;
	}
	else
	{
		complex* Y;
		complex* Z;
		complex W;
		Y=new complex[L/2];
		Z=new complex[L/2];
		for(int i=0;i<L/2;i++)
		{
			Y[i]=X[2*i];
			Z[i]=X[2*i+1];
		}
		Y=gui_dit_fft(zhishu-1,Y);
		Z=gui_dit_fft(zhishu-1,Z);
		for(int k=0;k<L/2;k++)
		{
			W.x=cos(2*pi*double(k)/double(L));
			W.y=-sin(2*pi*double(k)/double(L));
			X[k]=cadd(Y[k],cmul(Z[k],W));
			X[k+L/2]=csub(Y[k],cmul(Z[k],W));
		}
		return X;

	}
}
complex* gui_dif_fft(int zhishu,complex* X)
{
	double pi=3.141592654;
	int L=int(pow(double(2),double(zhishu)));
	if(L==2)
	{
		complex a=X[0];
		complex b=X[1];
		X[0]=cadd(a,b);
		X[1]=csub(a,b);
		return X;
	}
	else
	{
		complex* Y;
		complex* Z;
		complex W;
		Y=new complex[L/2];
		Z=new complex[L/2];
		for(int n=0;n<L/2;n++)
		{
			W.x=cos(2*pi*double(n)/double(L));
			W.y=-sin(2*pi*double(n )/double(L));
			Y[n]=cadd(X[n],X[n+L/2]);
			Z[n]=cmul(csub(X[n],X[n+L/2]),W);
		}
		Y=gui_dif_fft(zhishu-1,Y);
		Z=gui_dif_fft(zhishu-1,Z);
		for(int k=0;k<L/2;k++)
		{
			X[2*k]=Y[k];
			X[2*k+1]=Z[k];
		}
		return X;
	}
}


void DFT(int zhishu,complex*X)
{
    int L=int(pow(double(2),double(zhishu)));
    double pi=3.141592654;
	complex* Y=new complex[L];
	complex W;
	ofstream DFT;
	DFT.open("DFT.txt");
	for(int i=0;i<L;i++)
	{
		Y[i].x=0;
		Y[i].y=0;
		for(int j=0;j<L;j++)
		{
			W.x=cos(2*pi*double(i)*double(j)/double(L));
			W.y=-sin(2*pi*double(i)*double(j)/double(L));
			Y[i]=cadd(Y[i],cmul(X[j],W));
			if (abs(Y[i].x)<1e-5)
				Y[i].x=0;
			if (abs(Y[i].y)<1e-5)
				Y[i].y=0;
		}
		DFT<<Y[i].x<<"+"<<Y[i].y<<"i"<<endl;
	}
}