fft.rar_fft的C语言代码

#define PI 3.14159265358979 #define numlength 2048 #include <math.h> void radix2(int n, short *xy, short *w); void bitrev_cplx(int *x, short *index, int nx); void digitrev_index(short *index, int n, int radix); static short index1[64],w1[numlength]; static short x1[numlength*2]; static int x2[numlength]; static int radix; int nx1; int i; main() { double delta; radix=2; nx1=numlength; digitrev_index(index1, nx1, radix); delta=2*PI/nx1; for(i=0;i<nx1/2;i++) { w1[2*i]=32767*(-cos((double)i*delta)); w1[2*i+1]=32767*(-sin((double)i*delta)); } for(i=0;i<nx1;i++) { x1[2*i]=(short)((cos(PI*i/10.0)+cos(PI*i/20.0)+cos(PI*i/5.0))*0x80); x1[2*i+1]=0; } radix2(nx1, x1, w1); for(i=0;i<nx1;i++) { x2[i]=sqrt(x1[2*i]*x1[2*i]+x1[2*i+1]*x1[2*i+1]); } bitrev_cplx(x2, index1, nx1); } void digitrev_index(short *index, int n, int radix) { int i,j,k; short nbits, nbot, ntop, ndiff, n2, raddiv2; nbits=0; i=n; while(i>1) { i=i>>1; nbits++; } raddiv2=radix>>1; nbot=nbits>>raddiv2; nbot=nbot<<raddiv2-1; ndiff=nbits&raddiv2; ntop=nbot+ndiff; n2=1<<ntop; index[0]=0; for(i=1,j=n2/radix+1;i<n2-1;i++) { index[i]=j-1; for(k=n2/radix; k*(radix-1)<j; k/=radix) j-=k*(radix-1); j+=k; } index[n2-1]=n2-1; } void bitrev_cplx(int *x, short *index, int nx) { int i; short i0,i1,i3; short j0,j1,j3; int xi0,xi1,xi3; int xj0,xj1,xj3; short t; int a,b,ia,ib,ibs; int mask; int nbits, nbot, ntop, ndiff, n2, halfn; nbits=0; i=nx; while(i>1) { i=i>>1; nbits++; } nbot=nbits>>1; ndiff=nbits&1; ntop=nbot+ndiff; n2=1<<ntop; mask=n2-1; halfn=nx>>1; for(i0=0; i0<halfn; i0+=2) { b=i0&mask; a=i0>>nbot; if(!b)ia=index[a]; ib=index[b]; ibs=ib<<nbot; j0=ibs+ia; t=i0<j0; xi0=x[i0]; xj0=x[j0]; if(t) { x[i0]=xj0; x[j0]=xi0; } i1=i0+1; j1=j0+halfn; xi1=x[i1]; xj1=x[j1]; x[i1]=xj1; x[j1]=xi1; i3=i1+halfn; j3=j1+1; xi3=x[i3]; xj3=x[j3]; if(t) { x[i3]=xj3; x[j3]=xi3; } } } void radix2(int n, short xy[], short w[]) { short n1,n2,ie,ia,i,j,k,l; short xt,yt,c,s; n2=n; ie=1; for(k=n;k>1;k=(k>>1)) { n1=n2; n2=n2>>1; ia=0; for(j=0;j<n2;j++) { c=w[2*ia]; s=w[2*ia+1]; ia=ia+ie; for(i=j;i<n;i+=n1) { l=i+n2; xt=xy[2*l]-xy[2*i]; xy[2*i]=xy[2*i]+xy[2*l]; yt=xy[2*l+1]-xy[2*i+1]; xy[2*i+1]=xy[2*l+1]+xy[2*l+1]; xy[2*l]=(c*xt+s*yt)>>15; xy[2*l+1]=(c*yt-s*xt)>>15; } } ie=ie<<1; } }