急求 快速傅里叶变换的程序 会的大侠帮忙做一个~

以前用FFT都是直接用Matlab里面的，可是如果实际工程里面需要，还是得写一个C语言版本的。C++处理复数比较容易，但目前嵌入式开发还是C语言的天下，因此C语言的FFT应用起来更方便。写完贴出来，希望对大家有一些帮助。贴出来排版可能有点乱，那不是我的原因，我写的程序都是整整齐齐的，可以直接点击文章后面的目录下载源程序。
最近C程序写的比较多，C++好久不写，有点荒废了。。。
/**
* FFT - Fast Fourier transform. The length of X must be a power
* of two, for a fast radix-2 fast-Fourier transform algorithm
* is used. spadger@bmy 2007.9.2
*/

＃i nclude
＃i nclude

#define M_PI 3.14159265358979323846

typedef struct { double r,i; } cplx_t;

void cplx_mul(cplx_t *x, cplx_t *y, cplx_t *r)
{
r->r=x->r*y->r-x->i*y->i;
r->i=x->r*y->i+x->i*y->r;
}

void cplx_exp(cplx_t *x, cplx_t *r)
{
double expx=exp(x->r);
r->r=expx*cos(x->i);
r->i=expx*sin(x->i);
}

void bit_reverse(cplx_t *x, int N)
{
double t;
cplx_t tmp;
unsigned int i=0,j=0,k=0;
for(i=0; ik=i;
j=0;
t=log(0.0+N)/log(2.0);
while((t--)>0) {
j<<=1;
j|=k&1;
k>>=1;
}
if(j>i) {
tmp=x[i];
x[i]=x[j];
x[j]=tmp;
}
}
}

void fft(cplx_t *x, int N)
{
cplx_t u,d,p,W,tmp;
int i=0,j=0,k=0,l=0,M=floor(log(0.0+N)/log(2.0));
if(log(0.0+N)/log(2.0)-M > 0){
printf("The length of x (N) must be a power of two!!!\n");
return;
}
bit_reverse(x,N);
for(i=0; il=1<for(j=0; jfor(k=0; ktmp.r=0.0;
tmp.i=-2*M_PI*k/2/l;
cplx_exp(&tmp,&W);
cplx_mul(&x[j+k+l],&W,&p);
u.r=x[j+k].r+p.r;
u.i=x[j+k].i+p.i;
d.r=x[j+k].r-p.r;
d.i=x[j+k].i-p.i;
x[j+k]=u;
x[j+k+l]=d;
}
}
}
}

/**
* for test and demonstation, set '#if 0' to comment this out.
*/
#if 1
#define DATA_LEN 64
int main()
{
int i;
cplx_t x[DATA_LEN];
for(i=0;ix[i].r=i;
x[i].i=0;
}

printf("Before...\nReal\t\tImag\n");
for(i=0;iprintf("%f\t%f\n",x[i].r,x[i].i);

fft(x,DATA_LEN);

printf("After...\nReal\t\tImag\n");
for(i=0;iprintf("%f\t%f\n",x[i].r,x[i].i);

return 0;
}
#endif
源程序下载：
http://www.blog.com.cn/user88/spadger/upload/92271850.rar