数字信号处理 吴镇扬_CHP3快速傅立叶变换

1、 快速傅里叶变换快速傅里叶变换(FFT)(FFT)是求解离散傅里叶变换是求解离散傅里叶变换(DFT)(DFT)的快速算法。的快速算法。问题的提出问题的提出 直接计算直接计算N N点点DFTDFT需要的计算量是多少？需要的计算量是多少？ 计算一个计算一个X(k)X(k)需要需要N N次复数乘法和次复数乘法和N N一一1 1次复数次复数加法。算出全部加法。算出全部N N点点X(k)X(k)共需共需N N2 2次复数乘法和次复数乘法和N(NN(N一一1)1)次复数加法次复数加法. . 总运算量近似地正比于总运算量近似地正比于N N2 2 。当。当N N值很大（如值很大（如2-D2-D图像处理），运算

2、量将非常庞大；同时，对于图像处理），运算量将非常庞大；同时，对于实时性很强的信号处理来说，必将对计算速度有实时性很强的信号处理来说，必将对计算速度有十分苛刻的要求。为此，需要改进对十分苛刻的要求。为此，需要改进对DFTDFT的计算方的计算方法，以减少总的运算次数。法，以减少总的运算次数。在正交矩阵中，虽然有在正交矩阵中，虽然有N N2 2个元素，但只有个元素，但只有N N个不同个不同的值，且有些取值特别简单，主要由于旋转因子具的值，且有些取值特别简单，主要由于旋转因子具有如下的特点：有如下的特点：对称性对称性周期性周期性下面以四点下面以四点DFT为例来说明快速算法的思路。为例来说明快速算法的思

3、路。01.1W 2.1,1NmNWW3.N rrWW24.1NW 25.NrrWW 如何充分利用这些关系* ( )( )( )( )( )( )( )( )=jjjjjjjjjXxeeeXxXxeeeXxeeeWWWW2221 12 13 14442221 22 23 24442221 32 33 34441111111100111221331111111111111( )( ) ( )( )xxxx 0123( )( )( )( ) ( )( )( )( ) ( )( ) ( )( ) ( )( ) ( )( )= ( )( ) ( )( ) ( )( ) ( )( ) XxXWWxXxXW

4、WxxxxxxxxxWxxxxxxxxW111111011110111221111131130213021302130213交换矩阵第二列和第三列得交换矩阵第二列和第三列得从上面的结果可以看出从上面的结果可以看出, ,利用对称性和周期性，求利用对称性和周期性，求四点四点DFTDFT只需要一次复数乘法，称为只需要一次复数乘法，称为Coolkey-TukeyCoolkey-Tukey算法。算法。11(0) (0)(2) (1)(3)(1) (0)(2) (1)(3)(2) (0)(2) (1)(3)(3) (0)(2) (1)(3)XxxxxXxxxxWXxxxxXxxxxW11(0)(2)xx(

5、1)(3)xx(0)(1)XX(2)(3)XX111W(0)(2)xx(0)(2)xx(1)(3)xx(1)(3)xx概述概述u算法分类：算法分类：N N为为2 2的整次幂的整次幂：按基数分为按基数分为基基-2FFT-2FFT算法算法、基基-4FFT-4FFT算法算法、混合基混合基FFTFFT算法算法、分裂基分裂基FFTFFT算法算法; ;当当N N不是不是2 2的整次幂的整次幂: :典型的有典型的有Winograd Winograd 算法算法. .按抽取方法分：按抽取方法分：时间抽取时间抽取(Decimation(DecimationininTimeTime，简称，简称DIT)DIT)；频率

