fft.rar_C 傅里叶变换_visual c

#include "conio.h" #include "math.h" #include "stdio.h" #define N 64 #define PI 3.1415926 #define w0 (0.125*PI) #define Cmul(a,b,c) a.x=b.x*c.x-b.y*c.y;a.y=b.x*c.y+b.y*c.x; /*利用结构体做复数乘法*/ #define Cequal(a,b) a.x=b.x;a.y=b.y; /*结构体a的值赋给b*/ #define Cadd(a,b,c) a.x=b.x+c.x;a.y=b.y+c.y; /*结构体实现复数的加减*/ #define Csub(a,b,c) a.x=b.x-c.x;a.y=b.y-c.y; #define Wn(w,r) w.x=cos(2.0*PI*r/n);w.y=-sin(2.0*PI*r/n); /*旋转因子的分解*/ struct comp { float x; float y; }; void main() { int i,j,nu2,nm1,n,m,le,le1,k,ip,z; int flag,f,n1; struct comp a[N],t,t1,w,d; float a_ipx,m1; printf("\nThis program is about FFT by DIF way. "); printf("\nplease enter N : "); scanf("%d",&n1); n=n1; m1=log(n1)/log(2); m=log(n1)/log(2); /*级数 */ if (m!=m1) n=pow(2,m+1); for(i=0;i<n;i++) {a[i].x=a[i].y=0.0;} /*所有存储单元清零*/ printf("\n"); for(i=0;i<n1;i++) /*录入数据到存储单元，以结构体的形式分别存储实部和虚部*/ { printf("\nplease enter data(%d)_[Re]: ",i); scanf("%f",&a[i].x); printf("\nplease enter data(%d)_[Im]: ",i); scanf("%f",&a[i].y); } for(z=0;z<=2;z++) /*z=0,时计算FFT；z=1时计算新序列的IFFT；z=2时计算原序列的IFFT*/ { flag=-1; /*flag控制r的增加*/ for (m=(log(n)/log(2));m>=1;m--) /*此循环控制级数m*/ { le=pow(2,m); /*le控制同一级中具有相同 的碟形结的距离*/ flag++; le1=le/2; /*le1控制在同一级中有不同 的个数*/ for( j=0;j<le1;j++) /*此循环控制同一级中r的增加，实现一级不同的蝶形运算*/ { for (i=j;i<=(n-1);i+=le) /*此循环控制一级中具有相同 的碟形运算*/ { ip=i+le1; Cequal(t,a[i]); Cequal(t1,a[ip]); f=(int) (i*pow(2,flag))%n; Wn(w,f); Cadd(a[i],t,t1); /*计算 */ Csub(a[ip],t,t1); /*计算 a_ipx=a[ip].x; */ if (z==1) { w.y*=-1; } if (z==2) { w.y*=-1; } a[ip].x=a[ip].x*w.x-a[ip].y*w.y; a[ip].y=a_ipx*w.y+a[ip].y*w.x; /*计算 ，反变换则 取反*/ } } } nu2=n/2; /*用雷德算法将倒位序输出变为自然序输出*/ nm1=n-2; j=0;i=0; while(i<=nm1) { if (i<j) { Cequal(d,a[j]); Cequal(a[j],a[i]); Cequal(a[i],d); } k=nu2; while(k<=j) { j=j-k;k=k/2; } j=j+k; i=i+1; } if(z==0) { printf("\n input num fft:\n\n"); } else if(z==1) printf("\n sheng cheng xin num ifft(input):\n\n" ) ; else printf("\n input num ifft"); for(i=0;i<n;i++) if(z==0) { printf(" %7.4f",a[i].x); if (a[i].y>=0) printf(" + %7.4f j \n",a[i].y); else printf(" - %7.4f j \n",fabs(a[i].y)); } else if(z==1) { a[i].x=a[i].x/n; /*IFFT需将结果除以n*/ printf(" %7.4f",a[i].x); a[i].y=a[i].y/n; if (a[i].y>=0) printf(" + %7.4f j \n",a[i].y); else printf(" - %7.4f j \n",fabs(a[i].y)); } else { printf(" %7.4f",a[i].x/n); a[i].y=a[i].y/n; if (a[i].y>=0) printf(" + %7.4f j \n",a[i].y); else printf(" - %7.4f j \n",fabs(a[i].y)); } } printf("\n"); }