集美大学学士学位论文
void fft1(double A[2][N],int ifft) {
unsigned long int l,k,i,j,jl0,jl0k,jl1k,reseq,ij; double wrc,wis ,ar,ai,x[2][N]; if(ifirst==0){initial(),ifirst=1;} if(ifft=-1)
for(i=0;i for(j=0;j reseq=inverseq(j,l); ij=kup[l]*reseq; jl0=2*kup[l]*j; wrc=wr[ij]; wis=wi[ij]; for(k=0;k for(i=0;i { j=inverseq(i,M); x[0][i]=A[0][j]; x[1][i]=A[1][j];} for(i=0;i {A[0][i]=x[0][i];A[1][i]=x[1][i];} if(ifft==-1) for(i=0;i /**************************************************/ unsigned long int inverseq(unsigned long int i1,unsigned long int imax) { unsigned long int j,ii,i0,inv; i0=i1;inv=0; 26 叶东明:二维FFT的程序实现及其应用 for(j=0;j {ii=i0%2;inv=2*inv+ii;i0=(i0-ii)/2;} return(inv); } /**************************************************************************/ 附录2:实例2.3.1的主函数程序fft1m1.c #define N 1024 #include \ void functf(double A[2][N]); void main(void) { unsigned long int i; static double A[2][N]; FILE *fp; fp=fopen(\ system(\ functf(A); fprintf(fp,\ printf(\ for(i=0;i fprintf(fp,\ printf(\ for(i=0;i fprintf(fp,\ printf(\ for(i=0;i /**************************************************/ void functf(double A[2][N]) { unsigned long int i; double dt,tt,A0,A1; 27 集美大学学士学位论文 dt=0.1; tt=(double)N*dt; A0=exp(-tt); for(i=0;i A[0][i]=(A1+A0/A1)*dt; A[1][i]=0.0; } } /**************************************************************************/ 附录3:程序文件fft2.c #include unsigned int inverseq(unsigned int i1,unsigned int imax); void fft21(double A[2][N1],int ifft); void fft22(double A[2][N2],int ifft); void initial(void); static int ifirst=0;sign=-1; static unsigned M1,M2,ijup[16],ikup[16],jjup[16],jkup[16]; static double pi,wri[N1],wii[N1],wrj[N2],wij[N2]; /********************initial: *******************************/ void initial(void) { unsigned int i,j,ni; double unit,tt; ni=N1; M1=0; while(ni>1) { ni=ni/2;M1=M1+1;} ni=1; for(j=0;j { printf(\\\nStrike any key to exit!\\n\ getchar();exit(1);} ni=N2;M2=0; while(ni>1){ni=ni/2;M2=M2+1;} ni=1; for(j=0;j {printf(\ getchar();exit(1);} ikup[0]=N1/2;/* kup[l]=2^(M1-l-1)*/ 28 叶东明:二维FFT的程序实现及其应用 ijup[0]=1; for(j=1;j pi=4.0*atan((double)1.0); unit=2.0*pi/((double)N1); tt=0.0; for(i=0;i {wri[i]=cos(tt);wii[i]=sin(tt);tt=tt+unit;} jkup[0]=N2/2;/* jkup[l]=2^(M2-l-1)*/ jjup[0]=1; for(j=1;j unit=2.0*pi/((double)N2); tt=0.0; for(i=0;i {wrj[i]=cos(tt);wij[i]=sin(tt);tt=tt+unit;} } /*row:ifft=1<-->FFT,ifft=-1<-->IFFT**/ /*****************************************************/ void fft21(double A[2][N1],int ifft) { unsigned int l,k,i,j,jl0,jl0k,jl1k,reseq,ij; double wr,wi,ar,ai,x[2][N1]; if(ifft==-1) for(i=0;i for(j=0;j reseq=inverseq(j,l); ij=ikup[l]*reseq; jl0=2*ikup[l]*j; wr=wri[ij]; wi=sign*wii[ij]; for(k=0;k 29 集美大学学士学位论文 A[0][jl0k]=A[0][jl0k]+ar; A[1][jl0k]=A[1][jl0k]+ai;} } } for(i=0;i { j=inverseq(i,M1); x[0][i]=A[0][j]; x[1][i]=A[1][j];} for(i=0;i {A[0][i]=x[0][i];A[1][i]=x[1][i];} if(ifft==-1) for(i=0;i /*column:ifft=1<-->FFT,ifft=-1<-->IFFT*/ /**************************************************/ void fft22(double A[2][N2],int ifft) { unsigned int l,k,i,j,jl0,jl0k,jl1k,reseq,ij; double wr,wi,ar,ai,x[2][N2]; if(ifirst==0){initial();ifirst=1;} if(ifft==-1) for(i=0;i for(j=0;j reseq=inverseq(j,l); ij=jkup[l]*reseq; jl0=2*jkup[l]*j; wr=wrj[ij]; wi=sign*wij[ij]; for(k=0;k for(i=0;i 30