6、抽取；频率抽取(Decimation(DecimationininFrequencyFrequency，简称，简称DIF)DIF) 为了将大点数的为了将大点数的DFT DFT 分解为小点数的分解为小点数的DFTDFT运算，要求序列的长度运算，要求序列的长度N N为为N N2 2M M (M (M为正为正整数整数) )。该情况下的变换称为基。该情况下的变换称为基2 FFT2 FFT。 N点DFTN/2点 DFTN/4点 DFT 2点 DFT 1个 2个 4个 N/2个问题是如何分最有效？可以对时间变量分问题是如何分最有效？可以对时间变量分（DIT)，也可对频率变量分，也可对频率变量分(DIF）/

7、()/( )()()()()NNrkrkNNrrNNrkkrkNNNrrX kxr WxrWxr WWxrW2 12 1221002 12 12200221221()(), ,/xrxrrN2210 12 1和，基本思路：从时域将基本思路：从时域将N N点序列点序列x(n)x(n)按奇偶项分解为按奇偶项分解为两组，分别计算两组两组，分别计算两组N/2N/2点点DFTDFT，然后再合成一个，然后再合成一个N N点点DFTDFT，按此方法继续下去，直到，按此方法继续下去，直到2 2点点DFTDFT，从而减，从而减少运算量。少运算量。算法具体步骤：算法具体步骤：一、算法的推导一、算法的推导1 1、序

8、列、序列x(n)x(n)按奇偶项分解为两组，将按奇偶项分解为两组，将DFTDFT运算也运算也相应分为两组相应分为两组则则/( )()( )()( )( )( )NrkNrNrkNrkNA kxr WB kxrWX kA kW B k2 1202 1202212 2、两个、两个N/2N/2点的点的DFTDFT合成一个合成一个N N点点DFTDFT问题：问题：A(k)A(k)，B(k)B(k)都只有都只有N/2N/2个点，怎样得到个点，怎样得到X(k)X(k)的后的后N/2N/2点？利用周期性和对称性得点？利用周期性和对称性得 (/ )( )( )kNX kNA kW B k23.3.1 3.3.

9、1时间抽取（时间抽取（DIT）基）基2FFT算法算法3 3、继续分解（一直分解到两点、继续分解（一直分解到两点DFTDFT变换）变换） A(K)和和B(K)仍是高复合数仍是高复合数(N2)的的DFT，我们可，我们可按上述方法继续以分解。令按上述方法继续以分解。令r2l，r2l十十1，l0，1，N4-1，则，则A(K)和和B(K)可分别表示为可分别表示为3.3.1时间抽取（时间抽取（DIT）基）基2FFT算法算法3.3.1时间抽取（时间抽取（DIT）基）基2FFT算法算法令令则则3.3.1时间抽取（时间抽取（DIT）基）基2FFT算法算法同理，令同理，令则则按此方法一直分解下去直到按此方法一直分

10、解下去直到2 2点点DFTDFT，当，当N=8N=8时，如下：时，如下：3.3.1时间抽取（时间抽取（DIT）基）基2FFT算法算法下面通过讨论寻找下面通过讨论寻找FFTFFT的一般规律。的一般规律。二、算法的讨论二、算法的讨论1 1、“级级”的概念的概念 在分解过程中，每分一次，在分解过程中，每分一次，称为一级运算称为一级运算。因为。因为M=log2M=log2N N，所以，所以N N点点DFTDFT可以分解为可以分解为M M级，级，按抽取算法的信号流图中来定义，从左到右分按抽取算法的信号流图中来定义，从左到右分别称为别称为0 0级、级、1 1级到级到M-1M-1级。级。( )( )( )(

11、 )( )( )rmmNmrmmNmxpxpW xqxqxpW xq112 2、蝶形单元、蝶形单元 在算法的信号流图中，第在算法的信号流图中，第m m级存在这种运级存在这种运算，这种结构几何形状像蝴蝶，称为蝶形单元算，这种结构几何形状像蝴蝶，称为蝶形单元p p、q q是参于本蝶形单元运算的上、下节点的序号。由是参于本蝶形单元运算的上、下节点的序号。由于第于第m m级序号的两点只参于这一个蝶形单元的运算，级序号的两点只参于这一个蝶形单元的运算，其输出在第其输出在第m m十十l l级。且这一蝶形单元也不再涉及别的级。且这一蝶形单元也不再涉及别的点。由于这一特点，在计算机编程时，我们可将蝶形点。由于

12、这一特点，在计算机编程时，我们可将蝶形单元的输出仍放在输入数组中，这一特点称为单元的输出仍放在输入数组中，这一特点称为“同址同址运算运算”。3.3.1时间抽取（时间抽取（DIT）基）基2FFT算法算法loglogcaNMNMNMNNMN2222 由于一级都含有由于一级都含有N/2N/2个蝶形单元，每个蝶形单元个蝶形单元，每个蝶形单元需要需要1 1次复数乘法和两次复数加法，因此完成次复数乘法和两次复数加法，因此完成loglog2 2N N级级共需要的复数乘法和加法分别为共需要的复数乘法和加法分别为 直接计算直接计算DFTDFT时所需的复乘数与复加数都是与时所需的复乘数与复加数都是与N2N2成正比

13、的。所以采用成正比的。所以采用FFTFFT算法使运算量大大减少。显算法使运算量大大减少。显然，然，N N值愈大，节省的运算量愈多。值愈大，节省的运算量愈多。3 3、“组组”的概念的概念 在分解过程中，每一级的在分解过程中，每一级的N/2N/2个蝶形单元个蝶形单元可以分成若干组，每一组具有相同的结构和可以分成若干组，每一组具有相同的结构和W W因子分布。第因子分布。第m m级可分成级可分成N/2N/2m+1m+1组。组。例：例：N=8=23，分分3级。级。第一级的分组及第一级的分组及Wr因子如下：因子如下：m=0级级,分成四组：因子为分成四组：因子为m=1级级,分成二组分成二组,因子因子为为m=

14、2级级,分成一组分成一组,因子因子为为/NWWW000428/,NNNNWWW WW W0102022288,NNNNW W W WW W W W0123012388884 4、W Wr r因子的分布因子的分布 由上分析可知由上分析可知结论：结论：每由后向前（每由后向前（m由由M-1-0级）推进一级，则级）推进一级，则此系数为后级系数中偶数序号的那一半。此系数为后级系数中偶数序号的那一半。5 5、码位倒置、码位倒置 在在FFTFFT算法中，输出的频谱依照正常次序算法中，输出的频谱依照正常次序排列，但输入的序列排列，但输入的序列x(n)x(n)是按奇偶分开的，分是按奇偶分开的，分开的规律，以开的

15、规律，以N=8N=8为例，按如下方法进行排序为例，按如下方法进行排序（1 1）、将）、将x(n)x(n)的序号写成二进制的序号写成二进制 x(000),x(001),x(000),x(001), , x(110) ，x(111)x(111)。（2 2）将二进制的码进行翻转，得）将二进制的码进行翻转，得 x(000),x(100),x(000),x(100), , x(011) ， x(111)x(111)。（3 3）将二进制的翻转码转换为对应的十进制）将二进制的翻转码转换为对应的十进制 x(0),x(4),x(0),x(4), x(3), x(3)，x(7)x(7)。 这就是按奇偶抽取得到的顺

16、序。这就是按奇偶抽取得到的顺序。说明：说明： 在上述的基在上述的基2FFT2FFT算法中，由于每一算法中，由于每一步分解都是按输入序列步分解都是按输入序列x(n)x(n)在时域上的次在时域上的次序是属于偶数还是奇数来抽取的，所以称序是属于偶数还是奇数来抽取的，所以称为为“按时间抽取法按时间抽取法”或或“时间抽取时间抽取”。 上述的基上述的基2FFT2FFT算法中，抽取也可在算法中，抽取也可在频域进行，引出频率抽取（频域进行，引出频率抽取（DIFDIF）基）基2FFT2FFT算法。算法。 设输入序列长度为设输入序列长度为N=2M(M为正整数为正整数)，频率抽取法将输入序列不是按奇、偶分组，频率抽

17、取法将输入序列不是按奇、偶分组，而是按前后对半分开，这样可将而是按前后对半分开，这样可将N点点DFT写成前后两部分写成前后两部分;将该序列的频域的输出序将该序列的频域的输出序列列X(k)(也是也是N点序列，按其点序列，按其频域顺序的奇频域顺序的奇偶分解偶分解为越来越短的子序列，称为为越来越短的子序列，称为基基2按按频率抽取的频率抽取的FFT算法算法。也称为。也称为Sander-Tukey算法。算法。算法分析算法分析 现将输入现将输入x(n)按按n的顺序分前后两部分的顺序分前后两部分:前半子序列前半子序列x(n),0nN/2-1;后半子序列后半子序列x(n+N/2),0nN/2-1;例：例：N=

18、8时，前半序列为：时，前半序列为：x(0),x(1),x(2),x(3); 后半序列为：后半序列为： x(4),x(5),x(6),x(7); 考虑考虑N点的点的DFT,由由DFT定义定义得得()()( ) ( )()( )() ( )()NNNnknknkNNNnnNNNnknkNNnnNNknkNNnNkNNjkkNj kNNNkNNX kx n Wx n Wx nWNx n Wx nWNx nx nWWkWWeekW 112200112220012202222222211（偶又奇/( ) ( )(), ,() ( )(), ,(NnkNnNnrNnjnrjnrnrnrnrNNNNNNX

19、kx nx nWkNkrNNXrx nx nWrWWeeW2 102 12022222220 222220 1122代入又）按按k k的奇偶将的奇偶将X(k)X(k)分成奇偶两部分分成奇偶两部分, k=2r, k=2r和和k=2r+1,k=2r+1,考虑考虑k k为偶数情况为偶数情况/( )( )(), ,()( )NnrNnNNg nx nx nnXrg n W2 1200 1 21222令令/( ) ( )(), ,() ( )(), ,NnkNnNnr nNnNX kx nx nWkNkrNNXrx nx nWr2 102 1200 22221210 1122代入考虑考虑k k为奇数情况

20、为奇数情况/( ) ( )(), ,()( )nNNnrNnNNh nx nx nWnXrh n W2 1200 1 212221令令结论结论 一个一个N点的点的DFT被分解为两个被分解为两个N/2点点;与时间抽取法的推与时间抽取法的推演过程一样，由于演过程一样，由于N=2M,因此因此,N/2仍为偶数，所以可以将仍为偶数，所以可以将N/2点点DFT的输出的输出X(k)再分为偶数组和奇数组，这样就将一个再分为偶数组和奇数组，这样就将一个N/2点的点的DFT分成两个分成两个N/4点点DFT的输入，也是将的输入，也是将N/2点的点的DFT的输入上、下对半分后通过蝶形运算而形成，直至最后的输入上、下对

21、半分后通过蝶形运算而形成，直至最后为为2点点DFT。/( )( )(), ,( ) ( )(), ,()( ),()( )nNNNnrnrNNnnNNg nx nx nnNNh nx nx nWnXrg n WXrh n W2 12 122000 1 21220 1 21222218点点DIF基基2FFT算法流图算法流图 3.3.2 频率抽取（频率抽取（DIF）基）基2FFT算法算法3.3.2 频率抽取（频率抽取（DIF）基）基2FFT算法算法DITDIT与与DIFDIF的相同之处：的相同之处：（1 1）DIFDIF与与DITDIT两种算法均为原位运算。两种算法均为原位运算。（2 2）DIFD

22、IF与与DITDIT运算量相同。运算量相同。DITDIT与与DIFDIF的不同之处：的不同之处：（1 1）DIFDIF与与DITDIT两种算法结构倒过来。两种算法结构倒过来。DIFDIF为输入顺序，输出乱序。运算完毕再运行为输入顺序，输出乱序。运算完毕再运行“二进二进制倒读制倒读”程序。程序。DITDIT为输入乱序，输出顺序。先运行为输入乱序，输出顺序。先运行“二进制倒读二进制倒读”程序，再进行求程序，再进行求DFTDFT。（2 2）DIFDIF与与DITDIT根本区别：在于蝶形结不同。根本区别：在于蝶形结不同。DITDIT的复数相乘出现在减法之前。的复数相乘出现在减法之前。DIFDIF的复数

23、相乘出现在减法之后。的复数相乘出现在减法之后。 以上讨论的都是以以上讨论的都是以2为基数的为基数的FFT算法，即算法，即N=2M，这，这种情况实际上使用得最多。种情况实际上使用得最多。 优点：程序简单，效率高，使用方便。优点：程序简单，效率高，使用方便。 实际应用时，有限长序列的长度实际应用时，有限长序列的长度N很大程度上由人很大程度上由人为因素确定，因此多数场合可取为因素确定，因此多数场合可取 N=2M，从而直接使用，从而直接使用以以 2 为基数的为基数的FFT算法。算法。 如如N不能人为确定不能人为确定 , N的数值也不是以的数值也不是以2为基数的整为基数的整数次方数次方,处理方法有两种处

24、理方法有两种: 补零补零:将将x(n)补零补零 , 使使N=2M. 例如例如N=30,补上补上x(30)=x(31)=0两点两点,使使N=32=25,这这样可直接采用以样可直接采用以2为基数为基数M=5的的FFT程序。程序。 有限长度序列补零后并不影响其频谱有限长度序列补零后并不影响其频谱 X(ej ) ，只是，只是频谱的采样点数增加了频谱的采样点数增加了,上例中由上例中由30点增加到点增加到32点，所点，所以在许多场合这种处理是可接受的。以在许多场合这种处理是可接受的。 如要求准确的如要求准确的N点点DFT值值,可采用任意数为基数的可采用任意数为基数的 FFT 算法算法 , 其其 计算效率低

25、于以计算效率低于以2为基数为基数FFT算法。算法。 如如 N 为复合数，可分解为两个整数为复合数，可分解为两个整数p与与q的乘积的乘积,像前面以像前面以2为基数时一样为基数时一样,FFT的基本思想是将的基本思想是将DFT的的运算尽量分小运算尽量分小,因此因此,在在N=pq情况下情况下,也希望将也希望将N点的点的DFT分解为分解为p个个q点点DFT或或q个个p点点DFT,以减少计算量。以减少计算量。步骤：步骤：1010nnQnkk Pk01,n k分别为0, 1,，Q-1; 10,n k分别为0, 1,，P-1。 N点点DFT可以重新写成为可以重新写成为 101010( )(),( )NknNn

26、X kX k PkX k kx n W10100111()()1000QPk P kn Q nNnnx nQn W0 11 00 01 1010 11 00 00111101000111000,QPk n Qk n Pk nk n PQNNNNnnQPk n Qk n Pk nNNNnnX k kx n n WWWWx n n WWW考虑到 再令1010nkPQnkNWW令 1010010100100110,QnnkQnkNPnnkPWWWnnxkkXPnnkPWnnxnkX,01000,0011001nkQnkNQnWWnkXkkX00001001,nkNWnkXnkXQnnkQWnkXkk

27、X,(1) 先将先将 x(n)通过通过 x(n1Q+n0)改写成改写成 x(n1,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)以以P=3，Q=4, N=12为例为例 (2) 求个点的求个点的DFT 20301001110,nnkWnnxn

28、kX（）（）X1(k0,n0)乘以得到乘以得到X1(k0,n0)。（）求（）求P个个点的点的DFT，参变量是，参变量是k0（）将（）将X2(k0,k1)通过通过X(k0+k1P)恢复为恢复为X(k)N=12为组合数时的为组合数时的FFT（）求（）求个个P点点DFT需要需要P2次复数乘法和次复数乘法和P(P-1)次复数加法；次复数加法；（）乘（）乘个个W因子需要因子需要次复数乘法；次复数乘法；（）求（）求P个个点点DFT需要需要PQ2 次复数乘法和次复数乘法和P（-1)次复数加法。次复数加法。总的复数乘法量总的复数乘法量: QP2+N+PQ2=N(P+Q+1)；总的复数加法量总的复数加法量: Q

29、P(P-1)+PQ(Q-1)=N(P+Q-2) 例例: N=23*29=667, N2=444889, N(P+Q+1)=35351当组合数当组合数 N N= =P P1 1P P2 2P P3 3P Pm m中所有的中所有的P Pi i均为均为4 4时，就是时，就是基四基四FFTFFT算法算法 以上提出以上提出FFT算法，可以很快地求出全部算法，可以很快地求出全部DFT值。值。即求出有限长序列即求出有限长序列x(n)的的z变换变换X(z)在单位园上在单位园上N个等个等间隔抽样点间隔抽样点zk处的抽样值。它要求处的抽样值。它要求N为高度复合数。为高度复合数。即即N可以分解成一些因子的乘积。例可

30、以分解成一些因子的乘积。例N=2L 实际上：实际上：（1）也许对其它围线上也许对其它围线上z变换取样发生兴趣变换取样发生兴趣。（2）只需要计算单位圆上某一段的频谱只需要计算单位圆上某一段的频谱,即即M不等于不等于N。如窄带信号，希望在窄带频率内频率抽样能够非常如窄带信号，希望在窄带频率内频率抽样能够非常密集，提高分辨率，带外则不考虑。密集，提高分辨率，带外则不考虑。（3）若若N是大素数时，不能加以分解，又如何有效计是大素数时，不能加以分解，又如何有效计算这种序列算这种序列DFT。例。例N=311，若用基，若用基2则须补则须补N=28=512点，要补点，要补211个零点。个零点。问题提出问题提出

31、 为了提高为了提高DFT的灵活性，须用新的方法。线性调的灵活性，须用新的方法。线性调频频z变换变换(CZT)就是适用这种更为一般情况下，由就是适用这种更为一般情况下，由x(n)求求X(zk)的快速变换的快速变换。CZT 来自于雷达专业的专用词汇。来自于雷达专业的专用词汇。Z 变换采用螺线抽变换采用螺线抽样样,可计算单位圆上任一段曲线的可计算单位圆上任一段曲线的Z变换，适用于变换，适用于更一般情况下（更一般情况下（M不等于不等于N）由）由x(n)求求X(zr)的快速的快速算法算法, 达到频域细化的目的，这种变换称为线性调达到频域细化的目的，这种变换称为线性调频频Z变换变换(简称简称CZT )。

32、为适应为适应z可以沿平面内更一般的路径取值，可以沿平面内更一般的路径取值，我们我们沿沿z平面上的一段螺线作平面上的一段螺线作等分角的抽样等分角的抽样，则，则z的取样点的取样点Zr可表示为：可表示为：)( )NnnX zx n z10（, ,rrzAWrM0 11,jjAA eWW e0000 已已 知知 N点序列点序列x(n) ，0nN-1,其其z变换为变换为其中其中M：表示欲分析的复频谱的点数。：表示欲分析的复频谱的点数。M不一定等于不一定等于N。A，W 都为任意复数都为任意复数,令令 一、一、CZTCZT的定义的定义( ) ( )( )( )NnkrnNnnrnzX zCZT x nx n

33、 zx n A W1010代入中得, ,rM0 11上式即为上式即为CZTCZT的定义的定义. .现在讨论现在讨论A A0 0,W,W0 0, ,0 0, ,0 0的的含义含义: : 为输出为输出M M点的变换域值点的变换域值.r=0.r=0时的时的A A0 0e ej j0 0是是CZTCZT的起点的起点, ,随着随着r r的变化的变化,r,r0 0,r,r1 1, ,r,rM-1M-1构成构成CZTCZT的变化路径，对于的变化路径，对于M-1M-1点其极坐标为点其极坐标为()()jj MMQA eWe001100CZTCZT在现在现Z Z平面上的变换路径是一条平面上的变换路径是一条螺旋线螺

34、旋线。（1）A为起始样点位置为起始样点位置,jAA eA0000：起点半径，大于起点半径，大于1 1时，表示螺旋线在单位圆外时，表示螺旋线在单位圆外，反之，在单位圆内。，反之，在单位圆内。起点半相角。起点半相角。（2）当当W01,螺旋线内旋，反之外旋。螺旋线内旋，反之外旋。（3）当当A0=W0=1时，时，CZT的变换路径为单位圆上的变换路径为单位圆上的一段弧，起点为的一段弧，起点为P，终点为，终点为Q，且，且M不一定等不一定等于于N。变换求出该序列即由，为均匀分布在单位园上此时等分）时，（。即即（当满足下面特殊条件：DFTCZTzNNMWeeWWcAeAAbNMarNjjj221,)(01,

35、1)(: )4(002000000()()() ()( )( ) ( )nr nrNNnnrnrnnrnr nNnnnrnrrnX zx n A Wx n A WWWWx n A WW22222222211222001222012利用 我们希望做的是频谱分析，因此考虑我们希望做的是频谱分析，因此考虑A0=W0=1时，时，在单位圆上在单位圆上CZT，且，且M不一定等于不一定等于N。( )( ), ( )()( ) (), ,()( )( )() ( )( ), ,( )nnnrNrnrrrrrg nx n A Wh nWX zWg n h rn rMX zg nh nWX zWg rh rkMW

36、y r222222221202220 110 11令由上式可知：可以由与的离散卷并乘得到，即：()rX z( )x nnnA W22rW22( )nh nW22( )g n( )y rCZTCZT的线性滤波计算步骤的线性滤波计算步骤二、二、CZTCZT的计算方法的计算方法分析：分析：从上面的推导过程可以看出，计算从上面的推导过程可以看出，计算CZTCZT关键是计算一个线性卷积关键是计算一个线性卷积( )( )( )y rg rh r其中，其中，g(n)g(n)应为应为N N点序列，点序列，h(n)h(n)应为偶对称的应为偶对称的无限长序列，考虑到输出无限长序列，考虑到输出M M点序列，点序列，

37、h(n)h(n)的实的实际长度应为际长度应为M M点。因此，可用点。因此，可用DFTDFT来实现两者来实现两者的卷积，具体步骤如下：的卷积，具体步骤如下：()( ),( )( ),.jjnnnnmAA eWW eAWA WnNg nnNnNAx ng nNnLLNML0022000000220101011012根据已知和即由已知 、求出的各值, 并得且使（1 1）计算并设置）计算并设置g(n)g(n)（2 2）计算并设置）计算并设置h(n)h(n)将将h(n)h(n)设置成长度为设置成长度为L L的序列，考虑到其为偶对的序列，考虑到其为偶对称序列，且取称序列，且取M M点（输出点（输出M M点

38、序列），取点序列），取()( )nL nWnMh nMnLNLNnLW 222201011如下图所示如下图所示（3 3）计算）计算h h(n)(n)和和g g(n)(n)的的DFTDFT，得到，得到L L点序列点序列H(k)和和G(k)。（4 4）令）令Y(k) =H(k) G(k)（乘积）后作乘积）后作Y(k) 的的IDFTIDFT（反变换）得到时域输出序列反变换）得到时域输出序列y(r) 。（5 5）取）取y(r)y(r)的前的前M M点，并乘以点，并乘以W W-r2/2-r2/2, ,则得最后的输出则得最后的输出X(zX(zr r) )，即，即()( )rrX zWg r22 与标准与标

39、准FFT算法相比，算法相比，CZT算法有以下特点：算法有以下特点：（1）输入序列长）输入序列长度度N及输出序列长及输出序列长度度M不需要相等不需要相等，且且N及及M不必是高度合成数，二者均可为素数。不必是高度合成数，二者均可为素数。（2）Zk的角间隔是任意的，的角间隔是任意的，说明说明其频率分辨率也是其频率分辨率也是任意任意可控可控的的，角间隔小，分辨率高，反之，分辨，角间隔小，分辨率高，反之，分辨率低。率低。（3）周线不必是）周线不必是z平面上的圆，在语音分析中螺旋平面上的圆，在语音分析中螺旋周线具有某些优点。周线具有某些优点。（4）由于由于起始点起始点z0可任意选定，因此可以从任意频可任意

40、选定，因此可以从任意频率上开始对输入数据进行窄带高分辨率的分析。率上开始对输入数据进行窄带高分辨率的分析。总之，总之，CZT算法具有很大的灵活性算法具有很大的灵活性。nkNW1)IDFT1)IDFT的运算方法的运算方法101( )IDFT( )( )NnkNkx nX kX k WN10()D FT ( )( )NnkNnXkx nx n WIDFT： DFT：以上所讨论的以上所讨论的FFTFFT算法可用于算法可用于IDFTIDFT运算运算简称为简称为IFFTIFFT比较比较IDFTIDFT的定义式：的定义式：IDFTIDFT与与DFTDFT的差别：的差别： 1 1）把）把DFTDFT中的每一

