FFT的C語言算法實(shí)現(xiàn)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上FFT的C語言算法實(shí)現(xiàn)程序如下： /*FFT*/ #include <stdio.h> #include <math.h> #include <stdlib.h> #define N 1000 typedef struct double real; double img; complex; void fft(); /*快速傅里葉變換*/ void ifft(); /*快速傅里葉逆變換*/ void initW(); void change(); void add(complex ,complex ,complex *); /*復(fù)數(shù)加

2、法*/ void mul(complex ,complex ,complex *); /*復(fù)數(shù)乘法*/ void sub(complex ,complex ,complex *); /*復(fù)數(shù)減法*/ void divi(complex ,complex ,complex *);/*復(fù)數(shù)除法*/ void output(); /*輸出結(jié)果*/ complex xN, *W;/*輸出序列的值*/ int size_x=0;/*輸入序列的長度，只限2的N次方*/ double PI; int main() int i,method; system("cls"); PI=atan(

3、1)*4; printf("Please input the size of x:n"); /*輸入序列的長度*/ scanf("%d",&size_x); printf("Please input the data in xN:(such as:5 6)n"); /*輸入序列對(duì)應(yīng)的值*/ for(i=0;i<size_x;i+) scanf("%lf %lf",&xi.real,&xi.img); initW(); /*選擇FFT或逆FFT運(yùn)算*/ printf("Use

4、FFT(0) or IFFT(1)?n"); scanf("%d",&method); if(method=0) fft(); else ifft(); output(); return 0; /*進(jìn)行基-2 FFT運(yùn)算*/ void fft() int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i< log(size_x)/log(2) ;i+) l=1<<i; for(j=0;j<size_x;j+= 2*l ) for(k=0;k<l;k+)

5、mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; void ifft() int i=0,j=0,k=0,l=size_x; complex up,down; for(i=0;i< (int)( log(size_x)/log(2) );i+) /*蝶形運(yùn)算*/ l/=2; for(j=0;j<size_x;j+= 2*l ) for(k=0;k<l;k+) add(xj+k,xj+k+l,&

6、amp;up); up.real/=2;up.img/=2; sub(xj+k,xj+k+l,&down); down.real/=2;down.img/=2; divi(down,Wsize_x*k/2/l,&down); xj+k=up; xj+k+l=down; change(); void initW() int i; W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;i<size_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.img=-1*sin(2*PI/size_x*i)

7、; void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i<size_x;i+) k=i;j=0; t=(log(size_x)/log(2); while( (t-)>0 ) j=j<<1; j|=(k & 1); k=k>>1; if(j>i) temp=xi; xi=xj; xj=temp; void output() /*輸出結(jié)果*/ int i; printf("The result are as followsn");

8、 for(i=0;i<size_x;i+) printf("%.4f",xi.real); if(xi.img>=0.0001) printf("+%.4fjn",xi.img); else if(fabs(xi.img)<0.0001) printf("n"); else printf("%.4fjn",xi.img); void add(complex a,complex b,complex *c) c->real=a.real+b.real; c->img=a.img+b.img; void mul(complex a,complex b,complex *c) c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; void sub(complex a,complex b,complex *c) c->real=a.real-b.real; c->img=a.img-b.img; void divi(complex a,complex b,complex *c) c->r