幾種快速傅里葉變換(FFT)的C++實現

鏈接：http://blog.csdn.net/zwlforever/archive/2008/03/14/2183049.aspx
一篇不錯的FFT 文章，收藏一下。
 DFT的的正變換和反變換分別為（1）和（2）式。假設有N個數據，則計算一個頻率點需要N次復數乘法和N-1次復數加法，整個DFT需要N*N次復數乘法和N（N-1）次復數加法；由於一次的復數乘法需要進行4次的實數乘法和2次的復數加法，一次的復數加法需要兩次的實數加法，因此整個DFT需要4*N*N次的實數乘法和2*N（N-1）+2*N*N≈4*N*N次的復數加法。當N比較大時，所需的計算工作量相當大，例如N=1024時大約需要400萬次乘法運算，對於實時信號處理來說，將對計算設備提出非常苛刻的要求，於是就提出如何能夠減少計算DFT的運算量的問題。

1965年，庫力和圖基在《計算數學》雜志上發表《機器計算傅立葉級數的一種算法》，此文是最早提出FFT算法的。此后關於DFT的快速算法稱為人們研究的熱點課題，也正是FFT算法的出現，才使得數字信號處理能夠應用於實時場合並進一步推動數字信號處理的發展和應用。

 

大多數的FFT算法都是利用（3）式的周期性、共軛對稱性、可壓縮和擴展性的特點來壓縮計算量。。

1）、根據DFT定義進行計算的代碼

//Data為輸入數據指針，Log2N=log2(length),flag=-1時為正變換，flag=1時為反變換，變換結果同樣由指針Data指向的空間返回
void dft(complex<double>*Data,int Log2N,int flag)
{
        int i,j,length;
        complex<double> wn;
        length=1<<Log2N;
        complex<double>*temp=new complex<double>(length);
    for(i=0;i<length;i++)
   {
               temp[i]=0;
               for(j=0;j<length;j++)
       {
                    wn=complex<double>(cos(2.0*pi/length*i*j),sin(flag*2.0*pi/length*i*j));
                    temp[i]+=Data[j]*wn;    
                }           
       }
       if(flag==1)
           for(i=0;i<length;i++)
               Data[i]=temp[i]/length;
        delete[] temp;
}  

2）、基2時間抽選FFT

把時域的數字信號序列按照奇偶進行分組計算，可以進行如下的變換，從變換結果可以知道，一個長度為N的DFT可以變換成長度為N/2的兩個子序列的組合。依次類推，可以直到轉為N/2個2點的傅立葉變化的組合。不過這時的輸入應該為以2為基的倒位序。

由於經過多次的奇偶抽選，輸入數據要按照基2倒序的原則進行重排，輸出數據為正常順序，倒序算法另外敘述。下面首先用遞歸的形式進行算法的描述，由於遞歸方法沒有充分利用DIT2算法的優點---原位計算，因此遞歸形式只是為了描述的清晰。

 

void dit2rec(complex<double>*InData,complex<double>*OutData,int length,int sign)
{
   complex<double>*EvenData=new complex<double>(length/2);
   complex<double>*OddData  =new complex<double>(length/2);
   complex<double>*EvenResult=new complex<double>(length/2);
   complex<double>*OddResult=new complex<double>(length/2);
   int i,j;
   if(length==1)
   {
      OutData[0]=InData[0]/N;
      return;
   }
  for(i=0;i<length/2;i++)
  {
    EvenData[i]=InData[2*i];
    OddData[i]=InData[2*i+1];
  }
  dit2rec(EvenData,EvenResult,length/2,sign);
  dit2rec(OddData,EvenResult,length/2,sign);
  for(i=0;i<length/2;i++)
  {
    OutData[i]=EvenData+OddData*complex<double>(cos(2*pi*i/length),sin(sign*2*pi*i/length));
    OutData[i+length/2]=EvenData- OddData*complex<double>(cos(2*pi*i/length),sin(sign*2*pi*i/length));
  }
  delete[] EvenData,OddData,EvenResult,OddResult;
  return;
}

 

非遞歸實現如下（現不考慮輸入的倒序數問題）：