41、个系数中的每一个系数 改为改为 ， 2 2）再乘以常数）再乘以常数 1/N 1/N ，nkNW10)(1)(NknkNWkXNnx)(1)(1)(10kXDFTNWkXNnxNknkN)(kX 第二种方法，完全不需要改动第二种方法，完全不需要改动FFT程序，而是直接利用它作程序，而是直接利用它作 IFFT。 考虑到考虑到 故故 IFFT计算分三步：计算分三步： 将将X（k）取共轭（虚部乘以）取共轭（虚部乘以-1） 对对 直接作直接作FFT 对对FFT的结果取共轭并乘以的结果取共轭并乘以1/N，得，得x（n）。）。2）实数序列的）实数序列的FFT 以上讨论的以上讨论的FFT算法都是复数运算算法都

42、是复数运算,包括序列包括序列x(n)也认为是复数也认为是复数,但大多数场合但大多数场合,信号是实数序列信号是实数序列,任何实任何实数都可看成虚部为零的复数数都可看成虚部为零的复数,例如例如,求某实信号求某实信号x(n)的复的复谱谱,可认为是将实信号加上数值为零的虚部变成复信号可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用再用FFT求其离散傅里叶变换。这种作法很求其离散傅里叶变换。这种作法很不经济不经济,因为把实序列变成复序列因为把实序列变成复序列,存储器要增加一倍存储器要增加一倍,且计算机运行时且计算机运行时,即使虚部为零即使虚部为零,也要进行涉及虚部的也要进行涉及虚部的

43、运算运算,浪费了运算量。合理的解决方法是利用复数据浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算对实数据进行有效计算,下面介绍两种方法。下面介绍两种方法。 （1）用）用 一个一个N点点FFT同时计算两个同时计算两个N点实序列的点实序列的DFT 设设x (n)、y (n)是彼此独立的两个是彼此独立的两个N点实序列点实序列,且且 X (k)=DFTx (n) ， Y (k)=DFTy(n) 则则X (k)、Y(k)可通过一次可通过一次FFT运算同时获得。运算同时获得。 首先将首先将x (n)、y(n)分别当作一复序列的实部及虚部分别当作一复序列的实部及虚部,令令 g(n)=x

44、(n)+jy(n)通过通过FFT运算可获得运算可获得g(n)的的DFT值值 kjGkGkjYkjXkYkXkYkjYkjXkXkjYkXkGirriiririr kNGkGjkNGkGkNGkGkXngiirr212121ReDFT*利用离散傅里叶变换的共轭对称性利用离散傅里叶变换的共轭对称性通过通过g(n)的的FFT运算结果运算结果G(k),由上式也可得到由上式也可得到Y(k) 的值。的值。 kNGkGjkNGkGkNGkGkjYngjiirr212121ImDFT* kNGkGjkNGkGkYrrii2121 （2）用一个）用一个N点的点的FFT运算获得一个运算获得一个2N点实序列的点实序

45、列的DFT 设x(n)是2N点的实序列,现人为地将x(n)分为偶数组x1 (n)和奇数组x2 (n) x1 (n)= x(2n) n=0,1,N-1 x2 (n)= x(2n+1) n=0,1,N-1然后将x1 (n) 及x2 (n) 组成一个复序列： y(n)= x1 (n) +j x2 (n) 通过N点FFT运算可得到：Y(k)=X1(k)+jX2(k) ，N点根据前面的讨论，得到)()(2)()()(21)(*2*1kNYkYjkXkNYkYkX 为求为求 2N 点点 x(n)所对应所对应 X(k) ，需求出，需求出 X(k)与与 X1 (k) 、 X2 (k) 的的关系关系 ： 10)

