G-Spider @2011 初始版。
void butterfly(unsigned int l,int flag,double *x,double *y){ //flag=1 ,fft //flag=-1,ifft unsigned int i,j,r,tmpi,m1,m2,m3,m4,m5,N,N0,k1,k2; double N_1,u,v,tmptheta; N=(unsigned int)1<<l; N0=N>>2; N_1=1.0/(double)N; tmptheta=_2PI*N_1; if(flag==-1) { for(i=0;i<N;i++) { y[i]=-y[i]; } }//======================================================//得到正余弦表 Arrsin[0]=0.0; Arrcos[0]=1.0; Arrsin[1]=sin(tmptheta); Arrcos[1]=cos(tmptheta); if(l>1) { for(i=2;i<=N0;i+=2) {//i=2,4,6,.... tmpi=i>>1; Arrsin[i]= 2.0* Arrsin[tmpi] * Arrcos[tmpi]; Arrcos[i]=1.0- 2.0* Arrsin[tmpi] * Arrsin[tmpi]; //====================================================== Arrsin[i+1]= Arrcos[i] * Arrsin[1] + Arrsin[i] * Arrcos[1]; Arrcos[i+1]=-Arrsin[i] * Arrsin[1] + Arrcos[i] * Arrcos[1]; //printf(".16lf .16lf/n",Arrsin[i],Arrcos[i]); } N0<<=1; for(;i<N0;i+=2) {//i=N0'+2,N0'+4,N0'+6,....,N0 (N0'=N0/2) //sin(pi/2 -x)=cos(x); sin(N0-i)=cos(i); //cos(pi/2-x)=sin(x); cos(N0-i)=sin(i); tmpi=N0-i; Arrsin[i]=Arrcos[tmpi]; Arrcos[i]=Arrsin[tmpi]; //====================================================== tmpi--; Arrsin[i+1]=Arrcos[tmpi]; Arrcos[i+1]=Arrsin[tmpi]; } }//======================================================//二分蝶计算 for(j=0;j<l;) { m1=(unsigned int)1<<j++; m2=m1<<1; m3=N>>j; m4=0; for(i=0;i<m3;i++) { m5=0; for(r=0;r<=(m1>>1);r++) { k1=m4+r; k2=k1+m1;
u=x[k2]*Arrcos[m5]+y[k2]*Arrsin[m5]; v=y[k2]*Arrcos[m5]-x[k2]*Arrsin[m5]; x[k2]=x[k1]-u; x[k1]=x[k1]+u; y[k2]=y[k1]-v; y[k1]=y[k1]+v; m5+=m3; }
m5=m3; for(;r<m1;r++) { k1=m4+r; k2=k1+m1; u=y[k2]*Arrcos[m5]-x[k2]*Arrsin[m5]; v=-(y[k2]*Arrsin[m5]+x[k2]*Arrcos[m5]); x[k2]=x[k1]-u; x[k1]=x[k1]+u; y[k2]=y[k1]-v; y[k1]=y[k1]+v; m5+=m3; } m4+=m2; } }//======================================================//逆变换 if(flag==-1) { for(i=0;i<N;i++) { x[i]= x[i]*N_1; y[i]=-y[i]*N_1; } }
}
//======================
//动态分配
void butterfly(unsigned int l,int flag,double *x,double *y){ //flag=1 ,fft //flag=-1,ifft unsigned int i,j,r,tmpi,m1,m2,m3,m4,m5,N,N0,k1,k2; double N_1,u,v,tmptheta,*Arrsin=NULL,*Arrcos=NULL; N=(unsigned int)1<<l; N0=N>>2; N_1=1.0/(double)N; tmptheta=_2PI*N_1; if(flag==-1) { for(i=0;i<N;i++) { y[i]=-y[i]; } } //======================================================//得到正余弦表 Arrsin=(double *)malloc( ((N>>1)+1)*sizeof(double) ); Arrcos=(double *)malloc( ((N>>1)+1)*sizeof(double) ); if(Arrsin==NULL || Arrcos==NULL) printf("Memory Allocated Error/n"); Arrsin[0]=0.0; Arrcos[0]=1.0; Arrsin[1]=sin(tmptheta); Arrcos[1]=cos(tmptheta); if(l>1) { for(i=2;i<=N0;i+=2) {//i=2,4,6,.... tmpi=i>>1; Arrsin[i]= 2.0* Arrsin[tmpi] * Arrcos[tmpi]; Arrcos[i]=1.0- 2.0* Arrsin[tmpi] * Arrsin[tmpi]; //====================================================== Arrsin[i+1]= Arrcos[i] * Arrsin[1] + Arrsin[i] * Arrcos[1]; Arrcos[i+1]=-Arrsin[i] * Arrsin[1] + Arrcos[i] * Arrcos[1]; } N0<<=1; for(;i<N0;i+=2) {//i=N0'+2,N0'+4,N0'+6,....,N0 (N0'=N0/2) //sin(pi/2 -x)=cos(x); sin(N0-i)=cos(i); //cos(pi/2-x)=sin(x); cos(N0-i)=sin(i); tmpi=N0-i; Arrsin[i]=Arrcos[tmpi]; Arrcos[i]=Arrsin[tmpi]; //====================================================== tmpi--; Arrsin[i+1]=Arrcos[tmpi]; Arrcos[i+1]=Arrsin[tmpi]; } }//======================================================//二分蝶计算 for(j=0;j<l;) { m1=(unsigned int)1<<j++; m2=m1<<1; m3=N>>j; m4=0; for(i=0;i<m3;i++) { m5=0; for(r=0;r<=(m1>>1);r++) { k1=m4+r; k2=k1+m1;
u=x[k2]*Arrcos[m5]+y[k2]*Arrsin[m5]; v=y[k2]*Arrcos[m5]-x[k2]*Arrsin[m5]; x[k2]=x[k1]-u; x[k1]=x[k1]+u; y[k2]=y[k1]-v; y[k1]=y[k1]+v; m5+=m3; }
m5=m3; for(;r<m1;r++) { k1=m4+r; k2=k1+m1; u=y[k2]*Arrcos[m5]-x[k2]*Arrsin[m5]; v=-(y[k2]*Arrsin[m5]+x[k2]*Arrcos[m5]); x[k2]=x[k1]-u; x[k1]=x[k1]+u; y[k2]=y[k1]-v; y[k1]=y[k1]+v; m5+=m3; } m4+=m2; } }//======================================================//逆变换 if(flag==-1) { for(i=0;i<N;i++) { x[i]= x[i]*N_1; y[i]=-y[i]*N_1; } }
}
