第3章 快速算法FFT

1、3.3 快速傅里叶变换快速傅里叶变换 （FFT） 虽然频谱分析和DFT运算很重要，但在很长一段时间里，由于DFT运算复杂，并没有得到真正的运用。直到1965年库利（Cooley）和图基（Tukey）首次提出DFT运算的一种快速算法以后，情况才发生了根本变化，人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套一套高速有效的运算方法高速有效的运算方法快速付里变换（快速付里变换（FFT）算法）算法。3.3.1减少运算量的思路设设有限长序列 x（n）长度为N，则：101, 1 , 0)()()(NnnkNNkwnxnxDFTkX因x（n）和 都是复数，因此，每计算一个X（k）值，要进行

2、N次复数相乘，和N-1次复数相加，整个DFT运算需要: N*N次复数相乘 N*(N-1)次复数相加。nkNw可见，在DFT计算中，不论是乘法和加法，运不论是乘法和加法，运算量均与算量均与N2成正比成正比。因此，N较大时，运算量十分可观。例，计算N=10点的DFT，需要100次复数相乘，而N=1024点时，需要1048576（一百多万）次复数乘法，如果要求实时处理，则要求有很高的计算速度才能完成上述计算量。 反变换IDFT与DFT的运算结构相同，只是多乘一个常数1/N，所以二者的计算量相同。FFT算法的基本思想：算法的基本思想：两个特性可减少运算量：1）系数 是一个周期函数，它的周期性和对称性可

3、用来改进运算，提高计算效率。 例 又如 因此 利用这些周期性和对称性周期性和对称性，使DFT运算中有些项可合并；nkNjnkNew2nkNnNkNkNnNwww)()(, 12/NNwkNNkNww)2/( 2）利用 的周期性和对称性，把长度为N点的大点数的DFT运算依次分解为若干个小点数的DFT。 因为DFT的计算量正比于N2，N小，计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，但基本上可分为两类：时间抽取法和频率抽取法。nkNw3.3.2基2-FFT算法假定N是2的整数次方，N=2M，M：正整数 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项， r=0

4、,1，N/2-1 )() 12()()2(21rxrxrxrx（1）按时间抽取的按时间抽取的FFT2011)()(NnNnnkNnkNwnxwnx偶数奇数10)()()(NnnkNwnxnxDFTkx求x(n)的DFT：12/0)12(12/02) 12()2(NrkrNNrrkNwrxwrx12/0212/02) 12()2(NrrkNNrkNrkNwrxwwrxnNnNjnNjnNweew222222 /2 1/2 1002212(2 )(21)NNrkkrkNNNrrkNX kxr WWxrWXkW Xk /2 1102(2 )NrkNrXkxr W /2 1202(21)NrkNrXk

5、xrW注意到，X1（k）， X2（k）只有N/2个点，即k=0，1，N/2-1，还必须应用系数 wkN 的周期性和对称性表示 X（k）的 N/2 N-1点：rkNkNrNww22/2 12XX(),0,1,122kNNNX kkWkkkNNkNWW)2(可见，一个N点的DFT被分解为两个N/2点的DFT，这两个N/2点的DFT再合成为一个N点DFT 12,0,1,12kNNX kXkW Xkk 12(),0,1,122kNNNX kXkW Xkk结论：一个结论：一个N点的点的DFT的运算可以分解为两个的运算可以分解为两个N/2点点的的DFT运算运算 ，再经上述算法计算得到。，再经上述算法计算得

6、到。蝶形信号流图蝶形信号流图将X1（k）和X2（k）合成X(k)运算可归结为： bWabWaN=8为例子为例子 两个两个4点点DFT组合成一个组合成一个8点点DFT 按照这个办法，由于N=2M，仍然是偶数，可以被整除，因此可以对两个N/2点的DFT再分别作进一步的分解。即对X1（k）， X2（k）的计算，又可以分别通过计算两个长度为N/4=2点的DFT，进一步节省计算量。这样，一个点的DFT就可以分解为四个点的DFT。 4个个2点点DFT组合成组合成2个个4点点DFT 最后剩下的是2点DFT，它可以用一个蝶形结表示：) 1 ()0() 1 ()0() 1 () 1 ()0() 1 ()0()0

