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

第4章 FFT 4.1 引言 4.1.1 离散傅里叶变换的矩阵表示及其运算量 DFT在数字信号处理中起着非常重要的作用，这是与DFT存在着高效算法，即快速傅里叶变换(FFT)分不开的。快速运算的关键是减少运算量。,离散傅里叶变换对为：(4.1)(4.2)式中。下面要用矩阵来表示DFT关系。令：,并令TN、表示两个变换方阵，有：,如果用i(0iN-1)表示这两个N阶方阵的行号，用j(0jN-1)表示这两个N阶方阵的列号，那末很容易看出，TN 方阵中i行j列的元素为，而方阵中i行j列的元素为。于是（4.1）式可以写成：而（4.2）式可以写成：,一般情况下，信号序列x(n)及其频谱序列X(k)都是用复数来表示的，WN当然也是复数。因此，计算DFT的一个值X(k)需要进行N次复数乘法（与1相乘也包括在内）和N-1次复数加法；所以，直接计算N点的DFT需要进行N2 次复数乘法和N(N-1)复数加法。,显然，直接计算N点的IDFT所需的复乘和复加的次数也是这么多。当N足够大时，N2 N(N-1),因此，DFT与IDFT的运算次数与N2 成正比，随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，因此有必要对DFT与IDFT的计算方法予以改进。,4.1.2 因子的特性 DFT和IDFT的快速算法的导出主要是根据 因子的特性。1周期性：对离散变量n有同样的周期性。,2对称性：或 3.其它:,4.2 基2时间抽选的FFT算法 4.2.1 算法推导 已经知道：令DFT的长度N=2M，M为正整数。,令：于是有：,其中，是由x(n)的偶数抽样点形成的DFT；而,是由x(n)的奇数抽样点形成的DFT。但是这两个式子并不完全是N/2点的DFT，因为k的范围仍然是由0到N-1，因此，还应该进一步考虑k由N/2到N-1范围的情况。,现在令，故对于后半段有：同理：又知：,图 4.1 N点DFT分解为两个N/2点的DFT(N=8),图 4.2 N/2 点的DFT 分解为两个N/4点的DFT(N=8),综上所述，可以得到：其中G(k)、P(k)分别是x(n)的偶数点和奇数点的N/2点DFT。,这样，我们就将一个N点的DFT分解成了两个N/2点的DFT，由于DFT的运算量与其点数的平方成正比，因此使运算量减少了。但是，还应该将每一个N/2点的DFT再分解为两个N/4点的DFT，如此下去，直到分解为2点的DFT为止，总共需要进行log2N-1=log2(N/2)次分解。,图 4.3 2 点 DFT 信号流图（蝶形图）,对于2点DFT，有：所以2点DFT的运算只需一次加法和一次减法，这样的运算叫做蝶形运算，这样的信号流图叫做蝶形图。,该算法每次分解都是将时域序列按奇偶分为两组，因此要求N等于2的正整数幂，故将这种FFT算法叫做基2时间抽选法。,4.2.2 算法特点 1.倒序重排这种FFT算法的每次分解都是将输入序列按照奇偶分为两组，故要不断地将每组输入数据按奇偶重排，直到最后分解为2点的DFT，输入数据才不再改变顺序。这样做的结果，使得作FFT运算时，输入序列的次序要按其序号的倒序进行重新排列。,现在将图4.4中输入序号以及重排后的序号按二进制写出如下（注：下标“2”表示二进制数）。可以看出，将输入序号的二进制表示（n2n1n0）位置颠倒，得到（n0n1n2），就是相应的倒序的二进制序号。因此，输入序列按倒序重排，实际上就是将序号为（n2n1n0）的元素与序号为（n0n1n2）的元素的位置相互交换。,2.同址计算 从图4.4可以看出，整个算法流图可以分为四段，（0）段为倒序重排，后面3段为3(log28=3)次迭代运算：首先计算2点DFT，然后将2点DFT的结果组合成4点DFT，最后将4点DFT组合为8点DFT。因此，对于N点FFT，只需要一列存储N个复数的存储器。,开始时，输入序列的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个蝶型，因此共有N次复数加法；此外，除了第一次的2点DFT之外，每次迭代还包括了N/2次复数乘法（即乘WN的幂）。因此，一个N点的FFT共有复数乘法的次数为：,复数加法的次数为：因此，FFT算法比直接计算DFT大大减少了运算量，尤其是当N较大时，计算量的减少更为显著。比如，当N=1024时，采用FFT算法时复数乘法的次数，不超过直接计算DFT时复乘次数的千分之五。,4.2.3 关于FFT算法的计算机程序在一般情况下，进行FFT运算的序列至少都有几百点的长度，因此需要编制FFT算法程序以便能够利用计算机来快速进行计算。可以参照图4.4来编制计算N(N必须等于2的正整数幂)点FFT的计算机程序，整个程序可以分为两部分：一部分是倒序重排，另一部分是用三层嵌套的循环来完成M=log2N次迭代。,三层循环的功能是：最里的一层循环完成蝶形运算，中间的一层循环完成因子WNk的变化，而最外的一层循环则是完成M次迭代过程。倒序重排的程序是一段经典程序，它以巧妙的构思、简单的语句用高级编程语言来完成了倒序重排的功能。下面是一段用FORTRAN语言编写的倒序重排程序。,N=2*M（表示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),IF(I.GE.J)GOTO 5（如果IJ，就转移到语句5）（将输入序列中序号互为倒序 的两个数值交换位置）,5 K=NV2 6 IF(K.GE.J)GOTO 7(如果KJ，就转移到语句 7)J=J-kK=K/2 GOTO 6（转移到语句6）7 J=J+K,由于是同址计算，因此程序中只用了一个数组X()，这个数组共有N个元素。X()既表示输入序列，也表示按倒序重排后的这N个数值。程序中的字母都是大写的，这是FORTRAN语句的特点。,I既是循环变量，也代表输入序列的正序序号，J代表倒置后的序号。由于FORTRAN语言的循环变量不能从0开始，故以8点FFT为例，I的范围为18，下面是正序序号I与倒序序号J之间的对应关系：I:1 2 3 4 5 6 7 8 J:1 5 3 7 2 6 4 8,输入序列按倒序重排，实际上就是将X(I)和X(J)互换位置。显然，如果I=J，就不需要交换；而如果IJ，就表示X(I)与X(J)已经交换过了。因此，在循环语句下面的语句“IF(I.GE.J)GOTO 5”就表明了交换的条件：只有当IJ时才执行下面三条语句，使X(I)与X(J)互换位置。,正序序号I也是循环变量，从1开始，每次循环增加1。关键问题是对于每一个I，所对应的J是什么？程序中从语句5开始直到循环结束就是为下一次循环确定所对应的J的。虽然I-1和J-1的二进制表示是互为倒序的关系，但是，程序中并不是由I得到J，而是由上一个J来得到下一个J。思路是：I-1由0到7，用二进制来表示就是从0002开始，每次在最低位加1，逐次增加到1112。,既然J-1的二进制表示是I-1二进制表示的倒置，那末J-1的变化就应该是在最高位逐次加1，而最高位的二进制数1表示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 基2频率抽选的FFT算法 时间抽选法是在时域内将输入序列x(n)逐次分解为偶数点子序列和奇数点子序列，通过求子序列的DFT而实现整个序列的DFT。而频率抽选法则是在频域内将X(k)逐次分解成偶数点子序列和奇数点子序列，然后对这些分解得越来越短的子序列进行DFT运算，从而求得整个的DFT。当然，同样要求N为2的正整数幂。,设，则可以分别表示出k为偶数和奇数时的X(k)。其中，,其中，显然，X(2r)为g(n)的N/2点DFT，X(2r+1)为p(n)WNn 的N/2点DFT。这样，一个N点的DFT分解为两个N/2点的DFT。将分解继续下去，直到分解为2点的DFT为止。当N=8时，基2频率抽选的FFT算法的整个信号流图如图4.6所示。,将图4.6与图4.4比较，可知频率抽选法的计算量与时间抽选法相同，而且都能够同址计算。时间抽选法是输入序列按奇偶分组，故x(n)的顺序要按倒序重排，而输出序列按前后分半，故X(k)的顺序不需要重排；频率抽选法则是输出序列按奇偶分组，故X(k)的顺序要按倒序重排，而输入序列按前后分半，故x(n)不需要重排。,4.4 基 4时间抽选的FFT算法 上面讨论的FFT算法都要求N=2M，M是正整数，因此是基2的FFT算法。若N=4M，则还可以采用基4的FFT算法。与推导基2 FFT算法类似，可以推导基4 FFT算法。将x(n)相间地分为四组，所得的X(k)按前后分为四组。第一次分组的流图如图4.7所示。其中，中间的矩阵为变换方阵T4，,而 Di(i=0,1,2,3)则为N/4 阶的对角矩阵，即：,即变换方阵，即：,图 4.7 基 4 时间抽选的 FFT算法的第一次分解,为了理解基4 FFT算法，可以将基2时间抽选的FFT算法的第一次分解与基4的第一次分解进行比较。,基2的情形第一次分解为两组，基4为四组。基2分解中的G(k)和P(k)分别为x(n)的偶数点和奇数点的N/2点DFT，而基4分解的每一项中的因子,都代表一个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的分解相对应，又比之复杂。,基4的情形也要继续分解下去，直到分解为四点DFT为止，共要经过log4N-1=log4(N/4)次分解。估算基4 FFT的计算量：一个N点的FFT，要经过log4N次变换。复数乘法体现在与对角矩阵D1、D2、D3相乘，总的复乘次数为：；总的复加次数为：。,基2 FFT算法的复乘次数为：复加次数为：,将Mc4与Mc2比较，Ac4和Ac2比较，可知，基4的复乘次数约减少了1/4，但复加次数却增加了。因此，在乘法运算要转化为加法运算来实现的情况下，基4算法的计算量会比基2算法少；但是，近年来，在许多微处理器和DSP中，实现一次乘法运算与实现一次加法运算的时间完全相同，这时采用基4算法就不合算了，而且基4算法还比基2算法复杂。,4.5 快速傅里叶反变换（IFFT）IFFT是IDFT的快速算法。由于DFT的正变换和反变换的表达式相似，因此IDFT也有相似的快速算法。可以用三种不同的方法来导出IFFT算法。方法1 设，则有：，n=0,1,N-1,根据基2 FFT的时间抽选法的第一次分解的结果：,可以得到：,图 4.8 由 X(k)、X(k+N/2)得到 G(k)、P(k),再由N/2点的DFT求得N/4点的DFT，依次类推下去，就可推到求出x(n)的各点，如图4.9所示。将此流图与图4.4比较，相当于整个流向反过来，此外，因子WNk成为WN-k，还增加了因子1/2。但是，图4.9的IFFT算法不能直接利用按照图4.4编写的FFT算法程序，却可以利用图4.6的频率抽选FFT算法的程序，只需要将X(k)作为输入序列，因子WNk变为WN-k，并且将最后所得的输出序列的每个元素都除以N。,方法2 将DFT的正变换式：与其反变换式：,比较，很容易知道，可以利用FFT算法的程序来计算IFFT，只需要将X(k)作为输入序列，x(n)则是输出序列，另外将因子WNk 变为WN-k，当然，最后还必须将输出序列的每个元素除以N。,方法3 对DFT的反变换式取共轭，有：,与DFT的正变换式比较，可知完全可以利用FFT的计算程序，只需要将X*(k)作为输入序列，并将最后结果取共轭，再除以N就得到x(n)。,4.6 线性调频z变换（CZT）算法 如果所需要的频率抽样点并不均匀地分布在单位圆上，此时，常常采用Chirp z变换（CZT）算法,4.6.1 基本原理用CZT算法可以计算下列给定点zk上的X(z)（即在这些点处的复频谱）：(4.26)式中，W0、A0 为正实数。,这些z平面上的抽样点所沿的周线是一条螺旋线，参数W0控制周线盘旋的倾斜率。若W0大于1，则随着k的增加，周线向内盘旋，趋向原点；若W0小于1，则随着k的增加，周线向外盘旋。若W0=1，螺旋线实际上是一段圆弧，而如果又有A0=1，则这段圆弧是单位圆的一部分。A0和0分别为第一个抽样点（k=0）的模和幅角，其余抽样点沿螺旋周线按角度间隔0分布。,图 4.10 z 平面上一条螺旋周线上的抽样点,要计算：式中N为序列x(n)的长度。由恒等式：并且令，就可以得到,改变变量，将n换为r，k换为n，则有：n=0,1,M-1其中，可以看作线性时不变系统h(n)对输入信号y(n)的响应。,图 4.11 CZT 算法的计算过程,4.6.2 算法要点此算法主要是要计算y(n)与h(n)的线性卷积g(n)。由于序列x(n)长度为N，因此y(n)的长度也为N（0nN-1）。虽然 可以是无限长的，但是由于z 平面上的抽样点只有M 个，因此h(n)中只有有限个点值是所需要的。,由于只需要计算M个X(zn)之值，那么也只需要M个g(n)之值，也即对于线性卷积 只需要求n=0,1,M-1时的值。如图4.12所示，由于y(r)的范围是r由0到N-1，因此当求g(n)的第一个值即n=0时，也需要h(-r)中r由0到N-1这N个值。,为了计算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个值，而 L=(N-1)-(M-1)+1=N+M-1这也就是说，对于序列h(r)(或h(n)，我们只需要其中r（或n）由-(N-1)到M-1这一范围的L个值。,图 4.12 序列 y(n)和 h(n)的长度和所取范围,如果要用L点循环卷积来计算线性卷积g(n)，则序列y(n)后面应补M-1个0，使其长度为L；而用作循环卷积的有限长序列h(n)为：(4.32),然后计算L点循环卷积：其结果的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区间完全相同。而gL(n)的后N-1个值则是不需要的。,若计算循环卷积时需要利用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算法的特点 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次复乘。,(5)对G(k)进行L点IFFT运算以得到g(n)，需要 次复乘。(6)将g(n)乘上 求得M点输出X(zn)，需要M次复乘。这些复数乘法中主要是两次L点FFT(IFFT)的运算量，因此，CZT算法的运算量大致与 成正比，这里LN+M-1。,如果用z变换直接计算复频谱，即：则需要MN次复数乘法。实际上，当M、N较小时，直接用z变换式来计算会比较方便；乘积MN越大，CZT算法的运算量的减少越显著。,2与FFT算法比较与标准的FFT算法相比较，CZT算法有以下特点：(1)输入和输出序列长度不需要相等，即不需要M=N。(2)N与M都不需要是2的正整数幂。(3)zk 间的间隔可以任意选择，这样就可得到任意的分辨率；而且当所需间隔不同时，可以分段计算。,(4)zk 点所沿周线不必须是圆弧。(5)起始点的选择可以是任意的，因此便于从任意的频率上开始对输入数据进行频谱分析。(6)当(4.26)式中A=1，W=e-j2/N，并且N=M时，X(zk)即为x(n)的DFT，因此利用CZT算法也能快速计算x(n)的DFT，而且不要求N为2的正整数幂。,4.7 实序列的FFT的高效算法 上面讨论的FFT算法是将x(n)作为复数来对待的，如果x(n)是实序列，计算量还能够进一步减少。,4.7.1 两个长度相同的实序列 可以将两个长度相同的实序列组合成一个复序列来进行FFT运算，从而一次完成这两个实序列的FFT，减少了总的计算量。,设p(n)和g(n)是两个长度均为N的实序列，并设：又设，则由DFT的线性有：(4.36),P(k)和G(k)这两个复序列的实部都是周期性的偶序列，而其虚部都是周期性的奇序列。对复序列Y(k)又有(4.37)这里下标 r、i 分别表示实部和虚部，Y(k)与其实部、虚部的长度都为 N。现将(4.37)式中各序列作周期延拓，有(4.38),由周期性有：(4.39)(4.40),现在将序列 与 作如下分解：(4.41)(4.42),由（4.39）式和（4.40）式，容易证明在(4.41)和(4.42)这两个式子中，前一项都是偶序列，而后一项都是奇序列。将（4.41）式和（4.42）式代入（4.38）式，并将各项进行重新组合，得到：,令0kN-1，则上式为：这里的P(k)和G(k)的实部都是周期性偶序列，而它们的虚部都是周期性奇序列，此情况与（4.36）式中的复序列P(k)和G(k)的情况相同。因此有：,上两式中0kN-1。,4.7.2 一个2N点的实序列 将一个2N点的实序列x(n)按偶数点和奇数点分组形成两个N点实序列：,则有：k=0,1,N-1(4.48),其中：实序列p(n)和g(n)的DFT P(k)和G(k)可以采用4.7.1节所说的方法作一次N点复序列的FFT而同时得到，然后再按（4.48）式进行组合便得到了2N点实序列x(n)的DFT。,4.8 Matlab方法 4.8.1 利用Matlab计算FFT在第3章中介绍用Matlab方法计算信号的DFT时，提到了函数fft(x,N)和ifft(x.N)。对于这两个函数，如果N为2的正整数幂，则可以得到本章中介绍的基2 FFT快速算法；,如果N既不是2的正整数幂，也不是质数，则函数将N分解成质数，得到较慢的混合基 FFT算法；如果 N 为质数，则fft函数采用原来的 DFT 算法。这里我们不再详细介绍，具体用法参见3.6节。,4.8.2 用Matlab实现有限长序列的Chirp z变换在Matlab中实现有限长序列的CZT(Chirp z-transform)算法的函数为：y=czt(x,M,W,A)该函数的用法。（1）y=czt(x,M,W,A),该函数返回信号x的线性调频z变换（CZT）的值y，它是信号x沿着W 和A定义的螺旋线进行的z 变换。其中M是信号x的长度（点数），W是z平面上感兴趣的那部分螺旋线上抽样点之间的比值，A是螺旋线上的复数起始点。,（2）y=czt(x)同样返回信号x的线性调频z变换（CZT）的值y，其中M、W和A都取缺省值，即有：W=exp(2*j*pi/M),A=1可以看出，对于这些缺省值，y返回 x 信号在单位园上等间隔的M个点上的z变换，也就是x的离散傅立叶变换(DFT)。,例4.1 已知序列x(n)=(0.5)n,0n12，试计算序列x(n)在单位圆上的CZT,并与该序列的DFT 进行比较。解：运行结果如图4.14所示，从图中可以看出，此时序列在单位圆上的CZT就等于该序列的DFT。,例4.1 已知序列x(n)=(0.5)n,0n12，试计算序列x(n)在单位圆上的CZT,并与该序列的DFT进行比较。解：运行结果如图4.14所示，从图中可以看出，此时序列在单位圆上的CZT就等于该序列的DFT。,图4.14 序列 在单位圆上的CZT及其DFT,例4.2 假设序列x(n)由4个频率分别为6 Hz、6.3 Hz、9 Hz和8 Hz的正弦序列组合而成，抽样频率为40 Hz，时域抽样200点。（1）用 CZT 计算DFT；（2）直接计算DFT；（3）在510Hz的频段范围求CZT。,(a)CZT,(b)DFT,图4.15 CZT的应用,在图4.15中，（a）图是利用CZT求取DFT，（b）图是直接求DFT，两者结果相同；图中频率为6Hz和6.3Hz的两个正弦信号的频谱不易分辩出。图（c）是在510Hz这个频率范围内求出的CZT，它的起始点不在处而且，由于它的分点比较细，所以图中四个正弦信号的频谱都可以分辨出来。,