[code=C/C++][/code]
typedef struct
{
float re;
float im;
}COMPLEX;
//复数乘法
COMPLEX Mul(COMPLEX c1,COMPLEX c2)
{
COMPLEX c;
c.re=c1.re*c2.re-c1.im*c2.im;
c.im=c1.im*c2.re+c1.re*c2.im;
return c;
}
//复数减法
COMPLEX Sub(COMPLEX c1,COMPLEX c2)
{
COMPLEX c;
c.re=c1.re-c2.re;
c.im=c1.im-c2.im;
return c;
}
//复数加法
COMPLEX Add(COMPLEX c1,COMPLEX c2)
{
COMPLEX c;
c.re=c1.re+c2.re;
c.im=c1.im+c2.im;
return c;
}
/************快速傅立叶变换**************************
*****************************************************/
//TD为时域值，FD为频域值，power为变换点数对2取的对数
void FFT(COMPLEX *TD,COMPLEX *FD,int power)
{
/*******************定义数据*********************/
int count;//时域或频域变换的点数
int i,j,k,bfsize,p;
//k是叠代次数，j第几个分组，bfsize是蝶距，p是第j个
//分组在整个数据点数组中的相对位置
double angle; //相角
COMPLEX *W,*X1,*X2,*X;//存放数据点的空间地址
//W是存放加权系数的地址
//从X1到X2是进行了一次叠代，X1是叠代前数据点的地址，
//X2是叠代后的数据点的地址,X是临时存放数据点的地址
/******************计算fft点数*******************/
count=1<<power;
/*****************分配运算所需的存储器***********/
W=(COMPLEX*)malloc(sizeof(COMPLEX)*count/2);
X1=(COMPLEX*)malloc(sizeof(COMPLEX)*count);
X2=(COMPLEX*)malloc(sizeof(COMPLEX)*count);
/*********************计算加权系数****************/
for(i=0;i<count/2;i++)
{
angle=-i*pi*2/count;
W[i].re=cos(angle);
W[i].im=sin(angle);
}
/*********************写入时域点******************/
memcpy(X1,TD,sizeof(COMPLEX)*count);
/************计算蝶形*****************************/
for(k=0;k<power;k++) //最外的循环是叠代次数的循环
{
for(j=0;j<1<<k;j++) //在某一次叠代中，对j的循环是遍历整个数组的每一个分组
//分组指的是：一个N点DFT变为两个N/2点DFT,分组由一组变成了两组。
//对应第k次叠代分组数为2的k次方。
//另外，分组数乘以蝶距就是变换点的总数：2^k*bfsize=count
{
bfsize=1<<(power-k); //蝶距
for(i=0;i<bfsize/2;i++) //i是每一个分组内第几个点
{
p=j*bfsize; //p是每个分组起始点的序号
//蝶形计算
X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
}
}
//为下一次叠代调整存储空间
X=X1;
X1=X2;
X2=X;
}
/*********************重排顺序*******************************/
//j是要反序的数
for(j=0;j<count;j++)
{
p=0;
for(i=0;i<power;i++)//i代表反序的二进制位数
{
if(j&(1<<i)) p+=1<<(power-i-1);