#include "config.h"
#include<math.h>
struct fushu //fushu结构体表示复数
{ float re; //实部
float im; //虚部
};
void fsmul(struct fushu *p,struct fushu a,struct fushu b) //复数乘法
{ p->re=a.re*b.re-a.im*b.im;
p->im=a.re*b.im+a.im*b.re;
}
void fsadd(struct fushu * p,struct fushu a,struct fushu b) //复数加法
{ p->re=a.re+b.re;
p->im=a.im+b.im;
}
void fssub(struct fushu * p,struct fushu a,struct fushu b) //复数减法
{ p->re=a.re-b.re;
p->im=a.im-b.im;
}
/*****************************************/
/* DIT的基－2FFT算法 */
/* a: 待变换序列指针 */
/* l: 待变换序列长度 */
/* n: FFT点数 */
/*****************************************/
void FFT(struct fushu *a,uint16 l,uint16 n)
{
uint16 i,j,nv2,nm1,k,le,lei,ip,m,L;
struct fushu u,w,t;
float tmp;
float pai=atan(1)*4;
if(l<n) //若序列长度不够FFT变换点数则补零
for(i=l;i<n;i++)
{a[i].re=0;a[i].im=0;}
L=log(n)/log(2);
nv2=n/2;
nm1=n-1;
j=0;
for(i=0;i<nm1;i++) //将顺序输入的序列变换成倒位序
{if(i<j)
{t.re=a[j].re;t.im=a[j].im;
a[j].re=a[i].re;a[j].im=a[i].im;
a[i].re=t.re;a[i].im=t.im;
}
k=nv2;
while(k<=j)
{j=j-k;k=k/2;}
j+=k;
}
le=1; //进行FFT运算
for(m=1;m<=L;m++)
{lei=le;
le=le*2;
u.re=1;
u.im=0;
tmp=pai/lei;
w.re=cos(tmp);
w.im=-sin(tmp);
for(j=0;j<lei;j++)
{for(i=j;i<n;i+=le)
{ip=i+lei;
fsmul(&t,a[ip],u);
fssub(&(a[ip]),a[i],t);
fsadd(&(a[i]),a[i],t);
}
fsmul(&u,u,w);
}
}
}