46、12(21022120) 12()2()()(NnknNNnnkNNnnkWnxWnxWnxkX10210) 12()2(NnnkNkNNnnkNWnxWWnx 1011022101011) 12()()()2()()(NnnkNNnnkNNnNnnkNnkNWnxWnxkXWnxWnxkX而1）由）由x1(n)及及x2(n)组成复序列，经组成复序列，经FFT运算求得运算求得 Y(k)，2）利用共轭对称性求出）利用共轭对称性求出 X1(k)、X2(k)，3）最后利用上式求出）最后利用上式求出 X（k）, 达到用一个达到用一个N点的点的FFT计计算一个算一个2N点的实序列的点的实序列的DFT的目

47、的。的目的。 X(k)=X1(k)+W2Nk X2(k)所以所以 3) 线性卷积的线性卷积的FFT算法算法 线性卷积是求离散系统响应的主要方法之一线性卷积是求离散系统响应的主要方法之一,许多许多重要应用都建立在这一理论基础上重要应用都建立在这一理论基础上,如卷积滤波等。如卷积滤波等。 以前曾讨论了用圆周卷积计算线性卷积的方法归纳以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下如下: 将长为将长为N2的序列的序列x(n)延长到延长到L,补补L-N2个零，个零， 将长为将长为N1的序列的序列h(n)延长到延长到L,补补L-N1个零，个零， 如果如果LN1+N2-1,则圆周卷积与线性卷积相等则圆周卷

