FFT二分蝶计算

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("%18.16lf %18.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;
      }
   
    }

}