注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

鑫淼梦园的博客

圆你的梦想 从这里开始

 
 
 

日志

 
 

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

2014-10-18 01:24:56|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

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

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

 for(i=1;i<MaxValue;i++)
 {

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

(2)、基4的倒位序重排程序

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

void reverse4(complex<double> *Data,int Log4N)
{

 int i,MaxValue,length,MaxPos,CurPos,RevNum;
 complex<double> temp;
 length=1<<(2*Log4N);
 MaxPos=length/4;
 MaxValue=length-1;
 RevNum=0;
 for(i=1;i<MaxValue;i++)
 {
  CurPos=MaxPos;
  while(3*CurPos<(RevNum+1))
  {
   RevNum=RevNum-3*CurPos;
   CurPos=CurPos/4;
  }
  RevNum=RevNum+CurPos;
  if(RevNum<i)
  {
   temp=Data[i];
   Data[i]=Data[RevNum];
   Data[RevNum]=temp;

  }
 }

}


3)、基2时间抽选FFT(DIT2)

把时域的数字信号序列按照奇偶进行分组计算,可以进行如下的变换:

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

从变换结果可以知道,一个长度为N的DFT可以变换成长度为N/2的两个子序列的组合。依次类推,可以直到转为N/2个2点的傅立叶变换的组合。由于经过多次的奇偶抽选,输入数据要按照基2倒序的原则进行重排,输出数据为正常顺序。下面首先用递归的形式进行算法的描述,由于递归方法没有充分利用DIT2算法的优点---原位计算,因此递归形式只是为了描述的清晰,并不是使用的程序。

void dit2rec(complex<double>*InData,complex<double>*OutData,int length,int flag)

{

   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)//单个数的变换,不变,也是递归结束的条件

   {

      if (flag==1)

         OutData[0]=InData[0]/length;

      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(flag*2*pi*i/length));

    OutData[i+length/2]=EvenData- OddData*complex<double>(cos(2*pi*i/length),sin(flag*2*pi*i/length));

  }

  delete[] EvenData,OddData,EvenResult,OddResult;

  return;

}

非递归实现如下(现不考虑输入的倒序数问题):
//data为数据指针,返回值相同,Log2N=log2(length),flag=-1时为正变换,flag=1时为逆变换。

complex<double>* dit2(complex<double>*Data,int Log2N,int flag)
{
 int i,j,k;
 int length=1<<Log2N;
 complex<double> wn,temp;
 int step;
 int index0,index1;

 for(i=1;i<=Log2N;i++)
 {
  step=1<<i;

   for(j=0;j<step/2;j++)
   {
    wn=complex<double>(cos(2*pi/step*j),sin(flag*2*pi/step*j));

    for(k=0;k<length/step;k++)
    {
     index0=k*step+j;
     index1=k*step+step/2+j;
     temp=Data[index0];
     Data[index0]=temp+wn*Data[index1];
     Data[index1]=temp-wn*Data[index1];
    }
   }

 }

 if(flag==1)//如果为反变换,要除以序列的长度
  for(int i=0;i<length;i++)
   Data[i]/=length;
 return Data;

}

当i=1时,也就是第一次循环并没有必要进行复数运算,因为j只能取1,wn为实数,这个时间可以节省。因此可以改进为:

void dit2(complex<double>*Data,int Log2N,int flag)

{

   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(flag*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(flag==1)

   for(i=0;i<length;i++)

     Data[i]/=length;

}


4)、基2频率抽选FFT(DIF2)

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

//DIF2的递归版本实现:

void dif2rec(complex<double>*InData,complex<double>*OutData,int length,int flag)

{

   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)

   {

       if(flag==1)

          OutData[0]=InData[0]/length;

       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(flag*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;

}

 

非递归实现如下(现不考虑输入的倒序数问题):
//data为数据指针,返回值相同,r=log2(length),flag=-1时为正变换,flag=1时为逆变换。

complex<double>*dif2(complex<double>*Data,int Log2N,int flag)
{

 int i,j,k,step,length;
 int index0,index1;
 complex<double> wn,temp;
 length=1<<Log2N;
 for(i=1;i<=Log2N;i++)
 {
  step=1<<(Log2N-i+1);

  for(j=0;j<step/2;j++)
  {
   wn=complex<double>(cos(2*pi/step*j),sin(2*pi/step*j*flag));

   for(k=0;k<length/step;k++)
   {
    index0=k*step+j;
    index1=k*step+step/2+j;
    temp=Data[index0];
    Data[index0]=temp+Data[index1];
    Data[index1]=(temp-Data[index1])*wn;
   }
  }
 }
 if(flag==1)
  for(i=0;i<length;i++)
   Data[i]/=length;
 return Data;
}

和DIT一样,最外层的最后一个循环可以另外独立出来,因为最后一个循环没有必要进行复数运算,这样可以减少复数运算的次数。


5)、基4时间抽选FFT(DIT4)