48、积与线性卷积相等,此时此时,可用可用FFT计算线性卷积，方法如下计算线性卷积，方法如下: a.计算计算X(k)=FFTx(n)b. 求求H(k)=FFTh(n)c. 求求Y(k)=H(k)X(k) k=0-L-1d. 求求y(n)=IFFTY(k) n=0-L-1 可见，只要进行二次可见，只要进行二次FFT,一次一次IFFT就可完成就可完成线性卷积计算。线性卷积计算。 计算表明计算表明,L32时时,上述计算线性卷积的方法比上述计算线性卷积的方法比直接计算线卷积有明显的优越性直接计算线卷积有明显的优越性,因此因此,也称上述循也称上述循环卷积方法为快速卷积法。环卷积方法为快速卷积法。 上述结论适用

49、于上述结论适用于 x(n)、h(n) 两序列长度比较接近两序列长度比较接近或相等的情况或相等的情况,如果如果x(n)、h(n) 长度相差较多长度相差较多,例如例如, h(n) 为某滤波器的单位脉冲响应为某滤波器的单位脉冲响应,长度有限长度有限,用来处理一个用来处理一个很长的输入信号很长的输入信号 x(n)，或者处理一个连续不断的信号，或者处理一个连续不断的信号，按上述方法按上述方法,有三个问题：有三个问题：（1） h(n) 要补许多零再进行计算要补许多零再进行计算,计算量有很大的浪计算量有很大的浪费，或者根本不能实现。费，或者根本不能实现。（2）系统的存储量要求极高。）系统的存储量要求极高。（

