fft.zip_FFT 蝶形

#include<stdio.h> #include<math.h> #include<iostream.h> #define PI 3.1415926 struct compx {float real,imag;}; struct compx mul(compx a,compx b) { struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return c; } struct compx add(compx a,compx b) { struct compx c; c.real=a.real+b.real; c.imag=a.imag+b.imag; return c; } struct compx sub(compx a,compx b) { struct compx c; c.real=a.real-b.real; c.imag=a.imag-b.imag; return c; } struct compx div(compx a,compx b) { struct compx c; c.real=(a.real*b.real+a.imag*b.imag)/(b.real*b.real+b.imag*b.imag); c.imag=(a.imag*b.real-a.real*b.imag)/(b.real*b.real+b.imag*b.imag); return c; } void main() { int M,i,j,k,l; printf("please input M:\n"); scanf("%d",&M); int N=pow(2,M); int MM=M+1,NN=N+1; int **a=new int*[NN]; for(i=0;i<NN;i++) { a[i]=new int[MM]; } for(i=0;i<NN;i++) { for(j=0;j<MM;j++) { a[i][j]=0; } } int *b=new int[MM]; for (i=1;i<MM;i++) { b[i]=pow(2,M-i); } for (j=1;j<=M;j++) { int te=pow(2,j-1); for (k=1;k<=N/te;k++) { if (k%2==0) { for (l=(k-1)*te+1;l<=k*te;l++) { a[l][j]=1; } } } } int *y=new int[NN]; for (i=1;i<=N;i++) { y[i]=0; for (j=1;j<=M;j++) { y[i]+=a[i][j]*b[j]; } }////完成倒序 printf("倒序下标顺序为：\n"); for (i=1;i<=N;i++) { printf("%d\n",y[i]); } printf("********************\n"); int *c=new int[N]; printf("please input datas(数据数量需要为2的M次幂):\n"); for(i=0;i<N;i++) { scanf("%d",&c[i]); } compx *cc=new compx[N]; for (i=0;i<N;i++) { cc[i].real=c[y[i+1]]; cc[i].imag=0; printf("%f\n",cc[i].real); }//将输入数组倒序 int s=pow(2,M-1); compx *jj=new compx[s+1]; struct compx WN1; WN1.real=cos(2*PI/N); WN1.imag=0-sin(2*PI/N); struct compx ct=mul(WN1,WN1); struct compx q; q=ct; jj[0].real=1; jj[0].imag=0; jj[1].real=WN1.real; jj[1].imag=WN1.imag; jj[2].real=ct.real; jj[2].imag=ct.imag; for (int pp=3;pp<s;pp++) { jj[pp]=mul(q,WN1); q=jj[pp]; } printf("***************************\n"); compx t; for(i=1;i<=M;i++) { int B=pow(2,i-1); for(j=0;j<=B-1;j++) { int p=j*pow(2,M-i); for(k=j;k<=N-1;k=(k+pow(2,i))) { t=cc[k]; cc[k]=add(t,mul(cc[k+B],jj[p])); cc[k+B]=sub(t,mul(cc[k+B],jj[p])); } } } for (i=0;i<N;i++) { printf("%f +%fi \n",cc[i].real,cc[i].imag); } for(i=0;i<NN;i++) delete[] a[i]; delete[] a; delete[] b; delete[] y; delete[] cc; delete[] c; delete[] jj; }