数字信号处理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位置颠倒，得到n0n

9、1n2），就是相应的倒序的二进制序号。因此，输入序列按倒序重排，实际上就是将序号为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次迭代运算。每次迭代运算。每次迭代运算包含了次迭代运算包含

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

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

13、示N=2M ，M是输入的正整数） 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，就转移到语句5）n （将输入序列中序号互为倒序n 的两个数值交换位置） TIXIXJXJXT)( )()()( 5 K=NV2 6 IF(K.GE.J) GOTO 7 (如果KJ，就转移到语句 7)J=J-kK=K/2 GOTO 6 （转移到语句6）7 J=J+K n由于是同址计算，因此程序

14、中只用了一个数组X( )，这个数组共有N个元素。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)已经交换过了。因此，在循环语句下面的语

15、句“IF (I.GE.J) GOTO 5就表明了交换的条件：只有当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

16、表示N/2。所以，程序中将J与N/2比较来判断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而实现整而实现整

17、个序列的个序列的DFTDFT。而频率抽选法则是在频域内将。而频率抽选法则是在频域内将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(NnNNnrN

18、NNNnNNnNnrNnNnrNNnrnNWngWWNnxWnxWNnxWnxWnxrX11)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

19、点的DFT。将分解继续下去，直到分解为2点的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算法都要求算法都

20、要求N=2M N=2M ，M M是正整数，因此是是正整数，因此是基基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-

21、11 1 1 111111119464346444243424144WWWWWWWWWT而 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的第一次

22、分解进行比较。n )()()()(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分解

23、既与基2 FFT的分解相对应，又比之复杂。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 快速傅里叶反变换快速傅里叶反变换IFFTIFFT）IFFTIFFT是是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比较，相当于整个流向反过来，此外

26、，因子WNk成为WN-k，还增加了因子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， 当然，最后还必须将输出序

27、列的每个元素除以N。n 方法方法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 基本原理基本

28、原理用用CZT算法可以计算下列给定点算法可以计算下列给定点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分布。

29、图图 4.10 z 平面上一条螺旋周线上的抽样点平面上一条螺旋周线上的抽样点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)()()(nhny

30、ng 图 4.11 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时的值。如图

31、4.12所示，由于y(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

32、(r)-(N-1)(b)-(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

33、)到r=N-1区间完全相同。而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

35、) 对G(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

37、=e-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的正的正整数幂，则可以得到本章中介绍的基整数幂，则可以得到本章中介绍的基2 FFT快速