数字信号处理4人文科技ppt课件

1、 第第4 4章章 FFT FFT 4.1 4.1 引言引言 4.1.1 4.1.1 离散傅里叶变换的矩阵表示及其运算量离散傅里叶变换的矩阵表示及其运算量 DFT DFT在数字信号处置中起着非常重要的作用，在数字信号处置中起着非常重要的作用， 这是与这是与DFTDFT存在着高效算法，存在着高效算法， 即快速傅里叶变换即快速傅里叶变换(FFT) (FFT) 分不开的。快速运算的关键是减少运算量。分不开的。快速运算的关键是减少运算量。n离散傅里叶变换对为：n (4.1)n (4.2)n式中 。 下面要用矩阵来表示DFT关系。令： 101, 1 , 0,)()(:NnnkNNkWnxkXDFT1, 1

2、 , 0,)(1)(:10NnWkXNnxIDFTNknkNNjNeW2) 1() 1 ()0(Nxxxx) 1() 1 ()0(NXXXX并令TN、 表示两个变换方阵，有：1NT2)1(2)1(1)1(22221211111111NNNNNNNNNNNNNNNWWWWWWWWWT2)1(2)1()1()1(2222)1(2111111111nNNNNNNNNNNNNNNWWWWWWWWWTn假设用i (0iN-1)表示这两个N阶方阵的行号，用j(0jN-1)表示这两个N阶方阵的列号，那末很容易看出，TN 方阵中i行j列的元素为 ，而方阵中i行j列的元素为 。n 于是4.1式可以写成： n而4

3、.2式可以写成： jiNWjiNWxTXN 11XTNxNn 普通情况下，信号序列x(n) 及其频谱序列X(k) 都是用复数来表示的，WN当然也是复数。因此，计算DFT的一个值X(k) 需求进展N次复数乘法与1相乘也包括在内和N-1次复数加法；所以，直接计算N点的DFT需求进展N2 次复数乘法和N(N-1) 复数加法。 n显然，直接计算N点的IDFT所需的复乘和复加的次数也是这么多。当N足够大时，N2 N(N-1), 因此，DFT与IDFT的运算次数与N2 成正比，随着N的添加，运算量将急剧添加，而在实践问题中，N往往是较大的，因此有必要对DFT与IDFT的计算方法予以改良。 4.1.2 因子

4、的特性因子的特性 DFT和和IDFT的快速算法的导出主要是根据的快速算法的导出主要是根据 因子的特性。因子的特性。 1周期性：周期性： 对离散变量对离散变量n有同样的周期性。有同样的周期性。nkNnNNnkNNknNWWWW)(nkNWnkNW 2对称性：对称性： 或或 3. 其它其它: knNNknNnkNnkNWWWW)()(*nkNNnkNknNnkNWWWW)()(*kNNNjkNNNkNNkNWeWWWW222)2(kNkNjkNjkNWeeW 2222224.2 4.2 基基2 2时间抽选的时间抽选的FFTFFT算法算法 4.2.1 4.2.1 算法推导算法推导 曾经知道：曾经知道

5、： 令令DFTDFT的长度的长度N=2MN=2M，M M为正整数。为正整数。1-N,0,1,k )()()(10nkNNnWnxnxDFTkXn令： n于是有：12N,0,1,r ) 12()()2()(rxrprxrg1-N,0,1,k )()()()() 12()2()(120rk 2/120 2/120)12(1202kPWkGWrpWWrgWrxWrxkXkNNrNkNNrrkNNrkrNNrrkNn其中， n是由x(n)的偶数抽样点构成的DFT；而n kr 2/120120kr 2/)2()()(NNrNrNWrxWrgkG120kr 2/kr 2/120) 12()()(NrNNN

6、rWrxWnpkPn是由x(n)的奇数抽样点构成的DFT。但是这两个式子并不完全是N/2点的DFT，由于k的范围依然是由0到N-1，因此，还应该进一步思索k由N/2到N-1范围的情况。n如今令 ，故对于后半段有：n n同理： n又知： 12, 1 , 0Nk)()()()()2(120kr 21202120)2(2/22kGWrgWWrgWrgNkGNrNNrrkNrNrNkrNNN)()2(kPNkPk 2NNkNWW 图图 4.1 N点点DFT分解为两个分解为两个N/2点的点的DFT (N=8) 图图 4.2 N/2 点的点的DFT 分解为两个分解为两个N/4点的点的DFT (N=8)n综

7、上所述，可以得到： n其中G(k)、P(k) 分别是x(n)的偶数点和奇数点的N/2点DFT。 )()()2()()()(kPWkGNkXkPWkGkX kNkN12, 1 , 0Nkn 这样，我们就将一个N点的DFT分解成了两个N/2点的DFT，由于DFT的运算量与其点数的平方成正比，因此使运算量减少了。但是，还应该将每一个N/2点的DFT再分解为两个N/4点的DFT，如此下去，直到分解为2点的DFT为止，总共需求进展log2N-1=log2(N/2)次分解。 图图 4.3 2 点点 DFT 信号流图蝶形图信号流图蝶形图n对于2点DFT，有： n所以2点DFT的运算只需一次加法和一次减法，这

8、样的运算叫做蝶形运算，这样的信号流图叫做蝶形图。1111111122WTn该算法每次分解都是将时域序列按奇偶分为两组，因此要求N等于2的正整数幂，故将这种FFT算法叫做基2时间抽选法。 4.2.2 算法特点 1. 倒序重排这种FFT算法的每次分解都是将输入序列按照奇偶分为两组，故要不断地将每组输入数据按奇偶重排，直到最后分解为2点的DFT，输入数据才不再改动顺序。这样做的结果，使得作FFT运算时，输入序列的次序要按其序号的倒序进展重新陈列。 n如今将图4.4中输入序号以及重排后的序号按二进制写出如下注：下标“2表示二进制数。可以看出，将输入序号的二进制表示n2n1n0位置颠倒，得到n0n1n2

9、，就是相应的倒序的二进制序号。因此，输入序列按倒序重排，实践上就是将序号为n2n1n0的元素与序号为n0n1n2的元素的位置相互交换。 2. 同址计算同址计算 从图从图4.4可以看出，整个算法流图可以分为四段，可以看出，整个算法流图可以分为四段，0段段为倒序重排，后面为倒序重排，后面3段为段为3(log28=3)次迭代运算：首先计算次迭代运算：首先计算2点点DFT，然后将，然后将2点点DFT的结果组合成的结果组合成4点点DFT，最后将，最后将4点点DFT组合为组合为8点点DFT。因此，对于。因此，对于N点点FFT，只需求一列，只需求一列存储存储N个复数的存储器。个复数的存储器。 n开场时，输入

10、序列的N个数放于此N个存储器内，倒序重排后仍存于这N个存储器中，每一次迭代运算后的结果也依然存于这N个存储器中，整个运算完成后，这N个存储器中即为所求的频谱序列X(k) k = 0、1、.、 N-1。这就是所谓的同址计算，这样可以大大节约存储元件。 3. 运算量运算量察看图察看图4.4可知，图可知，图4.3所示的蝶形图实践上代表了所示的蝶形图实践上代表了FFT的根的根本运算，它实践上只包含了两次复数加法运算。一个本运算，它实践上只包含了两次复数加法运算。一个N(N=2M) 点的点的FFT，需求进展，需求进展M=log2N次迭代运算。每次迭代运算。每次迭代运算包含了次迭代运算包含了N/2个蝶型，

11、因此共有个蝶型，因此共有N次复数加法；次复数加法；此外，除了第一次的此外，除了第一次的2点点DFT之外，每次迭代还包括了之外，每次迭代还包括了N/2次复数乘法即乘次复数乘法即乘WN的幂。因此，一个的幂。因此，一个N点的点的FFT共共有复数乘法的次数为：有复数乘法的次数为： 复数加法的次数为： 因此，FFT算法比直接计算DFT大大减少了运算量，尤其是当N较大时，计算量的减少更为显著。比如，当N=1024时，采用FFT算法时复数乘法的次数，不超越直接计算DFT时复乘次数的千分之五。 2log2) 1(log2222NNNNMcNNNNAc222loglog22 4.2.3 关于FFT算法的计算机程

12、序在普通情况下，进展FFT运算的序列至少都有几百点的长度，因此需求编制FFT算法程序以便可以利用计算机来快速进展计算。可以参照图4.4来编制计算N (N必需等于2的正整数幂)点FFT的计算机程序，整个程序可以分为两部分：一部分是倒序重排，另一部分是用三层嵌套的循环来完成M=log2N次迭代。 n三层循环的功能是：最里的一层循环完成蝶形运算，中间的一层循环完成因子WNk的变化，而最外的一层循环那么是完成M次迭代过程。n倒序重排的程序是一段经典程序，它以巧妙的构思、简单的语句用高级编程言语来完成了倒序重排的功能。下面是一段用FORTRAN言语编写的倒序重排程序。 N=2*M 表示N=2M ，M是输

13、入的正整数 NV2=N/2 NV2是一个整数变量 NM1=N-1 NM1也是一个整数变量 J=1 对变量J赋初值 100 DO 7 I=1,NM1 (循环开场，到语句7终了； 循环变量I从1开场，到NM1终了，步长为1)n IF (I.GE.J) GOTO 5 假设IJ，就转移到语句5n 将输入序列中序号互为倒序n 的两个数值交换位置 TIXIXJXJXT)( )()()( 5 K=NV2 6 IF(K.GE.J) GOTO 7 (假设KJ，就转移到语句 7)J=J-kK=K/2 GOTO 6 转移到语句67 J=J+K n由于是同址计算，因此程序中只用了一个数组X( )，这个数组共有N个元素

14、。X( )既表示输入序列，也表示按倒序重排后的这N个数值。程序中的字母都是大写的，这是FORTRAN语句的特点。nI既是循环变量，也代表输入序列的正序序号，J代表倒置后的序号。由于FORTRAN言语的循环变量不能从0开场，故以8点FFT为例，I的范围为18，下面是正序序号I与倒序序号J之间的对应关系：n I: 1 2 3 4 5 6 7 8 n J: 1 5 3 7 2 6 4 8n 输入序列按倒序重排，实践上就是将X(I)和X(J)互换位置。显然，假设I=J，就不需求交换；而假设IJ，就表示X(I)与X(J)曾经交换过了。因此，在循环语句下面的语句“IF (I.GE.J) GOTO 5就阐明

15、了交换的条件：只需当IJ时才执行下面三条语句，使X(I)与X(J)互换位置。n正序序号I也是循环变量，从1开场，每次循环添加1。关键问题是对于每一个I，所对应的J是什么？程序中从语句5开场直到循环终了就是为下一次循环确定所对应的J的。虽然I-1和J-1的二进制表示是互为倒序的关系，但是，程序中并不是由I得到J，而是由上一个J来得到下一个J。思绪是：I-1由0到7，用二进制来表示就是从0002开场，每次在最低位加1，逐次添加到1112 。 n既然J-1的二进制表示是I-1二进制表示的倒置，那末J-1的变化就应该是在最高位逐次加1，而最高位的二进制数1表示N/2。所以，程序中将J与N/2比较来判别

16、J-1的最高位是0还是1，假设是0，就执行J+K=J+N/2，也就是将J-1的二进制最高位由0变为1。假设J-1的最高位是1，应该将它变为0，也就是将J减去N/2，然后调查次高位是0还是1，这应该用N/4来调查，于是执行“K=K/2。如此下去，直到得到下一个J，完成此次循环。4.3 4.3 基基2 2频率抽选的频率抽选的FFTFFT算法算法 时间抽选法是在时域内将输入序列时间抽选法是在时域内将输入序列x(n)x(n)逐次分解为偶数逐次分解为偶数点子序列和奇数点子序列，经过求子序列的点子序列和奇数点子序列，经过求子序列的DFTDFT而实现整而实现整个序列的个序列的DFTDFT。而频率抽选法那么是

17、在频域内将。而频率抽选法那么是在频域内将X(k)X(k)逐次逐次分解成偶数点子序列和奇数点子序列，然后对这些分解得分解成偶数点子序列和奇数点子序列，然后对这些分解得越来越短的子序列进展越来越短的子序列进展DFTDFT运算，从而求得整个的运算，从而求得整个的DFTDFT。当。当然，同样要求然，同样要求N N为为2 2的正整数幂。的正整数幂。 n 设 ，那么可以分别表示出k为偶数和奇数时的X(k)。n n其中， 12, 1 , 0Nr120rn 2/120rn 2/120rn 2/120120)2(22102)()2()()2()()()2(NnNNnrNNNNnNNnNnrNnNnrNNnrnN

18、WngWWNnxWnxWNnxWnxWnxrX11)2()()(2N,0,n Nnxnxng120rn 2/120rn 2/120rn 2/12022120rn 2/120120)2)(12(2)() 1()2()()2()()2()() 12 (NnNnNNnnNNNnnNNNnNNrNNnNnrNNnnNNNnNnNnrNnNnrNWWnpWWNnxWWnxWWWWNnxWWnxWNnxWWnxrXn其中， n 显然，X(2r) 为g(n) 的N/2点DFT，X(2r+1) 为p(n)WNn 的N/2点DFT。这样，一个N点的DFT分解为两个N/2点的DFT。将分解继续下去，直到分解为2点

19、的DFT为止。当N=8时，基2频率抽选的FFT算法的整个信号流图如图4.6所示。11)2()()(2N,0,n Nnxnxnpn将图4.6与图4.4比较，可知频率抽选法的计算量与时间抽选法一样，而且都可以同址计算。时间抽选法是输入序列按奇偶分组，故x(n)的顺序要按倒序重排，而输出序列按前后分半，故X(k) 的顺序不需求重排；频率抽选法那么是输出序列按奇偶分组，故X(k) 的顺序要按倒序重排，而输入序列按前后分半，故x(n) 不需求重排。4.4 4.4 基基 4 4时间抽选的时间抽选的FFTFFT算法算法 上面讨论的上面讨论的FFTFFT算法都要求算法都要求N=2M N=2M ，M M是正整数

20、，因此是是正整数，因此是基基2 2的的FFTFFT算法。假设算法。假设N=4M N=4M ，那么还可以采用基，那么还可以采用基4 4的的FFTFFT算算法。与推导基法。与推导基2 FFT2 FFT算法类似，可以推导基算法类似，可以推导基4 FFT4 FFT算法。算法。将将x(n) x(n) 相间地分为四组，所得的相间地分为四组，所得的X(k) X(k) 按前后分为四组。第按前后分为四组。第一次分组的流图如图一次分组的流图如图4.74.7所示。其中，中间的矩阵为变换所示。其中，中间的矩阵为变换方阵方阵T4T4， j- 1- j 11- 1 1- 1j 1- j- 11 1 1 111111119

21、464346444243424144WWWWWWWWWT而 Di (i = 0,1,2,3) 那么为N/4 阶的对角矩阵，即： 0,1,2,3i WWWDiNNiNiNi )14(21 即变换方阵，即： 4/NT2)14N(N/42)14N(N/41)14N(N/4)14N2(N/422N/412 N/414NN/42 N/41 N/44 W W W1 W W W1 W W W11 1 1 1NT 图图 4.7 基基 4 时间抽选的时间抽选的 FFT算法的第一次分解算法的第一次分解n为了了解基4 FFT算法，可以将基2时间抽选的FFT算法的第一次分解与基4的第一次分解进展比较。n )()()(

22、)(1- 11 1)2()(2kPWkGTkPWkGNkXkXkNkN12, 1 , 0Nkn基2的情形第一次分解为两组，基4为四组。基2分解中的G(k) 和P(k) 分别为x(n)的偶数点和奇数点的N/2点DFT，而基4分解的每一项中的因子n ) () () (4/xxxTNn都代表一个N/4点的DFT，x(n) 的间隔为4。基2中的前后两项用T2矩阵来衔接；而在基4的情形为4项，是用T4矩阵来衔接的。基2中的第二项的N/2点的DFT之前要乘上因子WNk；而基4的情况每个N/4点的DFT之前那么要乘上对角矩阵Di，其对角线上的元素为WNki。因此，基4的FFT分解既与基2 FFT的分解相对应

23、，又比之复杂。n基4的情形也要继续分解下去，直到分解为四点DFT为止，共要经过log4N-1= log4(N/4) 次分解。估算基4 FFT的计算量：一个N点的FFT，要经过log4N次变换。复数乘法表达在与对角矩阵D1、D2、D3相乘，总的复乘次数为： ；总的复加次数为： 。 ) 1(log4344NNMcNNAc44log3n 基2 FFT算法的复乘次数为：n n复加次数为： )21(log2log222log24422NNNNNNMcNNNNAc422log2logn将Mc4与Mc2比较，Ac4和Ac2比较，可知，基4的复乘次数约减少了1/4，但复加次数却添加了。因此，在乘法运算要转化为

24、加法运算来实现的情况下，基4算法的计算量会比基2算法少；但是，近年来，在许多微处置器和DSP中，实现一次乘法运算与实现一次加法运算的时间完全一样，这时采用基4算法就不合算了，而且基4算法还比基2算法复杂。 4.5 4.5 快速傅里叶反变换快速傅里叶反变换IFFTIFFTIFFTIFFT是是IDFTIDFT的快速算法。由于的快速算法。由于DFTDFT的正变换和反变换的表达的正变换和反变换的表达式类似，因此式类似，因此IDFTIDFT也有类似的快速算法。可以用三种不同也有类似的快速算法。可以用三种不同的方法来导出的方法来导出IFFTIFFT算法。算法。方法方法1 1 设设 ， 那么那么有：有： ，

25、 n=0,1,N-1 n=0,1,N-1)()(kXnxDFT 10)(1)(NknkNWkXNnxn根据基2 FFT的时间抽选法的第一次分解的结果：)()()2()()()(kPWkGNkXkPWkGkX kNkN12, 1 , 0Nkn可以得到：)2()(21)()2()(21)(NkXkXWkPNkXkXkG kN12, 1 , 0Nk 图图 4.8 由由 X(k)、X(k+N/2) 得到得到 G(k)、P(k)n再由N/2点的DFT求得N/4点的DFT，依次类推下去，就可推到求出x(n)的各点，如图4.9所示。将此流图与图4.4比较，相当于整个流向反过来，此外，因子WNk成为WN-k，

26、还添加了因子1/2。但是，图4.9的IFFT算法不能直接利用按照图4.4编写的FFT算法程序，却可以利用图4.6的频率抽选FFT算法的程序，只需求将X(k)作为输入序列，因子WNk变为WN-k，并且将最后所得的输出序列的每个元素都除以N。n方法方法2 2n 将将DFTDFT的正变换式：的正变换式： n与其反变换式：与其反变换式： 10)()(NnnkNWnxkX10)(1)(NknkNWkXNnxn比较，很容易知道，可以利用FFT算法的程序来计算IFFT，只需求将X(k) 作为输入序列，x(n) 那么是输出序列，另外将因子WNk 变为WN-k， 当然，最后还必需将输出序列的每个元素除以N。n

27、方法方法3n 对对DFT的反变换式取共轭，有：的反变换式取共轭，有：n 1, 1 , 0, )(1)(1)(10*10*NnWkXNWkXNnxNknkNNknkNn与DFT的正变换式比较，可知完全可以利用FFT的计算程序，只需求将X*(k) 作为输入序列，并将最后结果取共轭，再除以N就得到x(n)。4.6 4.6 线性调频线性调频z z变换变换CZTCZT算法算法 假设所需求的频率抽样点并不均匀地分布在单位圆上，此时，假设所需求的频率抽样点并不均匀地分布在单位圆上，此时，经常采用经常采用Chirp zChirp z变换变换CZTCZT算法算法 4.6.1 根本原理根本原理用用CZT算法可以计

28、算以下给定点算法可以计算以下给定点zk上的上的X(z)即在这些点处即在这些点处的复频谱：的复频谱： (4.26)式中，式中， ， ，W0 、A0 为正实数。为正实数。1M,0,1,k AWzkk,00jeWW00jeAA n这些z平面上的抽样点所沿的周线是一条螺旋线，参数W0控制周线盘旋的倾斜率。假设W0大于1，那么随着k的添加，周线向内盘旋，趋向原点；假设W0小于1，那么随着k的添加，周线向外盘旋。假设W0=1，螺旋线实践上是一段圆弧，而假设又有A0=1，那么这段圆弧是单位圆的一部分。A0和0分别为第一个抽样点k = 0的模和幅角，其他抽样点沿螺旋周线按角度间隔0分布。 图图 4.10 z

29、平面上一条螺旋周线上的抽样点平面上一条螺旋周线上的抽样点n要计算： n式中N为序列x(n) 的长度。由恒等式： n并且令 ， n就可以得到 1, 0,)()()()(1010Mk WAnxznxzXnxCZTNnnknNnnkk)(21222nknknk22)()(nnWAnxny22)(nWnh)()()(1022nkhnyWzXNnkkn改动变量，将n换为r，k换为n，那么有：n n=0,1,M-1n其中， 可以看作线性时不变系统h(n) 对输入信号y(n) 的呼应。 )()()()()()(22102222ngWnhnyWrnhryWzXnnNrnn)()()(nhnyng 图 4.11

30、 CZT 算法的计算过程 4.6.2 算法要点算法要点此算法主要是要计算此算法主要是要计算y(n)与与h(n) 的线性卷积的线性卷积g(n)。由于序列。由于序列x(n) 长度为长度为N，因此，因此y(n)的长度也为的长度也为N0nN-1。虽然。虽然 可以是无限长的，但是由于可以是无限长的，但是由于z 平面上的抽样点只需平面上的抽样点只需M 个，个，因此因此h(n)中只需有限个点值是所需求的。中只需有限个点值是所需求的。 2/2)(nWnhn由于只需求计算M个X(zn) 之值，那么也只需求M个g(n) 之值，也即对于线性卷积n n只需求求n = 0,1, M-1时的值。如图4.12所示，由于y(

31、r) 的范围是r由0到N-1，因此当求g(n) 的第一个值即n = 0时，也需求h(-r) 中r由0到N-1这N个值。 10)()()()()(Nrrnhrynhnyngn为了计算g(1)到g(M-1), 还应该将h(-r)向右移M-1次，也就是说，h(-r)中由r =-1到r =-(M-1)这M-1个值也是所需求的。这样，我们需求并且也只需求h(-r)中的由r =-(M-1)到r =N-1这L个值，而n L=(N-1)-(M-1)+1 = N+M-1n这也就是说，对于序列h(r)(或h(n)，我们只需求其中r或n由 -(N-1) 到M-1这一范围的L个值。M-1r 0h(r)-(N-1)(b

32、)-(M-1)r0h(-r)N-1(c) 图图 4.12 序列序列 y(n) 和和 h(n) 的长度和所取范围的长度和所取范围n假设要用L点循环卷积来计算线性卷积g(n)，那么序列y(n) 后面应补M-1个0，使其长度为L；而用作循环卷积的有限长序列h(n)为：n (4.32) 1 )(10 )()( 2/)(2/22LnMWLnhMnWnhnhLnnn然后计算L点循环卷积：n n其结果的L个值中，gL(0)到gL(M-1)这M个值正好与所需求的线性卷积g(0) 到g(M-1) 对应相等，这是由于，图4.13(c)中的序列 与图4.12(c中的序列h(-r)在r=-(M-1)到r=N-1区间完

33、全一样。而gL(n)的后N-1个值那么是不需求的。)( )( )()(10nRrnhryngLLrL)( rh n假设计算循环卷积时需求利用FFT来进展快速计算，那么要求L = 2s，s为正整数。此时，假设M+N-1 2s，那么要在y(n)后面补上M-1+K个零，使L=M+N-1+K=2s。当然，(4.32)式所示的序列h(n) 中也要补K个零，这K个零应该插在图4.13(b)中虚线所示的地方，即在序列h(r)中r=M-1之后插入K个零实践上，可以是K个恣意的数，这样所得到的y(n)与h(n)的循环卷积依然只取前面的M个值，即为所要求的线性卷积的M个值。 4.5.4 CZT算法的特点算法的特点

34、 1计算量计算量 按照图按照图4.11所示的计算过程，可以大致统计出用所示的计算过程，可以大致统计出用CZT算法算法计算计算N点长的时域序列点长的时域序列x(n) 在在z平面的一段螺旋状周线上的平面的一段螺旋状周线上的M个点处的复频谱所需求的乘法次数。个点处的复频谱所需求的乘法次数。(1) 计算 需求N 次复乘。(2) y(n)的L点FFT运算需求 次复乘。(3) 假设坚持M和N不变，那么DFT H(k) 可以事先算好存起来，而不需求现场运算操作。(4) 求G(k) = Y(k)H(k) 需求L次复乘。2/2)()(/)()(nnnWAnxnhAnxny)2/(log)2/(2LL(5) 对G

35、(k) 进展L点IFFT运算以得到g(n)，需求 次复乘。(6) 将g(n) 乘上 求得M点输出X(zn)，需求M次复乘。 这些复数乘法中主要是两次L点FFT(IFFT)的运算量，因此，CZT算法的运算量大致与 成正比，这里LN+M-1。)2/(log)2/(2LL2/2nW2log22LLn假设用z变换直接计算复频谱，即：n n那么需求MN次复数乘法。实践上，当M、N较小时，直接用z变换式来计算会比较方便；乘积MN越大，CZT算法的运算量的减少越显著。1M,0,1,k znxzXNnnkk10)()( 2 2与与FFTFFT算法比较算法比较与规范的与规范的FFTFFT算法相比较，算法相比较，

36、CZTCZT算法有以下特点：算法有以下特点：(1) (1) 输入和输出序列长度不需求相等，即不需求输入和输出序列长度不需求相等，即不需求M = NM = N。(2) N(2) N与与M M都不需求是都不需求是2 2的正整数幂。的正整数幂。(3) zk (3) zk 间的间隔可以恣意选择，这样就可得到恣意的分辨间的间隔可以恣意选择，这样就可得到恣意的分辨率；而且当所需间隔不同时，可以分段计算。率；而且当所需间隔不同时，可以分段计算。 (4) zk 点所沿周线不用需是圆弧。(5) 起始点的选择可以是恣意的，因此便于从恣意的频率上开场对输入数据进展频谱分析。(6) 当(4.26)式中A=1，W=e-

37、j2/N ，并且N=M时，X(zk)即为x(n)的DFT，因此利用CZT算法也能快速计算x(n)的DFT，而且不要求N为2的正整数幂。 4.7 4.7 实序列的实序列的FFTFFT的高效算法的高效算法 上面讨论的上面讨论的FFTFFT算法是将算法是将x(n)x(n)作为复数来对待的，假设作为复数来对待的，假设x(n)x(n)是实序列，计算量还可以进一步减少。是实序列，计算量还可以进一步减少。4.7.1 两个长度一样的实序列两个长度一样的实序列 可以将两个长度一样的实序列组合成一个复序列来进展可以将两个长度一样的实序列组合成一个复序列来进展FFT运算，从而一次完成这两个实序列的运算，从而一次完成

38、这两个实序列的FFT，减少了总，减少了总的计算量。的计算量。n 设p(n) 和g(n) 是两个长度均为N的实序列，并设：n n又设 ， ， n那么由DFT的线性有： (4.36)()()(njgnpny)()(kPnpDFT )()(kGngDFT )()(kYnyDFT )()()(kjGkPkYnP(k) 和G(k) 这两个复序列的实部都是周期性的偶序列，而其虚部都是周期性的奇序列。n 对复序列Y(k) 又有 (4.37)n这里下标 r、i 分别表示实部和虚部，Y(k) 与其实部、虚部的长度都为 N。现将 (4.37) 式中各序列作周期延拓，有 (4.38)()()(kjYkYkYir)(

39、)()(kYjkYkYirn由周期性有：n (4.39)n (4.40)()(),()(kNYkY kNYkYrrrr)()(),()(kNYkY kNYkYiiiin如今将序列 与 作如下分解：n (4.41)n (4.42)(kYr)(kYi)()(21)()(21)(kNYkYkNYkYkYrrrrr)()(21)()(21)(kNYkYkNYkYkYiiiiin 由4.39式和4.40式，容易证明在(4.41)和(4.42)这两个式子中，前一项都是偶序列，而后一项都是奇序列。n 将4.41式和4.42式代入4.38式，并将各项进展重新组合，得到：)( )( )()(21)()(21)(

40、)(21)()(21)(kGjkPkNYkYjkNYkYjkNYkYjkNYkYkYrriiiirrn令0kN-1，那么上式为： n这里的P(k) 和G(k) 的实部都是周期性偶序列，而它们的虚部都是周期性奇序列，此情况与4.36式中的复序列P(k)和G(k)的情况一样。因此有：)()()( )( )(kjGkPkjGkPkY 上两式中0kN-1。)()(21)()(21)( )(NiiNrrkNYkYjkNYkYkPkP)()(21)()(21)( )(NrrNiikNYkYjkNYkYkGkG 4.7.2 4.7.2 一个一个2N2N点的实序列点的实序列 将一个将一个2N2N点的实序列点的

41、实序列x(n)x(n)按偶数点和奇数点分组构成两按偶数点和奇数点分组构成两个个N N点实序列：点实序列： 1, 1 , 0 ) 12()()2()(Nnnxngnxnpn那么有：n k=0,1,N-1 (4.48)()()()()()( k 2k 2kGWkPNkXkGWkPkXNNn其中：n n实序列p(n) 和g(n) 的DFT P(k) 和G(k) 可以采用4.7.1节所说的方法作一次N点复序列的FFT而同时得到，然后再按4.48式进展组合便得到了2N点实序列x(n) 的DFT。10)()(NnnkNWnpkP1, 1 , 0)()(10NkWngkGNnnkN4.8 Matlab方法方法 4.8.1 利用利用Matlab计算计算FFT在第在第3章中引见用章中引见用Matlab方法计算信号的方法计算信号的DFT时，提到了函时，提到了函数数fft(x,N) 和和ifft(x.N)。对于这两个函数，假设。对于这两个函数，假设N为为2的正的正整数幂，那么可以得到本章中引见的基整数幂，那么可以得到本章中