基四时间抽选快速傅立叶算法

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客 

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客 

DIT4源程序:data为数据指针,返回值相同,Log4N=log4(length),flag=-1时为正变换,flag=1时为逆变换。

complex<double>*dit4(complex<double>*Data,int Log4N,int flag)
{

 int i,j,k,length,step,index0,index1,index2,index3;
 complex<double> wn1,wn2,wn3,img,temp0,temp1,temp2,temp3;
 length=1<<(2*Log4N);
 img=complex<double>(0,1);
 /////////////////////////////////////////////////////////
 for(i=1;i<=Log4N;i++)
 {
  step=1<<(2*i);

  for(j=0;j<step/4;j++)
  {
  
   wn1=complex<double>(cos(2*pi/step*j),sin(flag*2*pi/step*j));
   wn2=complex<double>(cos(4*pi/step*j),sin(flag*4*pi/step*j));
   wn3=complex<double>(cos(6*pi/step*j),sin(flag*6*pi/step*j));

   for(k=0;k<length/step;k++)
   {
    index0=k*step+j;index1=k*step+step/4+j;index2=k*step+step/2+j;index3=k*step+3*step/4+j;
    temp0=Data[index0]+Data[index1]*wn1+Data[index2]*wn2+Data[index3]*wn3;
    temp1=Data[index0]-Data[index1]*wn1*img-Data[index2]*wn2+Data[index3]*wn3*img;
    temp2=Data[index0]-Data[index1]*wn1+Data[index2]*wn2-Data[index3]*wn3;
    temp3=Data[index0]+Data[index1]*wn1*img-Data[index2]*wn2-Data[index3]*wn3*img;
    Data[index0]=temp0;
    Data[index1]=temp1;
    Data[index2]=temp2;
    Data[index3]=temp3;
   }
  }
 }
 if(flag==1)
 {
  for(i=0;i<length;i++)
  {
   Data[i]/=length;
  }
 }
 return Data;

}



6)、基4频率抽选FFT(DIF4)

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

DIF4源程序:data为数据指针,返回值相同,Log4N=log4(length),flag=-1时为正变换,flag=1时为逆变换。

complex<double>* dif4(complex<double>*Data,int Log4N,int flag)
{
 int i,j,k,step,length,index0,index1,index2,index3;
 complex<double> wn1,wn2,wn3,img,temp0,temp1,temp2,temp3;

 length=1<<(2*Log4N);
 img=complex<double>(0,1);
 
 for(i=1;i<=Log4N;i++)
 {
  step=1<<2*(Log4N-i+1);
  for(j=0;j<step/4;j++)
  {
   complex<double> wn1=complex<double>(cos(2*pi/step*j),sin(flag*2*pi/step*j));
   complex<double> wn2=complex<double>(cos(2*pi/step*2*j),sin(flag*2*pi/step*2*j));
   complex<double> wn3=complex<double>(cos(2*pi/step*3*j),sin(flag*2*pi/step*3*j));
   for(k=0;k<length/step;k++)
   {
    index0=k*step+j,index1=k*step+step/4+j;index2=k*step+step/2+j;index3=k*step+step*3/4+j;
    temp0=Data[index0]+Data[index1]+Data[index2]+Data[index3];
    temp1=(Data[index0]-Data[index1]*img-Data[index2]+Data[index3]*img)*wn1;
    temp2=(Data[index0]-Data[index1]+Data[index2]-Data[index3])*wn2;
    temp3=(Data[index0]+Data[index1]*img-Data[index2]-Data[index3]*img)*wn3;

    Data[index0]=temp0;
    Data[index1]=temp1;
    Data[index2]=temp2;
    Data[index3]=temp3;

   }
  }
 }

if(flag==1)
 for(i=0;i<length;i++)
  Data[i]/=length;
 return Data;
}

 


