#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;
}
傅里叶变换
最新推荐文章于 2024-05-14 15:17:16 发布