傅里叶变换

#include <iostream>
#include <complex>
#include <math.h>
#include <ctime>
#define pi 3.14159265358979323846
#define NUM 1000
using namespace std;

void dft(complex<double>*Data,int length,int flag)
{
    int i,j;
    complex<double> wn;
    complex<double>*temp=new complex<double>[length];
    for(i=0;i<length;i++)
    {
		temp[i]=0;
		for(j=0;j<length;j++)
		{
			wn=complex<double>(cos(2.0*pi/length*i*j),sin(flag*2.0*pi/length*i*j));
			temp[i]+=Data[j]*wn;
		}
    }
	
    for(i=0;i<length;i++)
	{
		Data[i]=temp[i];
		if(flag==1)Data[i]/=length;
	}
	delete[] temp;
}
void reverse2(complex<double> *data,int Log2N)
{
	int i,j;
	int RevNum;
	int MaxPos,CurPos,MaxValue;
	complex<double> temp;
	MaxValue=(1<<Log2N)-1;
	MaxPos=1<<(Log2N-1);
	RevNum=0;
	
	for(i=1;i<MaxValue;i++)
	{
		
		CurPos=MaxPos;
		while((CurPos&RevNum)!=0)
		{
			RevNum=RevNum&(~CurPos);
			CurPos=CurPos>>1;
		}
		RevNum=RevNum|CurPos;
		if (RevNum<i)
		{
			temp=data[RevNum];
			data[RevNum]=data[i];
			data[i]=temp;
		}
	}
}

complex<double>* dit2(complex<double>*Data,int Log2N,int flag)
{
	int i,j,k;
	int length=1<<Log2N;
	complex<double> wn,temp;
	int step;
	int index0,index1;
	
	for(i=1;i<=Log2N;i++)
	{
		step=1<<i;
		
		for(j=0;j<step/2;j++)
		{
			wn=complex<double>(cos(2*pi/step*j),sin(flag*2*pi/step*j));
			
			for(k=0;k<length/step;k++)
			{
				index0=k*step+j;
				index1=k*step+step/2+j;
				temp=Data[index0];
				Data[index0]=temp+wn*Data[index1];
				Data[index1]=temp-wn*Data[index1];
			}
		}
		
	}
	
	if(flag==1)//如果为反变换，要除以序列的长度
		for(int i=0;i<length;i++)
			Data[i]/=length;
		return Data;
		
}
int main()
{
	complex<double> x[NUM],y[NUM],ix[NUM],ox[NUM],x2[NUM];
	double d[NUM];
	time_t beg,end;
	int i;
	for(i=0;i<NUM;++i)
	{
		d[i]=sin(0.1*i);
		x[i].real(sin(0.1*i));
		x[i].imag(0);
		y[i].real(sin(0.1*i));
		y[i].imag(0);
		ix[i].real(sin(0.1*i));
		ix[i].imag(0);
		x2[i].real(sin(0.1*i));
		x2[i].imag(0);
	}
	beg=clock();
	reverse2(x,6);
	end=clock();
	cout<<end-beg<<endl;
	
	beg=clock();
	dit2(x,6, -1);
	end=clock();
	cout<<end-beg<<endl;
	
	beg=clock();
	reverse2(x,6);
	end=clock();
	cout<<end-beg<<endl;
	
	beg=clock();
	dit2(x,6, 1);
	end=clock();
	cout<<end-beg<<endl;
	beg=clock();
	dft(y,NUM,-1);
	end=clock();
	cout<<end-beg<<endl;

	beg=clock();
	dft(y,NUM,1);
	end=clock();
	cout<<end-beg<<endl;
	
	return 0;
	for (i = 0; i < NUM; i++)
	{
		printf("%.2lf %.2lf %.2lf %.2lf\n",y[i].real(),x[i].real(),d[i]);
	}
	return 0;
}