#include

#include

#define pi (double) 3.14159265359

/*复数的定义*/

typedef struct

{

double re;

double im;

}COMPLEX;

/*复数的加运算*/

COMPLEX Add(COMPLEX c1, COMPLEX c2)

{

COMPLEX c;

c.re = c1.re + c2.re;

c.im = c1.im + 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 Mul(COMPLEX c1, COMPLEX c2)

{

COMPLEX c;

c.re = c1.re*c2.re - c1.im*c2.im;

c.im = c1.re*c2.im + c1.im*c2.re;

return c;

}

/*快速傅立叶变换

TD为时域值，FD为频域值，power为2的幂数*/

void FFT(COMPLEX *TD, COMPLEX *FD, int power)

{

int count;

int i,j,k,bfsize,p;

double angle;

COMPLEX *W,*X1,*X2,*X;

/*计算傅立叶变换点数*/

count=1<

/*分配运算器所需存储器*/

W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);

X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);

X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);

/*计算加权系数*/

for(i=0;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

{

for(j=0;j<1<

{

bfsize=1<

for(i=0;i

{

p=j*bfsize;

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<

}

}

X=X1;

X1=X2;

X2=X;

}

/*重新排序*/

for(j=0;j

{

p=0;

for(i=0;i

{

if(j&(1<

p+=1<

}

FD[j]=X1[p];

}

/*释放存储器*/

free(W);

free(X1);

free(X2);

}

/*快速傅立叶反变换，利用快速傅立叶变换

FD为频域值，TD为时域值，power为2的幂数*/

void IFFT(COMPLEX *FD, COMPLEX *TD, int power)

{

int i,count;

COMPLEX *x;

/*计算傅立叶反变换点数*/

count=1<

/*分配运算所需存储器*/

x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);

/*将频域点写入存储器*/

memcpy(x,FD,sizeof(COMPLEX)*count);

/*求频域点的共轭*/

for(i=0;i

{

x[i].im=-x[i].im;

}

/*调用快速傅立叶变换*/

FFT(x,TD,power);

/*求时域点的共轭*/

for(i=0;i

{

TD[i].re/=count;

TD[i].im=-TD[i].im/count;

}

/*释放存储器*/

free(x);

}