void fft(double Xr[],double Xi[],int n,int inv)
{
int i=0,j=0; //n为fft的点数,m为2的次方数
int m=1,nn=n;
int I=0,J=0,k=0;
int L=0,B=0,P=0;
double Tr=0,Ti=0,Rr=0,Ii=0;
while(nn/2!=1){m=m+1;nn=nn/2;} //求m的值
for(i=0;i<n;i++)Xi[i]=inv*Xi[i]; //如果inv=-1,取共轭
Inverse(Xr,n); //调用倒序子函数Inverse,倒序
Inverse(Xi,n);
for(L=1;L<=m;L++) //DIF-FFT Decimation-In-Time 基2-FFT
{
B=(int)pow(2.0,double(L-1));
for(J=0;J<=B-1;J++)
{
k=0;
P=(int)(J*pow(2.0,double(m-L))); //计算旋转因子
for(k=J;k<n-1;)
{
Tr=Xr[k+B]*cos(2*pi*P/n)+Xi[k+B]*sin(2*pi*P/n);
Ti=Xi[k+B]*cos(2*pi*P/n)-Xr[k+B]*sin(2*pi*P/n);
Rr=Xr[k]+Tr;Ii=Xi[k]+Ti; //Xr[k]为实部,Xi[k]为虚部
Xr[k+B]=Xr[k]-Tr;Xi[k+B]=Xi[k]-Ti;
Xr[k]=Rr;Xi[k]=Ii;
k=k+(int)pow(2.0,double(L));
}
}
}
for(i=0;i<n;i++)Xi[i]=inv*Xi[i]; //如果inv=-1,取共轭
if(inv==-1) //如果inv=-1,要在X[N]前乘以1/N
for(i=0;i<n;i++)
{
Xr[i]=(1.0/n)*Xr[i];
Xi[i]=(1.0/n)*Xi[i];
}
}
void Inverse(double X[],int n) //倒序子程序
{
int I=0,J=0,LH=0,N1=0,K=0;
double T=0;
LH=n/2,N1=n-2;J=LH; //初始化
for(I=1;I<=N1;I++)
{
if(I<J){T=X[I];X[I]=X[J];X[J]=T;}
K=LH;
if(J<K)J=J+K;
else
{
while(J>=K)
{
J=J-K;
K=K/2;
}
J=J+K;
}
}
}
