几个快速傅立叶变换算法  转

 

几个快速傅立叶变换算法   

2008-03-09 19:25:47|  分类： 默认分类 |  标签：快速卷积  快速相关  混合基  素因子  dft  fft  分裂基   |字号大中小 订阅

离散傅里叶变换(DFT)

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定义进行计算的代码

直接利用DFT的定义进行计算的算法计算量非常大。

//Data为输入数据指针，Log2N=log2(length),flag=-1时为正变换，flag=1时为反变换，变换结果由指针Data指向的原空间返回

void dft(complex*Data,int Log2N,int flag)

{

    int i,j,length;

    complex wn;

    length=1<<Log2N;

    complex*temp=new complex[length];

    for(i=0;i

    {

      temp[i]=0;

      for(j=0;j

      {

        wn=complex(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

       Data[i]=temp[i]/length;

    delete[] temp;

}  

---

2）、倒位序重排

基2、基4和分裂基的DIT、DIF都需要进行倒位序重排。DIT输入为倒位序，输出为正常顺序，DIT输入为正常顺序，输出为倒位序。因此使用DIT要对输入先进行倒位序重排，DIT要对输出进行倒位序重排，才能得到正常顺序的结果。倒序序重排的实现主要是利用加1是从高位开始加，并且进位是向低位进的特点，这刚好和正位序相反，因此可以得到正位序相应的倒位序。还有一种方法，就是得到各个位上的数码，然后倒转再计入权值，就得到倒位序。
（1）、基2的倒位序重排程序

void reverse2(complex *data,int Log2N)
{
 int i,j;
 int RevNum;
 int MaxPos,CurPos,MaxValue;
 complex temp;
 MaxValue=(1<<Log2N)-1;
 MaxPos=1<<(Log2N-1);
 RevNum=0;

 for(i=1;i
 {

  CurPos=MaxPos;
  while((CurPos&RevNum)!=0)
  {
   RevNum=RevNum&(~CurPos);
   CurPos=CurPos>>1;
  }
  RevNum=RevNum|CurPos;
  if (RevNum
  {
   temp=data[RevNum];
   data[RevNum]=data[i];
   data[i]=temp;
  }
 }
}

（2）、基4的倒位序重排程序

//data为数据指针,Log4N=log4(length)。

void reverse4(complex *Data,int Log4N)
{

 int i,MaxValue,length,MaxPos,CurPos,RevNum;
 complex temp;
 length=1<<(2*Log4N);
 MaxPos=length/4;
 MaxValue=