FFT算法的C语言编程

#includestdio.h#includemath.h#includestdlib.h#definePI3.1415926/*ThisprogramisusedtoconductFFTorIFFTtransformbasedon2discompositiondatacanbeloadedbytypingonthekeyboardorfromfile.*/typedefstruct{doublereal;doubleimag;}complex;longNN;longt2t(longn,intM){longnn=0;for(inti=0;iM;i++){nn=nn+n%2*pow(2,M-1-i);n=n/2;}returnnn;}intdefM(longn){intM=0;while(pow(2,M)n){M++;}returnM;}voiddatain(complexx[],longN){intflag;longi;chars[20];FILE*fpi;doublea,b;printf(Loaddatafromkeybord(input0)orfile(input1)\n);scanf(%d,&flag);if(flag==0){printf(Pleaseenterthedataintheformofa,b--a+j*bendingwith'*'\n);for(i=0;getchar()!='*';i++)scanf(%lf%lf,&x[i].real,&x[i].imag);while(iN){x[i].real=0;x[i].imag=0;i++;}}elseif(flag==1){printf(PleaseinputthenameofthefilewithitsFiletypesuffix-(eg:data.txt)\n);gets(s);gets(s);//Toavoidthequickresponsetoenterfpi=fopen(s,r);if(!fpi){printf(Openfileerror);}i=0;while(!feof(fpi)){fscanf(fpi,%lf,&a);fscanf(fpi,%lf,&b);x[i].real=a;x[i].imag=b;i++;}fclose(fpi);while(iN){x[i].real=0;x[i].imag=0;i++;}}}voiddataout(complexx[],longN){longi;for(i=0;iN;i++)printf(%.3lf+j*%.3lf\n,x[i].real,x[i].imag);}voidreverse(complexx[],complexxi[],longn,intm){longi;longtemp;for(i=0;in;i++){temp=t2t(i,m);xi[i].real=x[temp].real;xi[i].imag=x[temp].imag;}}complexcom_add(complexx1,complexx2){complexy;y.real=x1.real+x2.real;y.imag=x1.imag+x2.imag;returny;}complexcom_sub(complexx1,complexx2){complexy;y.real=x1.real-x2.real;y.imag=x1.imag-x2.imag;returny;}complexcom_mul(complexx1,complexx2){complexy;y.real=x1.real*x2.real-x1.imag*x2.imag;y.imag=x1.real*x2.imag+x2.real*x1.imag;returny;}complexcom_div(complexx1,complexx2){complexy;doubledown=x2.real*x2.real+x2.imag*x2.imag;y.real=(x1.real*x2.real+x1.imag*x2.imag)/down;y.imag=(x2.real*x1.imag-x1.real*x2.imag)/down;returny;}voidfft(complexx[],longNN,intM){intm,i,j;longNw,Ni;complexWm,temp;for(m=0;mM;m++){Nw=pow(2,m+1);Ni=pow(2,m);for(i=0;iNi;i++){for(j=i;jNN;j=j+Nw){Wm.real=cos(2*i*PI/Nw);Wm.imag=-sin(2*i*PI/Nw);temp=com_mul(x[j+Ni],Wm);x[j+Ni]=com_sub(x[j],temp);x[j]=com_add(x[j],temp);}}}}voidifft(complexx[],longNN,intM){intm,i,j;longNw,Ni;complexWm,temp;for(m=0;mM;m++){Nw=pow(2,m+1);Ni=pow(2,m);for(i=0;iNi;i++){for(j=i;jNN;j=j+Nw){Wm.real=cos(2*i*PI/Nw);Wm.imag=-sin(2*i*PI/Nw);temp=com_mul(x[j+Ni],Wm);x[j+Ni]=com_sub(x[j],temp);x[j]=com_add(x[j],temp);}}}for(i=0;iNN;i++){x[i].real=x[i].real/NN;x[i].imag=-x[i].imag/NN;}}intmain(){longN,i;intM,flag;doublea,b;printf(Pleaseenterthelengthofthesequence\n);scanf(%ld,&N);M=defM(N);NN=pow(2,M);complexxy[NN];complexxyi[NN];datain(xy,N);for(i=N;iNN;i++){xy[i].real=0;xy[i].imag=0;}reverse(xy,xyi,NN,M);printf(PleasechooseFFT(1)orIFFT(-1)\n);scanf(%d,&flag);if(flag==1)fft(xyi,NN,M);elseif(flag==-1)ifft(xyi,NN,M);printf(Theinitialsequenceis:\n);dataout(xy,N);printf(TheFFTresultsareasfollows:\n);dataout(xyi,N);return0;}