50、3）带来了很大的系统延迟。）带来了很大的系统延迟。 为了克服上述三个问题，保持快速卷积法的优越为了克服上述三个问题，保持快速卷积法的优越性性,可将可将 x(n) 分为许多段分为许多段,每段的长度与每段的长度与 h(n) 接近接近 , 处处理方法有两种：理方法有两种：h(n)x(n)01) 1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnhnxny)()(*)()(*)()()(*)()(nhnxnyii则输入序列可表为：于是输出可分解为： 其中 假定 xi(n) 表示 x（n）序列的第i段 ： 1）先对 h（n）及 xi（n）补零，补到具有N点长度，N=N1+N2

51、-1。 一般选 N=2M。 iinyny)()(由于 yi（n）的长度为N1，而xi（n）的长度为N2，因此相邻两段 yi（n）序列必然有N-N2=N1-1点发生重叠。2）用基2 FFT计算 yi（n）=xi（n）*h（n）。3）重叠部分相加构成最后的输出序列。计算步骤：计算步骤：a. 事先准备好滤波器参数事先准备好滤波器参数 H（k）=DFTh(n)，N点点b. 用用N点点FFT计算计算Xi（k）=DFTxi(n)c. Yi(k)=Xi(k)H(k)d. 用用N点点IFFT求求 yi（n）=IDFTYi(k)e. 将重叠部分相加将重叠部分相加 这种方法和第一种方法稍有不同，即将上面分段序列中

52、补零的部分不是补零，而是保留原来的输入序列值，这时，如利用DFT实现h(n)和xi(n)的循环卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。 重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。 y0(n)中的中的N1 1, L 1点对应于线性卷积点对应于线性卷积 x(n) h(n)中的中的0 , N2-1点点y1(n)中的中的N1 1, L 1点对应于线性卷积点对应于线性卷积 x(n) h(n)中的中的N2,2N2-1点点y2(n)中的中的N1 1, L 1点对应于线性卷积点对应于线性卷积 x(n) h(n)中的中的2N2,3N2-1点点依此类推依此类推 ，

53、并将，并将yi(n)拼接起来构成拼接起来构成y (n) 4）用）用FFT计算相关函数计算相关函数 相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。 两个长为N的实离散时间序列 x（n）与y（n）的互相关函数定义为 1010)()()()()(NnNnxynynxnynxr则可以证明，则可以证明，rxy()的离散傅里叶变换为的离散傅里叶变换为 Rxy（k）=X*(k)Y(k) 其中其中 X(k)=DFTx(n), Y(k)=DFTy(n), Rxy(k)=DFTrxy() ， 0kN-1互相关函数定义为 x(n)及y(n)的卷积公式 相比

54、较，我们可以得到相关和卷积的时域关系： 1010)()()()()(NnNnxymnynxnymnxmr)()()()()(10mymxnynmxmfNn)()()()()()()(1010mymxnynmxnymnxmrNnNnxy10102)()(1)()()(NnNkkNjxyekYkXNnynxr因 得 证毕。)(*)()(DFTkXnRnxNN 当x（n）=y（n）时，得到x（n）的自相关函数为： 维纳辛钦定理： 自相关函数与信号功率谱互为傅立叶变换对。 10)()()(NnxxnxnxrkNjNkekXN2102| )(|1 上面的推导表明，互相关和自相关函数的计算可利用FFT实现

55、。 由于离散傅里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。 直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。 实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与循环卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。11)(mNmrxy10)(Nmmrxy(5) 对R（k）作 IFFT，取后N-1项，得取前N项，得 （1） 设x(n)和y(n)的长均为N，求线性相关；（2）为了使两个有限长序列的线性相关可用其循环相关代替而不产生混淆，选择周期 L2N-1，且L=2m，以使用FFT，将x(n)，y(n)补零至长为L。 （3） 用FFT计算X(k)，Y(k)（k=0，1，L-1）（4） R（k）=X*（k）y（k）幅度例8-20-1001020-50050100150200250300m-20-1001020-50050100150200250300m幅度幅度(a) (b)延迟序列的互相关函数(a)和自相关函数(b) 5）用）用FFT计算二维离散的傅里叶变换计算二维离散的傅里叶变换 二维信号有图象信号、时空信号、时频信号等。二维离散傅里叶变换可用于处