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

第4章FFT4.1引言

4.1.1离散傅里叶变换的矩阵表示及其运算量

DFT在数字信号处理中起着非常重要的作用，

这是与DFT存在着高效算法，

即快速傅里叶变换(FFT)分不开的。快速运算的关键是减少运算量。离散傅里叶变换对为： (4.1) (4.2)式中

。

下面要用矩阵来表示DFT关系。令：

并令TN、表示两个变换方阵，有：如果用i(0≤i≤N-1)表示这两个N阶方阵的行号，用j(0≤j≤N-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的计算方法予以改进。2．对称性：

或 3.其它:

4.2基2时间抽选的FFT算法

4.2.1算法推导

已经知道：

令DFT的长度N=2M，M为正整数。令：

于是有：是由x(n)的奇数抽样点形成的DFT。但是这两个式子并不完全是N/2点的DFT，因为k的范围仍然是由0到N-1，因此，还应该进一步考虑k由N/2到N-1范围的情况。现在令，故对于后半段有：

同理：

又知：

图4.1N点DFT分解为两个N/2点的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.32点

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）。这就是所谓的同址计算，这样可以大大节约存储元件。

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赋初值）100DO7I=1,NM1(循环开始，到语句7结束；循环变量I从1开始，到NM1结束，步长为1)

IF(I.GE.J)GOTO5（如果I≥J，就转移到语句5）（将输入序列中序号互为倒序的两个数值交换位置）5K=NV26IF(K.GE.J)GOTO7(如果K≥J，就转移到语句7)J=J-kK=K/2GOTO6（转移到语句6）7J=J+K┇由于是同址计算，因此程序中只用了一个数组X()，这个数组共有N个元素。X()既表示输入序列，也表示按倒序重排后的这N个数值。程序中的字母都是大写的，这是FORTRAN语句的特点。I既是循环变量，也代表输入序列的正序序号，J代表倒置后的序号。由于FORTRAN语言的循环变量不能从0开始，故以8点FFT为例，I的范围为1～8，下面是正序序号I与倒序序号J之间的对应关系：

I:12345678J:15372648输入序列按倒序重排，实际上就是将X(I)和X(J)互换位置。显然，如果I=J，就不需要交换；而如果I>J，就表示X(I)与X(J)已经交换过了。因此，在循环语句下面的语句“IF(I.GE.J)GOTO5”就表明了交换的条件：只有当I<J时才执行下面三条语句，使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算法。与推导基2FFT算法类似，可以推导基4FFT算法。将x(n)相间地分为四组，所得的X(k)按前后分为四组。第一次分组的流图如图4.7所示。其中，中间的矩阵为变换方阵T4，

而

Di(i=0,1,2,3)则为N/4阶的对角矩阵，即：

即变换方阵，即：

图4.7基4时间抽选的

FFT算法的第一次分解为了理解基4FFT算法，可以将基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分解既与基2FFT的分解相对应，又比之复杂。基4的情形也要继续分解下去，直到分解为四点DFT为止，共要经过log4N-1=log4(N/4)次分解。估算基4FFT的计算量：一个N点的FFT，要经过log4N次变换。复数乘法体现在与对角矩阵D1、D2、D3相乘，总的复乘次数为：；

总的复加次数为：。

基2FFT算法的复乘次数为：

复加次数为：

将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根据基2FFT的时间抽选法的第一次分解的结果：可以得到：

图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）算法

如果所需要的频率抽样点并不均匀地分布在单位圆上，此时，常常采用Chirpz变换（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.10z平面上一条螺旋周线上的抽样点要计算：

式中N为序列x(n)的长度。由恒等式：

并且令

，

就可以得到

改变变量，将n换为r，k换为n，则有：

n=0,1,…,M-1其中，可以看作线性时不变系统h(n)对输入信号y(n)的响应。

图4.11CZT算法的计算过程4.6.2算法要点此算法主要是要计算y(n)与h(n)的线性卷积g(n)。由于序列x(n)长度为N，因此y(n)的长度也为N（0≤n≤N-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个值。M-1r0h(r)-(N-1)(b)-(M-1)r0h(-r)N-1(c)

图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.4CZT算法的特点

1．计算量

按照图4.11所示的计算过程，可以大致统计出用CZT算法计算N点长的时域序列x(n)在z平面的一段螺旋状周线上的M个点处的复频谱所需要的乘法次数。(1)计算

需要N次复乘。(2)y(n)的L点FFT运算需要

次复乘。(3)若保持M和N不变，则DFTH(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算法的运算量大致与成正比，这里L≥N+M-1。如果用z变换直接计算复频谱，即：

则需要M×N次复数乘法。实际上，当M、N较小时，直接用z变换式来计算会比较方便；乘积M×N越大，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）式，并将各项进行重新组合，得到：令0≤k≤N-1，则上式为：这里的P’(k)和G’(k)的实部都是周期性偶序列，而它们的虚部都是周期性奇序列，此情况与（4.36）式中的复序列P(k)和G(k)的情况相同。因此有：

上两式中0≤k≤N-1。4.7.2一个2N点的实序列

将一个2N点的实序列x(n)按偶数点和奇数点分组形成两个N点实序列：

则有：

k=0,1,┅,N-1(4.48)其中：

实序列p(n)和g(n)的DFTP(k)和G(k)可以采用4.7.1节所说的方法作一次N点复序列的FFT而同时得到，然后再按（4.48）式进行组合便得到了2N点实序列x(n)的DFT。4.8Matlab方法4.8.1利用Matlab计算FFT在第3章中介绍用Matlab方法计算信号的DFT时，提到了函数fft(x,N)和ifft(x.N)。对于这两个函数，如果N为2的正整数幂，则可以得到本章中介绍的基2FFT快速算法；如果N既不是2的正整数幂，也不是质数，