#include "math.h"
typedef struct
{
double real;
double imag;
} compx ;//复数定义
compx compx_mul(compx b1,compx b2)
{
compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}//复数乘法
void FFT(compx *xin,int N)
{
int f,m,LH,nm,i,k,j,L;
double p , ps ;
int le,B,ip;
float pi;
compx w,t;
LH=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++)//f为总点数，m为级数
nm=N-2;
j=LH;
/*变址运算, 倒序*/
for(i=1;i<=nm;i++)
{
if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}
for(L=1;L<=m;L++)//第一重级数循环
{
le=(int)pow(2,L);
B=le/2;
pi=(float)3.14159;
for(j=0;j<=B-1;j++)//第二重循环计算旋转因子
{
p=pow(2,m-L)*j;
ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);
for(i=j;i<=N-1;i=i+le)// 第三重循环计算下一级的序列
{
ip=i+B;
t=compx_mul(xin[ip],w);
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
}
}
}
void IFFT(compx *X,int N)
{
int i;
for(i=0;i<Num;i++)X[i].imag=-X[i].imag;
FFT(X,N);
for(i=0;i<Num;i++)
{
X[i].imag=-X[i].imag/Num;
X[i].real=X[i].real/Num;
}
}
评论0