C＋＋编写的FFT算法函数（两种方法）

hyl2323

读研的时候用C++写的两种方法的FFT算法，FFT和IFFT是浙江大学出版社《数值分析》书的方法，记不得名字了。FFT2和IFFT2是一般信号处理书上讲的蝶形算法，都经过调试与matlab校验过。需要说明的是函数中用到的MatrixC是我自己用C++构造的一个复数矩阵类，Complex也是自己构造的复数类，各自封装了一些方法以便简化主程序，提高可读性。
void FFT(MatrixC &fk,MatrixC &Fn,int N){
  const double pi=3.1415926;
  int k,n,P,q;
  P=(int)(log10(N)/log10(2));
  MatrixC temp(fk);
  for(q=1;q<=P;q++)
  {
    for(k=0;k<(int)(pow(2,(P-q)));k++)
    {
      for(n=0;n<(int)(pow(2,(q-1)));n++)
      {
        Fn((int)(k*pow(2,q)+n),0)=temp((int)(k*pow(2,q-1)+n),0)
          +temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0);
        Fn((int)(k*pow(2,q)+n+pow(2,q-1)),0)=(temp((int)(k*pow(2,q-1)+n),0)
          -temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0))
          *Complex(cos(2*pi*k*pow(2,q-1)/N),-sin(2*pi*k*pow(2,q-1)/N));
      }
    }
    temp=Fn;
  }
}

void IFFT(MatrixC &fk,MatrixC &Fn,int N){
  const double pi=3.1415926;
  int k,n,P,q;
  P=(int)(log10(N)/log10(2));
  MatrixC temp(Fn);
  for(q=1;q<=P;q++)
  {
    for(k=0;k<(int)(pow(2,(P-q)));k++)
    {
      for(n=0;n<(int)(pow(2,(q-1)));n++)
      {
        fk((int)(k*pow(2,q)+n),0)=temp((int)(k*pow(2,q-1)+n),0)
          +temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0);
        fk((int)(k*pow(2,q)+n+pow(2,q-1)),0)=(temp((int)(k*pow(2,q-1)+n),0)
          -temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0))
          *Complex(cos(2*pi*k*pow(2,q-1)/N),sin(2*pi*k*pow(2,q-1)/N));
      }
    }
    temp=fk;
  }
  fk=fk*(1.0/N);
}

void FFT2(MatrixC &fk,MatrixC &Fn,int N){
  MatrixC ifk(fk);
  int n=(int)(log10(N)/log10(2));
  int i;
  //inverse the order of the index of fk
  int order,inverse,temp;
  for(order=0;order<N;order++){
   inverse=0;
   for(i=0;i<n;i++){
  temp=(order>>i)&1;
  temp=temp<<(n-1-i);
  inverse=inverse|temp;
   }
      ifk(order,0)=fk(inverse,0);
  }
  //FFT algorithm like butterfly
  Fn=ifk;
  int j,k,kk,d,w;
  const double pi=3.1415926;
  Complex W,Fn_temp;
  for(i=1;i<=n;i++){
   w=(int)pow(2,i);
   d=w/2;
   for(j=0;j<d;j++){
    k=j;
    W=Complex(cos(-2*pi/pow(2,i)*j),sin(-2*pi/pow(2,i)*j));
    while(k<N){
     kk=k+d;
     Fn_temp=Fn(k,0);
        Fn(k,0)=Fn(k,0)+W*Fn(kk,0);
     Fn(kk,0)=Fn_temp-W*Fn(kk,0);
     k+=w;
    }
   }
  }   
}

void IFFT2(MatrixC &fk,MatrixC &Fn,int N){
  MatrixC iFn(Fn);
  int n=(int)(log10(N)/log10(2));
  int i;
  //inverse the order of the index of fk
  int order,inverse,temp;
  for(order=0;order<N;order++){
   inverse=0;
   for(i=0;i<n;i++){
  temp=(order>>i)&1;
  temp=temp<<(n-1-i);
  inverse=inverse|temp;
   }
      iFn(order,0)=Fn(inverse,0);
  }
  //FFT algorithm like butterfly
  fk=iFn;
  int j,k,kk,d,w;
  const double pi=3.1415926;
  Complex W,fk_temp;
  for(i=1;i<=n;i++){
   w=(int)pow(2,i);
   d=w/2;
   for(j=0;j<d;j++){
    k=j;
    W=Complex(cos(-2*pi/pow(2,i)*j),-sin(-2*pi/pow(2,i)*j));
    while(k<N){
     kk=k+d;
     fk_temp=fk(k,0);
        fk(k,0)=fk(k,0)+W*fk(kk,0);
     fk(kk,0)=fk_temp-W*fk(kk,0);
     k+=w;
    }
   }
  }
  fk=fk*(1.0/N);
}

[ 本帖最后由 花如月 于 2007-12-17 18:13 编辑 ]

hyl2323

其实现在数学工具强大，我当时自己编的原因：一是对算法本身感兴趣，二是考验下自己的C++应用能力。收获是：我现在对FFT的理论和物理意义都比较了解了。呵呵。
     有点难以看懂，我现在都看不懂了，当时注释写太少了（不好的编程习惯），不好意思！

[ 本帖最后由 hyl2323 于 2007-12-17 18:06 编辑 ]

ysqcysqc

学习了，多谢分享！！