傅立叶变换与傅立叶反变换的C语言实现

#include <math.h>

#include <malloc.h>

#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<<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++)

{

bfsize=1<<power-k;

for(i=0;i<bfsize/2;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<<k)]);

}

}

X=X1;

X1=X2;

X2=X;

}

/*重新排序*/

for(j=0;j<count;j++)

{

p=0;

for(i=0;i<power;i++)

{

if(j&(1<<i))

p+=1<<power-i-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<<power;

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

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

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

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

/*求频域点的共轭*/

for(i=0;i<count;i++)

{

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

}

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

FFT(x,TD,power);

/*求时域点的共轭*/

for(i=0;i<count;i++)

{

TD[i].re/=count;

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

}

/*释放存储器*/

free(x);

}