void dit2(complex<double>*Data,int Log2N,int sign)
{
   int i,j,k,step,length;
   complex<double> wn,temp,deltawn;
   length=1<<Log2N;
   for(i=1;i<=Log2N;i++)
   {
      wn=1;step=1<<i;deltawn=complex<double>(cos(2*pi/step),sin(sign*2.0*pi/step))；
      for(j=0;j<step/2;j++)
      {        
        for(i=0;i<length/step;i++)
        {
           temp=Data[i*step+step/2+j]*wn;
           Data[i*step+step/2+j]=Data[i*step+j]-temp;
           Data[i*step+j]=Data[i*step+j]+temp;
         }
        wn=wn*deltawn;
      }
    }
    if(sign==-1)
       for(i=0;i<length;i++)
            Data[i]/=length;
 }

 

當i=1時，也就是第一次循環並沒有必要進行復數運算，因為j只能取1，wn為實數，這個時間可以節省。因此可以改進為:

void dit2(complex<double>*Data,int Log2N,int sign)
{
   int i,j,k,step,length;
   complex<double> wn,temp,deltawn;
  length=1<<Log2N;
   for(i=0;i<length;i+=2)
   {
      temp=Data[i];
      Data[i]=Data[i]+Data[i+1];
      Data[i+1]=temp-Data[i+1];
   }
   for(i=2;i<=Log2N;i++)
   {
      wn=1;step=1<<i;deltawn=complex<double>(cos(2.0*pi/step),sin(sign*2.0*pi/step));；
      for(j=0;j<step/2;j++)
      {        
        for(i=0;i<length/step;i++)
        {
           temp=Data[i*step+step/2+j]*wn;
           Data[i*step+step/2+j]=Data[i*step+j]-temp;
           Data[i*step+j]=Data[i*step+j]+temp;
         }
         wn=wn*deltawn;
      }
   }
   if(sign==1)
   for(i=0;i<length;i++)
     Data[i]/=length;
}

3）、基2頻率抽選FFT

//DIF2的遞歸版本實現：
void dif2rec(complex<double>*InData,complex<double>*OutData,int length,int sign)
{
   complex<double>* LData=new complex<double>(length/2);
   complex<double>* LResult=new complex<double>(length/2);
   complex<double>* RData=new complex<double>(length/2);
   complex<double>* RResult=new complex<double>(length/2);
   complex<double> temp;
   int i;
if(length==1)
   {
       OutData[0]=InData[0];
       return;
}
for(i=0;i<length/2;i++)
{
   LData[i]=InData[i];
   RData[i]=InData[i+length/2];
}
for(i=0;i<length/2;i++)
{
   temp=LData[i];
   LData[i]=LData[i]+RData[i];
   RData[i]=(temp-RData[i])* complex<double>(cos(2*pi*i/length),sin(sign*2*pi*i/length))
}
dit2rec(LData,LResult,length/2,sign);
dit2rec(RData,RResult,length/2,sign);
   for(i=0;i<length/2;i++)
   {
      OutData[2*i]=LResult[i];
      OutData[2*i+1]=RResult[i];
}
return;
}

 

 

//非遞歸實現如下（現不考慮輸入的倒序數問題）：
void dif2(complex<double>*InData,int r,int sign)
{
int length=1<<r;
int i,j,k,step;
complex<double> temp,wn;
for(i=1;i<=r;i++)
{
   step=1<<(r-i+1);
   wn=1;
   for(j=0;j<step/2;j++)
   {
      for(k=0;k<step/2;k++)
      {
         temp=InData[k*step+j];
         InData[k*step+j]=InData[k*step+j]-InData[k*step+step/2+j];
         InData[k*step+step/2+j]=(temp-InData[k*step+step/2+j])*wn;
}
wn=wn*complex<double>(cos(2*pi/step*j),sin(sign*2*pi/step*j));
}
}
}

 

和DIT一樣，最外層的最后一個循環可以另外獨立出來，因為最后一個循環沒有必要進行復數運算，這樣可以減少復數運算的次數。

基四時間抽選快速傅立葉算法