//========================================
//函数名:_FFT
//功能:FFT变换
//用法:void _FFT(double *fr, double *fi, int n, int flag)
//参数:fr:采样点的实数表,fi:采样点的虚数表,n:采样点个数,
// flag:flag=0表示求Fourier变换,flag=1表示求逆Fourier变换
//返回:无
//========================================
void _FFT(double *fr, double *fi, int n, int flag)
{
int mp,arg,q,cntr,p1,p2;
int i,j,a,b,k,m;
double sign,pr,pi,harm,x,y,t;
double *ca,*sa;
ca=(double *)calloc(n,sizeof(double));
sa=(double *)calloc(n,sizeof(double));
if(ca==NULL||sa==NULL) Error("calloc error in _FFT!");
j=0;
if(flag!=0)
{
sign=1.0;
for(i=0;i<=n-1;i++)
{
fr[i]=fr[i]/n;
fi[i]=fi[i]/n;
}
}
else
sign=-1.0;
for(i=0;i<=n-2;i++)
{
if(i<j)
{
t=fr[i];
fr[i]=fr[j];
fr[j]=t;
// t=fi[i];
// fi[i]=fi[j];
// fi[j]=t;
}
k = n / 2;
while(k <= j)
{
j -= k;
k /= 2;
}
j+=k;
}
mp=0;
i=n;
while(i!=1)
{
mp+=1;
i/=2;
}
harm=2*M_PI/n;
for(i=0;i<=n-1;i++)
{
sa[i]=sign*sin(harm*i);
ca[i]=cos(harm*i);
}
a=2;
b=1;
for(cntr=1;cntr<=mp;cntr++)
{
p1=n/a;
p2=0;
for(k=0;k<=b-1;k++)
{
i=k;
while(i<n)
{
arg=i+b;
if(k==0)
{
pr=fr[arg];
pi=fi[arg];
}
else
{
pr=fr[arg]*ca[p2]-fi[arg]*sa[p2];
pi=fr[arg]*sa[p2]+fi[arg]*ca[p2];
}
fr[arg]=fr[i]-pr;
fi[arg]=fi[i]-pi;
fr[i]+=pr;
fi[i]+=pi;
i+=a;
}
p2+=p1;
}
a*=2;
b*=2;
}
free(ca);
free(sa);
}