第24章-快速傅里叶变换原理(FFT)

安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第11页页共共1111页页第第2244章章快快速速傅傅里里叶叶变变换换原原理理（（FFFFTT））在数字信号处理中常常需要用到离散傅立叶变换(DFT)，以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。因此至DFT被发现以来，在很长的一段时间内都不能被应用到实际的工程项目中，直到一种快速的离散傅立叶计算方法——FFT，被发现，离散是傅立叶变换才在实际的工程中得到广泛应用。需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。特别声明：FFT原理的讲解来自网络和书籍。24.1FFT由来24.3直接计算DFT的问题及改进路径24.3改善DFT运算效率的基本途径24.4按时间抽选的基2-FFT算法24.5按频率抽选的基2-FFT算法24.6总结2244..11FFFFTT的的由由来来离散傅里叶变换（DiscreteFourierTransform，DFT）是数字信号处理最重要的基石之一，也是对信号进行分析和处理时最常用的工具之一。在200多年前法国数学家、物理学家傅里叶提出后来以他名字命名的傅里叶级数之后，用DFT这个工具来分析信号就已经为人们所知。历史上最伟大的数学家之一。欧拉是第一个使用“函数”一词来描述包含各种参数的表达式的人，例如：y=f(x)。他是把微积分应用于物理学的先驱者之一。给出了一个用实变量函数表示傅立叶级数系数的方程；用三角级数来描述离散声音在弹性媒介中传播，发现某些函数可以通过余弦函数之和来表达。但在很长时间内，这种分析方法并没有引起更多的重视，最主要的原因在于这种方法运算量比较大。直到1965年，Cooley和Tukey在《计算机科学》发表著名的《机器计算傅立叶级数的一种算法》论文，FFT才开始大规模应用。那个年代，有个肯尼迪总统科学咨询委员会。其中有项研究主题是，对苏联核测试进行检测，Tukey就是其中一员。美国/苏联核测试提案的批准，主要取决于不实地访问核测试设施而做出检测的方法的发展。其中一个想法是，分析离海岸的地震计情况，这种计算需要快速算法来计算DFT。其它应用是国家安全，如用声学探测远距离的核潜艇。所以在军事上，迫切需要一种快速的傅立叶变换算法，这也促进了FFT的正式提出。安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第22页页共共1111页页FFT的这种方法充分利用了DFT运算中的对称性和周期性，从而将DFT运算量从N2减少到N*log2N。当N比较小时，FFT优势并不明显。但当N大于32开始，点数越大，FFT对运算量的改善越明显。比如当N为1024时，FFT的运算效率比DFT提高了100倍。在库利和图基提出的FFT算法中，其基本原理是先将一个N点时域序列的DFT分解为N个1点序列的DFT，然后将这样计算出来的N个1点序列DFT的结果进行组合，得到最初的N点时域序列的DFT值。实际上，这种基本的思想很早就由德国伟大的数学家高斯提出过，在某种情况下，天文学计算（也是现在FFT应用的领域之一）与等距观察的有限集中的行星轨道的内插值有关。由于当时计算都是靠手工，所以产生一种快速算法的迫切需要。而且，更少的计算量同时也代表着错误的机会更少，正确性更高。高斯发现，一个富氏级数有宽度N=N1*N2，可以分成几个部分。计算N2子样本DFT的N1长度和N1子样本DFT的N2长度。只是由于当时尚欠东风——计算机还没发明。在20世纪60年代，伴随着计算机的发展和成熟，库利和图基的成果掀起了数字信号处理的革命，因而FFT发明者的桂冠才落在他们头上。之后，桑德（G.Sand）-图基等快速算法相继出现，几经改进，很快形成了一套高效运算方法，这就是现在的快速傅立叶变换（FFT）。这种算法使DFT的运算效率提高1到2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了良好的条件，大大推进了数学信号处理技术。1984年，法国的杜哈梅（P.Dohamel）和霍尔曼（H.Hollamann）提出的分裂基块快速算法，使运算效率进一步提高。库利和图基的FFT算法的最基本运算为蝶形运算，每个蝶形运算包括两个输入点，因而也称为基-2算法。在这之后，又有一些新的算法，进一步提高了FFT的运算效率，比如基-4算法，分裂基算法等。这些新算法对FFT运算效率的提高一般在50%以内，远远不如FFT对DFT运算的提高幅度。从这个意义上说，FFT算法是里程碑式的。可以说，正是计算机技术的发展和FFT的出现，才使得数字信号处理迎来了一个崭新的时代。除了运算效率的大幅度提高外，FFT还大大降低了DFT运算带来的累计量化误差，这点常为人们所忽略。2244..22直直接接计计算算DDFFTT的的问问题题及及改改进进路路径径2244..22..11问问题题的的提提出出设有限长序列x(n)，非零值长度为N，若对x(n)进行一次DFT运行，共需要多大的运算工作量。2244..22..22DDFFTT的的运运算算量量DFT和IDFT的变换式：X(k)=DFT[x(n)]=x(n)W0≤k≤N−1x(n)=IDFT[X(k)]=1NX(k))W0≤k≤N−1注意：安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第33页页共共1111页页1.x(n)为复数，W=e也为复数。2.DFT与IDFT的计算量相当。下面以DFT为例说明计算量：计算机运算时（编程实现）：k=0X(0)=x(0)W+x(1)W+….+x(N-1)W()k=1X(1)=x(0)W+x(1)W+….+x(N-1)W()k=2X(2)=x(0)W+x(1)W+….+x(N-1)W()..k=N-1X(N-1)=x(0)W()+x(1)W()+….+x(N-1)W()()由上面的结算可得DFT的计算量如下：复数乘法复数加法一个X(k)NN-1N个X(k)(N点DFT)N2N(N-1)复数乘法的计算量：(a+jb)(c+jd)=(ab-bd)+j(bc+ad)实数乘法复数加法一次复乘42一次复加2一个X(k)4N2N+2(N-1)=2(2N-1)N个X(k)(N点DFT)4N下面通过两个实例来说明计算量：例一：计算一个N点DFT，共需N次复乘。以做一次复乘1μs计算，若N=4096，所需时间为（4096）2=16777216μs≈17s例二：石油勘探，有24个通道的记录，每通道波形记录长度为5秒，若每秒抽样500点/秒。1.每道总抽样点数：500*5=2500点。2.24道总抽样点数：24*2500=6万点。3.DFT复乘运算时间：N=(60000)2=36*108次。(60000)2=36*108μs=3600s由于计算量大，且要求相当大的内存，难以实现实时处理，限制了DFT的应用，人们一直在寻求一种能提高DFT运算速度的方法。N个点N次复乘，N-1次复加安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第44页页共共1111页页FFT便是Cooley和Tukey在1995年提出来的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理成为一个新兴的应用学科。2244..33改改善善DDFFTT运运算算效效率率的的基基本本途途径径1、利用DFT运算的系数W的固有对称性和周期性，改善DFT的运算效率。1）对称性2）周期性3）可约性W的特性W=e对称性（W）*=W=W()=W()W∙WW∙W周期性W=W()=W()可约性W=WW=W//ee∙=e=-1特殊点W=1W/=−1W(/)=−W2、将长序列DFT利用对称性和周期性分解为短序列DFT的思路因为DFT的运算量与N2成正比，如果一个大点数N的DFT能分解为若干小点数DFT的组合，则显然可以达到减少运算工作量的效果。安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第55页页共共1111页页复乘：FFT算法的基本思想：利用DFT系数的特性，合并DFT运算中的某些项。把长序列DFT短序列DFT，从而减少运算量。FFT算法分类：时间抽选法DIT:Decimation-In-Time频率抽选法DIF:Decimation-In-Frequency2244..44按按时时间间抽抽选选的的基基22--FFFFTT算算法法2244..44..11算算法法原原理理设输入序列长度为N=2M（M为正整数），将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey算法。其中基2表示：N=2M，M为整数。若不满足这个条件，可以人为地加上若干零值（加零补长）使N2N4N2N点DFTN/2点DFTN/2点DFTN/4点DFTN/4点DFTN/4点DFTN/4点DFT……N2+N2N4+N4N4+N4安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第66页页共共1111页页其达到N=2,M。2244..44..22算算法法步步骤骤分组，变量置换X(k)=DFT[x(n)]=x(n)W0≤k≤N−1先将x(n)按n的奇偶分为两组，作变量置换：当n=偶数时，令n=2r；当n=奇数时，令n=2r+1；得到：x(2r)=x1(r)x(2r+1)=x2(r)r=0,……N2-1分组，变量置换X(k)=DFT[x(n)]=x(n)W=x(n)W为偶数+x(n)W为奇数=x(2r)W⁄+x(2r+1)W()⁄=x(2r)W⁄+Wx(2r+1)W()⁄由于W=e=e/=W/X(k)=x(2r)W⁄+Wx(2r+1)W()⁄=x(2r)W⁄+Wx(2r+1)W/⁄=X(k)+WX(k)其中k=0，1，….N/2–1。X(k)和X(k)只有N/2个点，以N/2为周期；而X(k)却有N个点，以N安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第77页页共共1111页页为周期。要用X(k)和X(k)表达全部的X(k)值。还必须利用W系数的周期特性。因