chap4-快速付里叶变换FFT课件

第四章快速付里叶变换(FFT)FastFourietTransformer第一节引言一、快速付里叶变换FFTn n有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT).n nFFT并不是一种新的变换形式,它只是DFT的一种快快 速速 算算 法法.并且根据对序列分解与选取方法的不同而产生了FFT的多种算法.n nFFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。二、二、FFT的的产生产生 1965年库利-图基在计算数学（MathematicofComputation）杂志上发表了著名的“机器计算付里级数的一种算法”文章，提出一种快速计算DFT的方法和计算机程序-揭开了FFT发展史上的第一页。促使FFT算法产生原因还有1967年至1968年间FFT的数字硬件制成，电子数字计算机的条件，使DFT的运算大简化了。FFT出现后使DFT的运算大大简化，运算时间一般可缩短一二个数量级之多，从而使DFT的运算在实际中真正得到了广泛的应用。三、本章主要内容三、本章主要内容n n1.直接计算DFT算法存在的问题及改进途径。n n2.多种DFT算法(时间抽取算法DIT算法，频率抽取算法DIF算法，线性调频Z变换即CZT法)n n3.FFT的应用第二节第二节直接计算直接计算DFT算法存在的算法存在的问题及改进途径问题及改进途径一、直接计算DFT计算量n n问题提出：设有限长序列x(n),非零值长度为N,计算对x(n)进行一次DFT运算，共需多大的运算工作量?1.比较比较DFT与与IDFT之间的运算量之间的运算量其中x(n)为复数，也为复数所以DFT与IDFT二者计算量相同。2.以以DFT为例，计算为例，计算DFT复数运算量复数运算量n n计算一个X(k)(一个频率成分)值，运算量为例k=1则要进行N次复数乘法+(N-1)次复数加法所以，要完成整个DFT运算，其计算量为：n nN*N次复数相乘次复数相乘+N*(N-1)次复数加法次复数加法3.一次一次复数乘法换算成实数运算量复数乘法换算成实数运算量n n复数运算要比加法运算复杂，需要的运算时间长。n n一个复数乘法包括4个实数乘法和个实数乘法和2个实数相个实数相法法。(a+jb)(c+jd)=(ac-bd)+j(bc+ad)n n4次实数乘法2次实数加法4.计算计算DFT需要的实数运算量需要的实数运算量n n每运算一个X(k)的值，需要进行4N次实数相乘和2N+2(N-1)=2(2N-1)次实数相加.整个DFT运算量为：4N2次实数相乘和次实数相乘和2N(2N-1)次次实数相加实数相加.由此看出：直接计算直接计算DFT时，乘法次数与加法次数时，乘法次数与加法次数都是和都是和N2成比例的。当成比例的。当N很大时，所需工作量非很大时，所需工作量非常可观。常可观。例子例子n n例1:当N=1024点时，直接计算DFT需要：N2=1048576次，即一百多万次的复乘运算这对实时性很强的信号处理(如雷达信号处理)来讲，它对计算速度有十分苛刻的要求-迫切需要改进DFT的计算方法，以减少总的运算次数。n n例2：石油勘探，24道记录，每道波形记录长度5秒，若每秒抽样500点/秒，每道总抽样点数=500*5=2500点24道总抽样点数=24*2500=6万点DFT运算时间=N2=(60000)2=36*108次二、改善二、改善DFT运算效率的基本途径运算效率的基本途径能否减少运算量，从而缩短计算时间呢？仔细观察DFT的运算就可看出，利用系数Wnk的以下固有特性，就可减少运算量：（1）WNnk的对称性（2）WNnk的周期性（3）WNnk的可约性将长序列将长序列DFT利用对称性和周期性利用对称性和周期性分解为短序列分解为短序列DFTn n因为DFT的运算量与N2成正比的n n如果一个大点数N的DFT能分解为若干小点数DFT的组合，则显然可以达到减少运算工作量的效果。把N点数据分成二半：其运算量为：再分二半：+=+=这样一直分下去，剩下两点的变换。n n快速付里时变换(FFT)就是在此特性基础上发展起来的，并产生了多种FFT算法，其基本上可分成两大类：n n按抽取方法分：时间抽取法时间抽取法(DIT);频率抽取法(DIF)n n按“基数”分：基基-2FFT算法算法；基-4FFT算法；混合基FFT算法；分裂基FFT算法n n其它方法：线性调频Z变换(CZT法)第三节第三节基基-2按时间抽取的按时间抽取的FFT算法算法Decimation-in-Time(DIT)(Coolkey-Tukey)一、算法原理一、算法原理n n设输入序列长度为N=2M(M为正整数，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey算法。n n其中基数2-N=2M，M为整数.若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到N=2Mn n设一序列设一序列x(n)的长度为的长度为L=9,应加零补长为应加零补长为n nN=24=16 应补应补7个零值。个零值。0123456789101112131415nx(n)二、算法步骤二、算法步骤1.分组，变量置换DFT变换：先将x(n)按n的奇偶分为两组，作变量置换:当n=偶数时，令n=2r;当n=奇数时，令n=2r+1;得到：x(2r)=x1(r);x(2r+1)=x2(r);r=0N/2-1;则其DFT可化为两部分：x(0),x(2)x(2r)偶数点x(1),x(3)x(2r+1)奇数点2.代入代入DFT中中代入代入DFT变换式：变换式：3.求出子序列的求出子序列的DFT上式得：n n一个N点的DFT被分解为两个N/2点DFT。X1(k),X2(k)这两个N/2点的DFT按照：n n再应用W系数的周期性，求出用X1 1(k),X2 2(k)表达的后半部的X(k+N/2)的值。频域中的N个点频率成分为：通过这样分解后运算工作量差不多节省了一半。结论：只要求出（0N/2-1）区间内的各个整数k值所对应的X1(k),X2(k)值，既可以求出(0N-1)整个区间内全部X(k)值，这就是FFT能大量节省计算的关键。时间抽取法蝶形运算流图符号时间抽取法蝶形运算流图符号求求 N=8点点FFT变换变换 按N/2=4，做4点的DFT：将N=8DFT分解成2个4点DFT：可知：时域上：x(0),x(2),x(4),x(6)为偶子序列x(1),x(3),x(5),x(7)为奇子序列频域上：X(0)X(3)，由X(k)给出X(4)X(7),由X(k+N/2)给出按时间抽取将一个按时间抽取将一个N点点DFT分解分解 为两个为两个N/2点点DFT（N=8）偶数序列x1(r)奇数序列x2(r)X(4)X(7)同学们自已写既然这样分解是有效的，由于N=2M，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。按这种方法不断划分下去，直到最后剩下的是2点DFT，两点DFT实际上只是加减运算。且式中：下下图图给给出出N=8时时，将将一一个个N/2点点DFT分分解解成成两两个个N/4点点DFT，由这两个由这两个N/4点点DFT组合成一个组合成一个N/2点点DFT的流图。的流图。N/2点DFT分解为两个N/4点DFT一个一个一个一个2 2点的点的点的点的DFTDFT蝶形流图蝶形流图蝶形流图蝶形流图X2(k)也可进行同样的分解：也可进行同样的分解：式中：一个N=8点DFT分解为四个2点DFT根据上面同样的分析知道，利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，比只用一次分解蝶形组合方式的计算量又减少了大约一半。最后，剩下的是2点DFT，对于此例N=8，就是四个N/4=2点DFT，其输出为X3(k)，X4(k)，X5(k)，X6(k)，k=0,1。故上式不需要乘法。类似地,可求出X4(k)，X5(k)，X6(k)。这种方法的每一步分解，都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。N=8按时间抽取的FFT运算流图三三、FFT算法中一些概念算法中一些概念（1）“级”概念n n将N点DFT先分成两个N/2点DFT，再是四个N/4点DFT直至N/2个两点DFT.每分一次称为“一”级运算。n n因为N=2m，所以N点DFT可分成m级（2）“组”概念每一级都有N/2个蝶形单元，例如：N=8,则每级都有4个蝶形单元。每一级的N/2个蝶形单元可以分成若干组，每一组具有相同的结构，相同的因子分布，第m级的组数为：例：N=8=23，分3级。m=1级,分成四组,每组系数为m=2级,分成二组,每组系数为m=3级,分成一组,每组系数为(3)因子的分布因子的分布四、按时间抽取四、按时间抽取FFT算法的特点算法的特点根据DIT基2-FFT算法原理，能得出任何N=2m点的FFT信号流图。总结按时间抽取法解过程的规律。1.原位运算（同址运算）原位运算（同址运算）从蝶形图可以看出这种运算是很有规律的，其每级（每列）计算都是由N/2个蝶形运算构成的，每一个蝶形结构完成下述基本迭代运算：式中，m表示第m列迭代，k,j为数据所在行数。该式表示的蝶形运算下图所示，由一次复乘和两次复加（减）组成。蝶形运算单元蝶形运算单元X1(0)X2(0)X1(2)X2(2)由N8按时间抽取法FFT运算流图看出，某一列的任何两个节点k和j的节点变量进行蝶形运算后，得到结果为下一列k,j两节点的节点变量，而和其他节点变量无关，因而可以采用原位运算，即某一列的N个数据送到存储器后，经蝶形运算，其结果为下一列数据，它们仍存储在这同一组存储器中，直到最后输出，中间无需其他存储器。也就是蝶形的两个输出值仍放回蝶形的两个输入所在的存储器中。这样只需N个存储单元。下一级的运算仍采用这种原位方式，只不过进入蝶形结的组合关系有所不同。这种原位运算结构可以节省存储单元，降低设备成本。2.倒位序规律倒位序规律观察N8按时间抽取法FFT运算流图，发现当运算完成后，FFT的输出X(k)按正常顺序排列在存储单元中，即按X(0)，X(1)，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，而是按x(0)，x(4)，x(7)的顺序存入存储单元，看起来好像是“混乱无序”的，实际上是有规律的，我们称之为倒位序。造成倒位序的原因是输入x(n)按标号n的偶奇的不断分组。如果n用二进制数表示为(n2n1n0)2（当N=8=23时，二进制为三位），第一次分组，n为偶数（相当于n的二进制数的最低位n0=0）在上半部分，n为奇数（相当于n的二进制数的最低位n0=1）在下半部分。下一次则根据次最低位n1为“0”或是“1”来分偶奇（而不管原来的子序列是偶序列还是奇序列），如此继续分下去，直到最后N个长度为1的子序列。下面的树状图描述了这种分成偶数子序列和奇数子序列的过程。倒位序的形成表表 N=8时的自然顺序二进制数和相应的倒位序二进制时的自然顺序二进制数和相应的倒位序二进制数数 自然顺序（自然顺序（I I）二进制数二进制数 倒位序二进制数倒位序二进制数 倒位序（倒位序（J J）0 01 12 23 34 45 56 67 70000000010010100100110111001001011011101101111110000001001000100101101100010011011010110111111110 04 42 26 61 15 53 37 7看出：码位倒读后的顺序刚好是数据送入计算机内的顺序。看出：码位倒读后的顺序刚好是数据送入计算机内的顺序。N=8 倒位序的变址处理倒位序的变址处理3.蝶形运算两节点的蝶形运算两节点的“距离距离”以8点FFT为例，其输入是倒位序的，输出是自然顺序的。其第一级（第一列）每个蝶形的两节点间“距离”为1，第二级每个蝶形的两节点“距离”为2，第三级每个蝶形的两节点“距离”为4。由此类推得，对N=2m点FFT，当输入为倒位序，输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为2m-1。N=8按时间抽取的FFT运算流图X1(k)X4(k)X2(k)X3(k)X5(k)X6(k)五、按时间抽取的五、按时间抽取的FFT算法与直算法与直接计算接计算DFT运算量的比较运算量的比较n n由前面介绍的按时间抽取的由前面介绍的按时间抽取的FFTFFT运算流图可见：运算流图可见：n n复乘次数：复乘次数：n n复加次数：复加次数：可以得出如下结论：可以得出如下结论：按时间抽取法所需的复乘数和复加数都是与按时间抽取法所需的复乘数和复加数都是与成正比。而直接计算成正比。而直接计算DFTDFT时所需的复乘数与复加数则都时所需的复乘数与复加数则都是与是与N N2 2成正比成正比。例子例子n n看N=8点和N=1024点时直接计算DFT与用基2-按时间抽取法FFT的运算量。看出：当N较大时，按时间抽取法将比直接法快得多。即当点数N越大时，FFT的优点更为明显。对于任何流图，只要保持各节点所连的支路及传输系数不变，则不论节点位置怎么排列所得流图总是等效的，所得最后结果都是x(n)的DFT的正确结果，只只是是数数据据的的提提取取和和存存放放的的次次序序不不同同而而已已。这样就可得到按时间抽取的FFT算法的若干其他形式流图。六、按时间抽取的六、按时间抽取的FFT算法的其他算法的其他形式流图形式流图数据存放的方式不同取用系数的顺序不同图时间抽取、输入自然顺序、输出倒位序的FFT流图第四节第四节 按频率抽取（按频率抽取（DIF）的基的基-2 FFT算法算法 一、一、算法原理算法原理仍设序列点数为N=2M，M为正整数。在把输出X(k)按k的奇偶分组之前，先把输入序列按前一半、后一半分开（不是按偶数、奇数分开），把N点DFT写成两部分。k=0,1,N-1k=0,1,N-1当k为偶数时，(-1)k=1；k为奇数时，(-1)k=-1。因此，按k的奇偶可将X(k)分为两部分:所以为X(2r)为前一半输入与后一半输入之和的N/2点DFT，X(2r+1)为前一半输入与后一半输入之差再与WNn之积的N/2点DFT。令:图频率抽取法蝶形运算单元把把一一个个N点点DFT按按k的的奇奇偶偶分分解解为为两两个个N/2点点的的DFT了了。N=8时，上述分解过程示于下图。时，上述分解过程示于下图。图按频率抽取，将N点DFT分解为两个N/2点DFT的组合 与与时时间间抽抽取取法法的的推推导导过过程程一一样样，由由于于N=2M，N/2仍仍是是一一个个偶偶数数，因因而而可可以以将将每每个个N/2点点DFT的的输输出出再再分分解解为为偶偶数数组组与与奇奇数数组组，这这就就将将N/2点点DFT进进一一步步分分解解为为两两个个N/4 点点DFT。这这两两个个N/4点点DFT的的输输入入也也是是先先将将N/2点点DFT的的输输入入上上下下对对半半分分开开后后通通过过蝶蝶形形运运算算而而形形成成的的，下下图图出出了了这这一步分一步分解的过程。解的过程。图图按频率抽取，将按频率抽取，将N点点DFT分解为两个分解为两个N/4点点DFT的组合的组合这样的分解可以一直进行到第M次（N=2M），第M次实际上是做两点DFT，这N/2个两点DFT的N个输出就是x(n)的N点DFT的结果X(k)。下图表示一个N=8完整的按频率抽取的基-2FFT运算结构。图图 按频率抽取的按频率抽取的FFT（N=8）信号流图信号流图初看起来，DIF法与DIT法的区别是：DIF输入是自然顺序，输出是倒位序的，这与DIT法正好相反。但这不是实质性的区别，因因为为DIF法法与与DIT法法一一样样，都都可可将将输输入入或或输输出出进进行行重重排排，使使二二者者的的输输入入或或输输出出顺顺序序变变成成自自然然顺顺序序或或倒倒位位序序顺顺序序。DIF的基本蝶形与DIT的基本蝶形则有所不同，这才是实质的不同，DIF的的复数乘法只出现在减法之后，复数乘法只出现在减法之后，DIT则是先作复乘后再作加减法。则是先作复乘后再作加减法。频率抽取法(DIF)蝶形运算单元时间抽取法(DIT)蝶形运算单元但是，DIF与DIT就运运算算量量来来说说则则是是相相同同的的，即都有M级（列）运算，每级运算需N/2个蝶形运算来完成，DIF法与DIT法都可进行原原位位运运算算。频率抽取FFT算法的输入是自然顺序，输出是倒位序的。因此运算完毕后，要通过变址计算将倒位序转换成自然位序，然后再输出。转换方法与时间抽取法相同。第五节第五节 离散傅立叶反变换（离散傅立叶反变换（IDFT）的）的快速计算方法快速计算方法0kN-10nN-1需要稍稍改动FFT的程序和参数才能实现IFFT算法。即先将X(k)取共轭，就可直接利用FFT子程序，最后再将运算结果取一次共轭，并乘以1/N,即得到x(n)值。因此FFT运算和IFFT运算就可以共用一个子程序块。4.10 线性卷积与线性相关的线性卷积与线性相关的FFT算法算法 线性卷积的线性卷积的FFT算法算法快速卷积快速卷积 以FIR滤波器为例，因为它的输出等于有限长单位脉冲响应h(n)与有限长输入信号x(n)的离散线性卷积。设x(n)为L点，h(n)为M点，输出y(n)为y(n)也是有限长序列，其点数为L+M-1点。下面首先讨论直接计算线性卷积的运算量。由于每一个x(n)的输入值都必须和全部的h(n)值相乘一次，因而总共需要LM次乘法，这就是直接计算的乘法次数，以md表示为md=LM用FFT算法也就是用圆周卷积来代替这一线性卷积时，为了不产生混叠，其必要条件是使x(n)，h(n)都补零值点，补到至少N=M+L-1，即:0nL-1LnN-10nM-1MnN-1然后计算圆周卷积N这时，y(n)就能代表线性卷积的结果。用FFT计算y(n)的步骤如下：求H(k)=DFTh(n)，N点；求X(k)=DFTx(n),N点；计算Y(k)=X(k)H(k);求y(n)=IDFTY(k)，N点。步骤、都可用FFT来完成。此时的工作量如下：三次FFT运算共需要次相乘，还有步骤的N次相乘，因此共需要相乘次数为比较直接计算线性卷积(简称直接法)和FFT法计算线性卷积(简称FFT法)这两种方法的乘法次数。分两种情况讨论如下：（1）x(n)与h(n)点数差不多。例如，M=L，则N=2M-12M，则这样可得下表：M=LM=L8 8161632326464128128256256512512102410242048204840964096K Kmm0.5720.5720.9410.9411.61.62.782.785.925.928.828.82161629.2429.2453.953.999.999.9当M=8时，FFT法的运算量大于直接法；当M=32时，二者相当；当M=512时，FFT法运算速度可快16倍；当M=4096时，FFT法约快100倍。可以看出，当M=L且M超过32以后，M越长，FFT法的好处越明显。因而将圆周卷积称为快速卷积。（2）当x(n)的点数很多时，即当LM。通常不允许等x(n)全部采集齐后再进行卷积;否则，使输出相对于输入有较长的延时。此外，若N=L+M-1太大，h(n)必须补很多个零值点，很不经济，且FFT的计算时间也要很长。这时FFT法的优点就表现不出来了，因此需要采用分段卷积或称分段过滤的办法。即将x(n)分成点数和h(n)相仿的段，分别求出每段的卷积结果，然后用一定方式把它们合在一起，便得到总的输出，其中每一段的卷积均采用FFT方法处理。有两种分段卷积的办法:重叠相加法和重叠保留法。2重叠相加法重叠相加法设h(n)的点数为M，信号x(n)为很长的序列。我们将x(n)分解为很多段，每段为L点，L选择成和M的数量级相同，用xi(n)表示x(n)的第i段：iLn(i+1)L-1其他ni=0,1,则输入序列可表示成这样，x(n)和h(n)的线性卷积等于各xi(n)与h(n)的线性卷积之和，即每一个xi(n)*h(n)都可用上面讨论的快速卷积办法来运算。由于xi(n)*h(n)为L+M-1点，故先对xi(n)及h(n)补零值点，补到N点。为便于利用基-2FFT算法，一般取N=2mL+M-1，然后作N点的圆周卷积：N由于xi(n)为L点，而yi(n)为（L+M-1）点（设N=L+M-1），故相邻两段输出序列必然有（M-1）个点发生重叠，即前一段的后（M-1）个点和后一段的前（M-1）个点相重叠，如下图所示。应该将重叠部分相加再和不重叠的部分共同组成输出y(n)。图重叠相加法图形图重叠相加法图形和上面的讨论一样，用FFT法实现重叠相加法的步骤如下:计算N点FFT，H(k)=DFTh(n)；计算N点FFT，Xi(k)=DFTxi(n)；相乘，Yi(k)=Xi(k)H(k)；计算N点IFFT，yi(n)=IDFTYi(k)；将各段 yi(n)（包括重叠部分）相加，。重叠相加的名称是由于各输出段的重叠部分相加而得名的。4.11 数字信号处理实现方法数字信号处理实现方法n n1.采用在通用计算机上用软件编程来实现。采用在通用计算机上用软件编程来实现。灵活性高，但一般不能做到实时处理。n n2.利用专用的利用专用的DSP芯片芯片。只针对某一种应用，通过加载数据、控制参数、加控制信号以使它具有有限的可编程能力；对某种专门的应用处理速度快，可做到实时性，但灵活性差。n n3.利用通用利用通用DSP芯片。芯片。它是可以软件编程的DSP，兼具上面两种方法的优点，即灵活，又可做到实时性。采用软件编程来实现采用软件编程来实现n n 变址变址n nMM级递推级递推P187 图图433利用专用的利用专用的DSP芯片芯片n n市场上推出专门用于FFT，FIR滤波器，卷积等专用数字芯片。n n如：BB公司：DF17XX系列n nMAXIM公司：MAXIM27X,MAXIM28Xn nNational公司：National-SEMI系列：MF系列。n n其软件算法已在芯片内部用硬件电路实现，使用者只需给出输入数据，可在输出端直接得到数据。利用通用利用通用DSP芯片芯片n n内部带有乘法器，累加器，采用流水线工作方式及并行结构，多总线速度快。配有适于信号处理的指令(如FFT指令)等。n n目前市场上的DSP芯片有：n n美国德州仪器公司(TI):TMS320CX系列占有90%n n还有AT&T公司dsp16,dsp32系列n nMotorola公司的dsp56x,dsp96x系列n nAD公司的ADSP21X,ADSP210X系列