7、(1202xWxxWxXxWxxWxXoNoN这样，一个8点的完整的按时间抽取运算的流图基基2-DIT的的8点点DFT运算流图运算流图 由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“时域抽取法”或“时时间抽取法间抽取法”。时间抽取法时间抽取法FFT的运算特点：的运算特点：1）蝶形运算）蝶形运算2）同址计算）同址计算 3）序数重排）序数重排4）蝶形类型随迭代次数成倍增加）蝶形类型随迭代次数成倍增加 1）蝶形运算）蝶形运算对于N=2M，总是可以通过M次分解最后成为2点的DFT运算。这样构成从x(n)到X(k)的M级运算过程。每一级运算都由N/2个蝶形运算个蝶形运算构

8、成。因此每一级运算都需要N/2次复乘和N次复加，经过时间抽取后M级运算总共需要的运算： 复乘 复加 NNMN2log22NNMN2log 而直接运算时则与N2 成正比。例N=2048，N2=4194304，(N/2)log2N=11264，N2 / (N/2)log2N=392.4。FFT显然要比直接法快得多。2）同址同址计算计算 当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫同址同址计算计算。每一级运算均可在同址存储器同址存储器进行，这种同址同址运算结构可节省存储单元，降低设备成本，还可节省寻址的时间。3）序数重排）序数重排 对按时

9、间抽取FFT的原位运算结构，当运算完毕时，正好顺序存放着正好顺序存放着 X（0），X（1），X（2），X（7），因此可直接按顺序输出按顺序输出。但这种原位运算的输入 x（n）却是按x(0)，x(4)，x(2)，x(6)，x(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是“码位倒置码位倒置”的顺序。例如，原来的自然顺序应是 x(1)的地方，现在放着 x(4)，用二进制码表示这一规律时， 则是在 x（0 0 1）处放着 x（1 0 0）， x（0 1 1）处放着 x（1 1 0）。对于原位运算结构，偶、奇分解过程如下： 表 码位倒置顺序自然顺

10、序二进码表示码位倒置码位倒序0000000010011004201001023011110641000011510110156110010371111117 在实际运算中，一般直接将输入数据 x（n）按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将通过变址运算将自然顺序的存储转换成码位倒置顺序的存储自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。Matlab的bitrevorder（x）函数就是实现自然序的数据排列调整成倒序排列。 目前有许多通用DSP芯片支持这种码位倒置码位倒置的寻址功能。也可以通过移动流图中的某些整行支路 ，把输入x(n

11、)调整成自然序，但是输出频谱序列则成为倒序排列。输出频谱是倒序排列的，可以通过倒序整理环节，成为自然序。但不利于编程和效率，意义不大。 4）蝶形类型随迭代次数成倍增加）蝶形类型随迭代次数成倍增加观察8点FFT的三次迭代运算：第一级迭代，有一种类型的蝶形运算系数W08，两个数据点间隔为1。第二级迭代，有二种类型的蝶形运算系数W08、W28，参加运算的两个数据点间隔为2。第三级迭代，有四类蝶形运算系数W08、W18、W28、W38，参加运算的两个数据点间隔为4。 结论：每迭代一次，蝶形类型增加一倍，数据点间隔也增大一倍。每一级的取数间隔和蝶形类型种类均为2i-1，i=1,2,M。 （2）基）基2频

12、率抽取的频率抽取的FFT（DIF）对于N=2M情况下，输入序列按前后对半分开，这样便将N点DFT写成前后两部分：12/012/)()()(NnNNnnkNnkNWnxWnxkX12/012/0)2()2()(NnNnkNnNnkNWNnxWnx12/0)2/()2/()(NnnkNkNNWNnxWnx把X（k）进一步分解为偶数组和奇数组：奇数为偶数1, 1) 1(, 1)2/(2/kWWkkNNNN12/0)2/()1()()(NnnkNkWNnxnxkX12/02)2/()()2(NnnrNWNnxnxrX12/02/)2/()(NnnrNWNnxnx12/0) 12()2/()(1)X(2

13、rNnrnNWNnxnx12/02/)2/()(NnnrNnNWWNnxnx令 a(n)=x(n)+x(n+N/2) b(n)=x(n)-x(n+N/2wnN这两个序列都是N/2点的序列，将其代入上两式，得12/02/)()2(NnnrNWnarXnrNNnWnbrX2/12/0)() 12(这正是两个N/2点的DFT运算，即将一个N点的DFT分解为两个N/2点的DFT，上式的运算关系可用下图表示. x(n)x(n+N/2)a(n)=x(n)+x(n+N/2)b(n)=x(n)-x(n+N/2)WNnWNn1比较发现，比较发现，DIF和和DIT的流图是流图转置关系，即把输入名和的流图是流图转置

14、关系，即把输入名和输出名对调，流向倒转但增益和结构不变。输出名对调，流向倒转但增益和结构不变。 因此，两者运算量完全一样。因此，两者运算量完全一样。3.3.3 N为组合数的为组合数的FFT算法算法 以上讨论的都是以2为基数的FFT算法，即N=2M。 实际应用时，有限长序列的长度N很大程度上由人为因素确定，因此多数场合可取 N=2M，从而直接使用以 2 为基数的FFT算法。 如N不能人为确定 , N的数值也不是以2为基数的整数次方，处理方法处理方法有两种: 补零: 将x(n)补零 , 使N=2M. 例如N=30，补上x(30)=x(31)=0两点,使N=32=25,这样可直接采用以2为基数M=5

15、的FFT程序。 有限长度序列补零后并不影响其频谱 X(ejw) ，只是频谱的采样点数增加了,上例中由30点增加到32点，所以在许多场合这种处理是可接受的。 如要求准确的N点DFT值，可采用任意数为基数的 FFT 算法 , 其 计算效率低于以2为基数FFT算法。 如 N 为复合数，可分解为两个整数两个整数P与与Q的乘积的乘积，像前面以2为基数时一样，FFT的基本思想是将DFT的点数尽量分小，因此，在N=PQ情况下，也希望将N点的DFT分解为P个Q点DFT或Q个P点DFT，以减少计算量。步骤：kPkknQnn10,kn分别为0, 1,，Q-1; 01,kn分别为0, 1,，P-1。 N点DFT可以

16、重新写成为 100101)(,)()(NnknNWnxkkXkPkXkX1010)(01010101QnPnnQnkPkNWnQnxQnPnnkNPnkNQnkNQnPnnkNPnkNQnkNPQnkNWWWnnxWWWWnnxkkX,考虑到 再令1010nkPQnkNWW令1010010100100110,QnnkQnkNPnnkPWWWnnxkkXPnnkPWnnxnkX,01000,0011001nkQnkNQnWWnkXkkX00001001,nkNWnkXnkXQnnkQWnkXkkX,以P=3，Q=4, N=12为例 (1) 先将 x(n)通过 x(n1Q+n0)改写成 x(n1,

17、n0)。因为 Q=4, n1=0,1,2, n0=0,1,2,3,故输入是按自然顺序的，即x(0,0)=x(0) x(0,1)=x(1) x(0,2)=x(2) x(0,3)=x(3)x(1,0)=x(4) x(1,1)=x(5) x(1,2)=x(6) x(1,3)=x(7)x(2,0)=x(8) x(2,1)=x(9) x(2,2)=x(10) x(2,3)=x(11)(2) 求个点的DFT （）X1(k0,n0)乘以 得到X1(k0,n0)。（）求P个点的DFT，参变量是k020301001110,nnkWnnxnkX304001102001,nnkWnkXkkX（）将X2(k0,k1)

18、通过X(k0+k1P)恢复为X(k)00nkNW运算量：（）求个P点DFT需要P2次复数乘法和P(P-1)次复数加法；（）乘个W因子需要次复数乘法；（）求P个点DFT需要PQ2 次复数乘法和P（-1)次复数加法。总的复数乘法量: QP2+N+PQ2=N(P+Q+1)；总的复数加法量: QP(P-1)+PQ(Q-1)=N(P+Q-2) N=667=23*29，乘法：667*（23+29+1）=35351 加法：667*（23+29-2）=33350N=667，乘法：667*667=444889； 12.58 加法：667*（667-1）=444222； 13.32N=1024；乘法：1024/2

19、*log2（1024）=5120； 加法： 1024*log2（1024） =10240； 上述分解原则可推广至任意基数的更加复杂的情况。例如, 如果N可分解为m个质数因子p1,p2,，pm，即 N=p1p2p3pm则：（1）可把N先分解为两个因子N=p1q1，其中q1=p2p3pm ，并用上述讨论的方法将DFT分解为p1个q1点DFT；（2）将q1分解为q1=p2q2，q2=p3p4pm，然后将每个q1 点DFT再分解为p2个q2点DFT；： 依此类推，：通过m次分解，一直分到最少点数的DFT运算，从而获得最高的运算效率。其运算量近似为N(p1+p2+pm)次复数乘法和复数加法。但计算效率的

20、提高是要以编程的复杂性为代价的，一般较少应用。 p1=p2 = =pm =2，为基基2 FFT 算法算法。当组合数 N=P1P2P3Pm中所有的Pi均为4时，就是基基4FFT算法算法 以N=43为例,第一级运算的一般形式为： ), 3(), 2(), 1 (), 0(111111111111), 3(), 2(), 1 (), 0(), 3(), 2(), 1 (), 0(212121210101010194643404644424043424140404040404011011011011nnxnnxnnxnnxjjjjnnxnnxnnxnnxWWWWWWWWWWWWWWWWnnXnnXnn

21、XnnX30,),(),(04300120101022kWnnnxnnkXknn以N=1024=210，基2FFT 5120 复乘， 基4FFT N*（M-1）=4096 复乘3.3.4 Chirp-z变换变换 采用FFT可以算出N点DFT值，即X（z）在z平面单位圆上的N点等间隔取样值。问题的提出： 不需要整个单位圆上的取样，如：对于窄带信号，只需对信号所在频段进行分析，希望频谱的采样集中在这一频带内，以获得较高的分辨率。 需要对其它围线上的z变换取样，例如：语音信号处理中，需要知道极点所在频率，如果极点离单位圆较远，则其单位圆上的频谱就很平滑，不易判断；如果采样可以沿一条接近这些极点的弧线

22、进行，则在极点所在频率上的频谱将出现明显的尖峰，由此可较准确地测定极点频率。 要求能有效地计算N是素数时，序列的DFT。 算法原理算法原理： 已知x（n），0nN-1令zk=AW-k ，k=0，M-1 ，M ：采样点数，A、W：任意复数00jojoeWWeAA其中：A0表示起始取样点的半径长度，通常A010表示起始取样点z0的相角0表两相邻点之间的等分角W0螺旋线的伸展率，W01则线内缩（反时针），W0=1则表示半径为A0的一段圆弧，若A0=1，这段圆弧则是单位圆的一部分。图 螺旋线采样 计算z变换在采样点 zk 的值 10)()(NnnkkznxzX nk可以用以下表示式来替换10)()(N

23、nnknkWAnxzX)(21222nknknk显然，按照以上公式计算出全部M点采样值需要 NM 次复乘和（N-1）M次复加，当N及M较大时，计算量迅速增加，以上运算可转换为卷积形式，从而可采用FFT进行，这样可大大提高计算速度。k=0,1, ,M-12222)(1022)()(nkNnnnkkWWAnxWzX则222)(,)()(2nnnWnhWAnxng定义：1, 1 , 0)()()()()(210222MkkhkgWnkhngWzXkNnkk则)()()()(2222202202njnonjoneWeWWnhx22)(nWnhxx(n)g(n)X(zk)22nnWA22nW系统的单位脉

24、冲响应频率随时间成线性增加的线性调频信号相似，因此称为Chirp -z变换变换。 算法实现算法实现： 由于输入信号 g（n）是有限长的，长为N，但序列 是无限长的，而计算 0M-1 点卷积 g（k）*h（k）所需要的 h（n）是取值在）是取值在n=-（N-1）M-1那一部分的值那一部分的值，因此，可认为h（n）是一个有限长序列，长为L=N+M-1。 所以，Chirp -z变换为两个有限长序列的线性卷积g（k）*h（k），可用通过FFT来实现。 h（n）的主值序列 可由h（n）作周期延拓后取0nL-1部分值获得，将 与g（n）作圆周卷积后，其输出的前M个值就是Chirp -z变换的M个值。22)

25、(nWnh)(nh)(nhChirp-z变换的计算步骤： (1) 求h（n）的主值序列 (2) 用FFT求 的傅里叶变换: H(k)=FFT , L点 (3) 对x（n）加权并补零1110)(2/)(2/22LnNLWMnWnhnLn)(nh)(nh（4）G（k）=FFT g（n） , L点（5）Y（k）=G（k）H（k） , L点（6）y（n）=IFFTY（k） , L点（7） , 0kM-11010)()(2/2LnNNnWAnxngnn)()(22kYWzXkk利用FFT计算Chirp-z变换计算量估算：乘法数估计 （1）（2）两步可以事先计算，不必实时计算。 （3）（7）两步两次加权共

26、计N+M次复乘,形成Y（k）需L次复乘，一个FFT与一个IFFT需Llog2L次乘，所以总乘法数：L+N+M+Llog2L，直接计算乘法数， NM N及M较大时，用FFT实现Chirp-Z变换，速度上有很大的改进。 Chirp-z变换的特点：变换的特点：1）输入序列的长度N 与 输出序列的长度 M 不需要相等；2）N 及 M 不必是高度复合数，二者均可为素数；3）相邻采样点 zk 之间的角间隔 0 是任意的，即频率分辨率是任意的；4）围线是任意的，不必是 Z 平面上的圆；5）起始点 z0 可任意选定，即可从任意频率上开始对输入数据进行窄带高分辨率分析；6）若 A=1 , M=N , , 可用

27、Chirp-z变换变换计算DFT（即使 N 为素数）。NjeW23.4 离散傅立叶变换的实际应用问题离散傅立叶变换的实际应用问题采样采样截短截短DFT)(ajXa()x t( )xn)(jeX( )( )( )NNx nxnR n)(jNeX)(kXNNkjNeX/2)( 混迭：对连续信号x(t)进行数字处理前，要进行采样nanTttxtx)()()( 采样序列的频谱是连续信号频谱的周期延拓，周期为fs，如采样率过低，不满足采样定理，fs32时，上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因此,也称上述循环卷积方法为快速卷积法。 上述结论适用于 x(n)、h(n) 两序列长度比较接近

28、或相等的情况，如果x(n)、h(n) 长度相差较多，例如， h(n) 为某滤波器的单位脉冲响应，长度有限，用来处理一个很长的输入信号 x(n)，或者处理一个连续不断的信号，按上述方法， h(n) 要补许多零再进行计算，计算量有很大的浪费，或者根本不能实现。 为了保持快速卷积法的优越性，可将 x(n) 分为许多段，每段的长度与 h(n) 接近 ，处理方法有两种：(1) 重叠相加法由分段卷积的各段相加构成总的卷积输出 假定 xi(n) 表示 x（n）序列的第i段 ： 则输入序列可表为：于是输出可分解为： 其中 01) 1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnh

29、nxny)()(*)()(*)()()(*)()(nhnxnyii 1）先对 h（n）及 xi（n）补零，补到具有L点长度，L=N1+N2-1。 一般选 L=2M。 2）用基2 FFT计算 yi（n）=xi（n）*h（n）。 3）重叠部分相加构成最后的输出序列。 由于 yi（n）的长度为L，而xi（n）的长度为N2，因此相邻两段 yi（n）序列必然有L-N2=N1-1点发生重叠。 iinyny)()(计算步骤：a. 事先准备好滤波器参数 H（k）=DFTh(n)，L点b. 用N点FFT计算Xi（k）=DFTxi(n)c. Yi(k)=Xi(k)H(k)d. 用N点IFFT求 yi（n）=IDF

30、TYi(k)e. 将重叠部分相加(2)重叠保留 这种方法和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列保留原来的输入序列值值，这时，如利用DFT实现h(n)和xi(n)的循环卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。 重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。 相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。 两个长为N的实离散时间序列 x（n）与y（n）的互相关函数定义为 1010)()()()()(NnNnxymnynxnymnxmr3.5.4相关函

31、数的FFT计算互相关不具有交换特性 110010()() ( )( ) ()() ( )()NNyxnnNxynrmy nm x ny n x nmx nmy nrm 则可以证明，rxy()的离散傅里叶变换为 Rxy（k）=X*(k)Y(k) 其中 X(k)=DFTx(n), Y(k)=DFTy(n), Rxy(k)=DFTrxy() ， 0kN-1证明：互相关函数定义为 x(n)及y(n)的卷积公式 相比较，我们可以得到相关和卷积的时域关系： 1010)()()()()(NnNnxymnynxnymnxmr)()()()()(10mymxnynmxmfNn)()()()()()()(1010

32、mymxnynmxnymnxmrNnNnxy)()()()(1)()()(10102kYkXIDFTekYkXNnynxrNnNkkNjxy得： Rxy（k）=X*(k)Y(k) )(*)()(DFTkXnRnxNN因 当x（n）=y（n）时，得到x（n）的自相关函数自相关函数为： 维纳辛钦定理： 自相关函数与信号功率谱互为傅立叶变换对。 10)()()(NnxxnxnxrkNjNkekXN2102| )(|1 上面的推导表明：(1)互相关和自相关函数的计算可利用FFT实现。(2) 由于离散傅里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。直接做N点FFT相当于对两个N

33、点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。 实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与循环卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。利用FFT求两个有限长序列的线性相关： 设x(n)和y(n)的长均为N，求线性相关；（1）为了使两个有限长序列的线性相关可用其循环相关代替而不产生混淆，选择周期L2N-1，且L=2m，以使用FFT，将x(n)，y(n)补零至长为L。 （2） 用FFT计算 X(k)，Y(k)（k=0，1，L-1）（3） 计算 R（k）=X*（k）Y（k）（4） 对R（k）作 IFFT，取后N-1项，得 取前N项，得 11)(mNmrxy10)(Nmmrxy解：matlab实现代码如下： x=1 3 -1 1 2 3 3 1;y=2 1 -1 1 2 0 -1 3;k=length(x);xk=fft(x,2*k);yk=fft(y,2*k);rm=real(ifft(conj(xk).*yk);rm=rm(k+2:2*k) rm(1:k);m=(-k+1):(k-1);stem(m,rm)xlabel(m); ylabel(幅度);【例3.5.1】假设x=1 3 -1 1 2 3 3 1，y=2 1 -1 1 2 0 -1 3，求它们的互相关函数。【例3.5.2】