#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);
}