FFT.rar_fft c_fft c code image_fft c语言_image FFT_傅里叶

//验证正确 #include "math.h" #include "malloc.h" #include "conio.h" #include "stdio.h" #include "string.h" #define pi (double)3.14159265359 #define m 4 //复数定义 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 Mul(COMPLEX c1,COMPLEX c2) { COMPLEX c; c.re=c1.re*c2.re-c1.im*c2.im; c.im=c1.re*c2.im+c2.re*c1.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; } //快速傅里叶变换 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); 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); } //逆变换 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); } void main() { int count; //YV为进行逆傅里叶变换后的时域值 COMPLEX *TV,*FV,*YV; int i; count=1<<m; TV=(COMPLEX *)malloc(sizeof(COMPLEX)*count); FV=(COMPLEX *)malloc(sizeof(COMPLEX)*count); YV=(COMPLEX *)malloc(sizeof(COMPLEX)*count); for (i=0;i<count;i++) { TV[i].re=i*3; TV[i].im=0; } TV[0].re=0.5751; TV[1].re=0.4514; TV[2].re=0.0439; TV[3].re=0.0272; TV[4].re=0.3127; TV[5].re=0.0129; TV[6].re=0.3840; TV[7].re=0.6831; TV[8].re=0.0928; TV[9].re=0.0353; TV[10].re=0.6124; TV[11].re=0.6085; TV[12].re=0.0158; TV[13].re=0.0164; TV[14].re=0.1901; TV[15].re=0.5869; FFT(TV,FV,m); IFFT(FV,YV,m); printf("傅里叶后的频域值为：\n"); for (i=0;i<count;i++) { printf("x[%d]=%f+(%f)i\n",i,FV[i].re,FV[i].im); } printf("\n"); printf("逆傅里叶后的时域值为：\n"); for (i=0;i<count;i++) { printf("Y[%d]=%f+(%f)i\n",i,YV[i].re,YV[i].im); } free(TV); free(FV); free(YV); // getch(); } /* //////////////////////////////////////////////// // sample: input.txt //////////////////////////////////////////////// 0.5751 0 0.4514 0 0.0439 0 0.0272 0 0.3127 0 0.0129 0 0.3840 0 0.6831 0 0.0928 0 0.0353 0 0.6124 0 0.6085 0 0.0158 0 0.0164 0 0.1901 0 0.5869 0 //////////////////////////////////////////////// // sample: output.txt //////////////////////////////////////////////// FFT: i real imag 0 4.6485 0.0000 1 0.0176 0.3122 2 1.1114 0.0429 3 1.6776 -0.1353 4 -0.2340 1.3897 5 0.3652 -1.2589 6 -0.4325 0.2073 7 -0.1312 0.3763 8 -0.1949 0.0000 9 -0.1312 -0.3763 10 -0.4326 -0.2073 11 0.3652 1.2589 12 -0.2340 -1.3897 13 1.6776 0.1353 14 1.1113 -0.0429 15 0.0176 -0.3122 ================================= IFFT: i real imag 0 0.5751 -0.0000 1 0.4514 0.0000 2 0.0439 -0.0000 3 0.0272 -0.0000 4 0.3127 -0.0000 5 0.0129 -0.0000 6 0.3840 -0.0000 7 0.6831 0.0000 8 0.0928 0.0000 9 0.0353 -0.0000 10 0.6124 0.0000 11 0.6085 0.0000 12 0.0158 0.0000 13 0.0164 0.0000 14 0.1901 0.0000 15 0.5869 -0.0000 */