7)、时间抽选分裂基FFT(srfftdit

基4FFT运算量比基2FFT运算量少,但如果不满足N=4^m就无法用基4FFT进行运算,又由于基2算法中偶序号部分FFT运算所需乘法次数比技术序号部分乘法次数少,因而对偶序数部分用基2FFT算法,而奇序数部分用基4算法,则运算量可以得到降低。1984年,法国的P.Dohamel和H.Hollmann把基2分解和基4分解融合在一起,提出“分裂基”FFT算法(split-radix FFT或SRFFT),这种算法在目前已知的N=2^m点FFT算法中乘、加法次数最少,已经非常接近理论上所需乘法次数的最小值,而且它的运算结构与基2FFT相似,没有更多的取指或存储数据的花销,特别适宜用ASIC实现。分裂基算法与混合基算法的不同之处在于,分裂基算法同一级内既有基2部分又有基4部分,而不像组合数FFT中同一级内只用同一种基运算,它们的相似之处在于可以有时间抽取算法,也可以有频率抽取算法。

将输入数据分成三个子序列:

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

做DFT运算有

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

又有上式k范围的限制,还应该补充下列式子才能完整计算N个频率点

【转载】几个快速傅立叶变换算法 - 你梦我圆 - 鑫淼梦园的博客

接着对N/2点DFT和N/4点DFT按以上方法分解成三部分,如此进行到最后不能分解位置。

 

 


8)、频率抽选分裂基FFT(srfftdif

 

 

void srfftdif(complex<double>*Data,int Log2N,int sign)
{
 int i,j,k,len,step2,step4,ix,id,i0,i1,i2,i3;
 complex<double> temp0,temp1,temp2,temp3;
 complex<double> wn1,wn3,img=complex<double>(0,1);
 len=1<<Log2N;

 for(i=1;i<Log2N;i++)
 {
  step2=1<<(Log2N-i+1);
  step4=step2/4;
  for(j=0;j<step4;j++)
  {
   wn1=complex<double>(cos(2.0*pi/step2*j),sin(sign*2.0*pi/step2*j));
   wn3=wn1*wn1*wn1;
   ix=j;
   id=2*step2;
   do
   {
    for(k=ix;k<len-1;k+=id)
    {
     i0=k;
     i1=k+step4;
     i2=k+2*step4;
     i3=k+3*step4;
     temp0=Data[i0]+Data[i2];
     temp1=Data[i1]+Data[i3];
     temp2=((Data[i0]-Data[i2])-img*(Data[i1]-Data[i3]))*wn1;
     temp3=((Data[i0]-Data[i2])+img*(Data[i1]-Data[i3]))*wn3;
     Data[i0]=temp0;Data[i1]=temp1;Data[i2]=temp2;Data[i3]=temp3;
    }
    ix=2*id-step2+j;
    id=4*id;
   } while(ix<(len-1));
  }
 }
 ix=0;
 id=4;
 do
 {
  for(i=ix;i<len-1;i+=id)
  {
   i0=i;
   i1=i+1;
   temp0=Data[i0]+Data[i1];
   temp1=Data[i0]-Data[i1];
   Data[i0]=temp0;
   Data[i1]=temp1;

  }
  ix=2*id-2;
  id=4*id;
 } while(ix<(len-1));
}


9)、高基FFT


10)、混合基FFT


11)、素因子FFT


12)、线性调频Z变换(chirp-z变换)算法

data为数据指针,Log2N=log2(length)

void czt(complex<double>*data,int Log2N,int numofsample,double initfreq,double finalfreq)
{
 int i,len,totallen,Log2Totallen,lenoffft;
 double deltaomiga=2.0*pi*(finalfreq-initfreq)/(numofsample-1.0);
 len=1<<Log2N;totallen=len+numofsample-1;
 ////////////////////////////////////////////
 i=1;Log2Totallen=0;
 while(i<totallen)
 {
  i=i*2;
  Log2Totallen++;
 }
 lenoffft=i;
 complex<double> *gn=new complex<double>[lenoffft];
 complex<double> *hn=new complex<double>[lenoffft];
 ///////////////////////////////////////////
 
 for(i=0;i<len;i++)
          gn[i]=complex<double>(cos(i*i/2.0*deltaomiga),sin(-i*i/2.0*deltaomiga))*complex<double>(cos(2.0*pi*initfreq*i),sin(-2.0*pi*initfreq*i))*data[i];
 for(i=len;i<lenoffft;i++) gn[i]=0;
 reverse2(gn,Log2Totallen);
 dit2(gn,Log2Totallen,-1);

 for(i=0;i<totallen;i++)
         hn[i]=complex<double>(cos(i*i/2.0*deltaomiga),sin(i*i/2.0*deltaomiga));
 for(i=totallen;i<lenoffft;i++) hn[i]=hn[lenoffft-i];
 reverse2(hn,Log2Totallen);
 dit2(hn,Log2Totallen,-1);

 for(i=0;i<lenoffft;i++)
  data[i]=gn[i]*hn[i];
 reverse2(data,Log2Totallen);
 dit2(data,Log2Totallen,1);
 
 for(i=0;i<numofsample;i++)
 {
         data[i]=complex<double>(cos(deltaomiga*i*i/2.0),sin(-deltaomiga*i*i/2.0))*data[i];
 }

        delete[] hn;
        delete[] gn;
}


13)、两个N点实数序列FFT


14)、单个2N点实数序列FFT


15)、共轭变换法IFFT


16)、两个有限长序列的快速卷积


17)、重叠相加法计算有限长和无限长序列的卷积


18)、重叠保留法计算有限长和无限长序列的卷积


19)、快速相关算法

  评论这张
 
阅读(261)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017