《数字信号处理-理论与实践》课件第5章

第5章DFT与FFT5.1

DFT的定义及物理意义5.2离散傅里叶级数（DFS）与DFT5.3

DFT的性质与定理5.4频域采样5.5基2时间抽取FFT算法5.6基2频率抽取FFT算法5.7

IDFT的快速实现——IFFT算法5.8

MATLAB实现习题5.1DFT的定义及物理意义

5.1.1DFT的定义

设序列x(n)长度为M，则它的N点离散傅里叶正变换（DFT）和N点离散傅里叶逆变换（IDFT）分别定义为(5.1-1)0≤k≤N－1(5.1-2)0≤n≤N－1式中，，N为DFT变换区间长度，N≥M。通常称式(5.1-1)和式(5.1-2)为离散傅里叶变换对。由DFT的定义可见，DFT使有限长时域离散序列与有限长频域离散序列之间建立起了对应关系，故可利用计算机完成两者间的变换，这是DFT的最大优点之一。5.1.2DFT的物理意义

设X(ejω)=DTFT［x(n)］,X(z)=Z［x(n)］,X(k)=DFT

［x(n)］N,则(5.1-3)(5.1-4)0≤k≤N－1(5.1-5)0≤k≤N－1

式(5.1-3)就是序列的Z变换与序列的傅里叶变换之间的关系，即单位圆上的Z变换就是序列的傅里叶变换。式(5.1-4)和式(5.1-5)指出了离散傅里叶变换(DFT)与离散时间序列的傅里叶变换(DTFT)以及序列的Z变换之间的关系，即序列x(n)的

N点DFT是对x(n)的频谱X(ejω)在［0,2π］上的N点等间隔采样，采样间隔为2π/N，也就是对序列频谱的离散化。这就是DFT的物理意义。既然离散傅里叶变换(DFT)就是对序列频谱的N点等间隔采样，那么当变换区间长度N变化时，所得的变换结果即DFT也是不同的。类似采样定理，当N越大时采样的谱线越密，

X(k)的包络线就越逼近X(ejω)的幅频特性曲线。5.2离散傅里叶级数（DFS）与DFT

5.2.1离散傅里叶级数(DFS)

我们知道，对于周期信号，通常都可以用傅里叶级数来描述。若连续时间周期信号为f(t)=f(t+mT)，

m为任意整数，则根据第1章的周期信号傅里叶级数展开公式，信号f(t)的指数形式傅里叶级数表示为(5.2-1)上式可以看成是信号f(t)被分解成不同频率的各次谐波的叠加，每个谐波都有一个幅值|F(jnΩ)|，表示该谐波分量所占的比重。其中ejΩt为基波，基频Ω=2π/T（T为信号周期）。同样，现设是周期为N的周期序列，即r为任意整数因为周期序列不满足绝对可和的条件，所以它的傅里叶变换和Z变换不存在。但是周期序列同连续时间周期信号一样，可以用离散傅里叶级数表示。对于周期序列，用指数形式的傅里叶级数表示应该为(5.2-2)式中，ω0=2π/N是基频，ejω0n为基频序列，ejkω0n

为k次谐波序列。虽然在表现形式上和连续周期函数是相同的，但是离散周期序列的傅里叶级数的谐波成分只有N个独立成分，这与连续傅里叶级数的无穷多个谐波成分相比是不同的。为了说明这一点，下面来分析一下第k+rN（r为任意整数）次谐波ej(k+rN)ω0n

和第k次谐波ejkω0n之间的关系。因为r为任意整数这说明第k+rN次谐波能够被第k次谐波代表，即复指数序列是k的周期函数。因此，在所有的谐波成分中，只有从k=0到

N－1的N个谐波成分是独立的，用这N个谐波就可完全地表示出。另外，为了计算方便，引入一个系数1/N，这就形成了周期序列的离散傅里叶级数，即(5.2-3)下面来讨论如何根据来求解系数。将式(5.2-3)两端同时乘以e－j(2π/N)rn，并对n在从0到N－1的一个周期内求和，得到利用复指数的正交性，即因此将r换成k，可得(5.2-4)式(5.2-4)就是计算k=0到N－1的N个谐波系数的公式，通常把称为的离散傅里叶级数。从式(5.2-4)中也可以看出，也是周期为N的周期序列，即有习惯上，我们常采用以下符号因此，将上式代入式(5.2-4)和式(5.2-3)，并重写如下(5.2-5)(5.2-6)式(5.2-5)和式(5.2-6)称为离散傅里叶级数变换对，式中的n和k都为离散变量。通常把n和k分别看做时间和频率变量，用DFS［·］表示时域到频域的变换，称为离散傅里叶级数正变换，用IDFS［·］表示由频域到时域的变换，称为离散傅里叶级数反变换。观察离散傅里叶级数的正反变换式，我们有以下结论：

（1）不论是还是，都是周期为N的无限长序列，因此，离散傅里叶级数是两个同周期的无限长序列之间的映射，即时域周期序列的离散傅里叶级数在频域也是一个周期序列。

（2）对于周期序列，我们只要知道其中一个周期的内容，则整个周期序列的信息也都知道了。因而周期序列的信息可以用N个序列值（一个周期）来代表。所以，不论是DFS还是IDFS，虽然是从无限长序列到无限长序列之间的变换，由于其周期特征，其计算中也仅需对序列的一个周期进行求和。从这一点来看，离散傅里叶级数也可以看做是建立起了从有限长序列到无限长序列的变换关系。实际上，任意一个长度为N的有限长序列x(n)都可以看做是一个周期为N的周期序列的一个周期；反过来说，任意一个周期为N的周期序列，都可以看做是它的一个周期所形成的N点序列进行以N为周期的周期延拓的结果。通常，我们将周期序列（假设周期为N）的第一个周期所处的区间［0,N－1］定义为“主值区间”，而处于主值区间的部分称为“主值序列”。

（3）对于周期是N的序列，因为不是绝对可和的，其Z变换不收敛，所以不能用Z变换表示。但若只取

的一个周期则其Z变换存在，即(5.2-7)对比式(5.2-5)与式(5.2-7)，可知(5.2-8)式(5.2-8)表明，是的一个周期形成的有限长序列x(n)的Z变换在z平面的单位圆上以等间隔角进行采样得到的。

同样，x(n)的傅里叶变换为(5.2-9)对比式(5.2-5)与式(5.2-9)，可知(5.2-10)式(5.2-10)表明，是的一个周期形成的有限长序列x(n)的傅里叶变换以为间隔进行等间隔采样得到的。因而，离散傅里叶级数实质上是建立了周期序列与有限长序列之间的映射关系。

例5-1设为周期冲激串,即求的DFS。解当0≤n≤N－1时，=δ(n)，因而利用式(5.2-5)可得到即对于所有的k值，均为1。因此，利用式(5.2-6)，表示成级数形式为例5-2设的周期N=10，在主值区间内，当0≤n≤4时，=1，当5≤n≤9时，=0。假设的主值序列为x(n)，其傅里叶变换为X(ejω)。试求和X(ejω)，并分别画出它们的幅度特性图。

解按照式(5.2-5)，有

按照式(5.2-9)，有图5-1例5-2的DTFT与DFS幅度图X(ejω)的幅度特性图（一个周期）和的幅度特性图（一个周期）如图5-1所示。从图5-1可以看出，相当于在ω=0到ω=2π的范围内，以2π/10的频率间隔在10个等间隔的频率上对x(n)的傅里叶变换进行采样。这里也说明了DFS实际上代表了序列各个频率分量的幅度。

注意，若x(n)不是取主值周期的序列，而是随便取

的一个周期序列，比如n取5~14，则再次计算离散傅里叶级数，得到的结果和取主值周期的序列计算得到的结果是一样的。5.2.2DFT与DFS

上面已经提过，对于有限长序列与周期序列x(n)可看做是的一个周期（主值序列），而可看做是x(n)的以N为周期的周期延拓。所以，不论n的取值如何，我们都有(5.2-11)通过引入求余符号((n))N来表示式(5.2-11)中的(n模N)，同时利用矩形序列的符号RN(n)，就可以将周期序列与其主值序列x(n)之间的关系表示为(5.2-12)同理，频域的周期序列也可看做是对频域的有限长序列X(k)的周期延拓，而有限长序列X(k)可看成是周期序列主值序列，即(5.2-13)我们回头再来看式(5.2-5)与式(5.2-6)的DFS及IDFS的表达式，其中的求和都只限定在主值区间内进行，故完全适用于主值序列。假设时域和频域的主值序列分别为x(n)与X(k)，则有(5.2-14)与(5.2-15)可以看出，式（5.2-15）恰为N点DFT的表达式，这说明DFT与DFS之间有着紧密的关系：有限长序列x(n)的N点离散傅里叶变换X(k)，就是以N为周期将x(n)进行周期延拓后的周期序列x((n))N的离散傅里叶级数的主值序列。式（5.2-14）说明，周期序列的离散傅里叶数，就是的主值序列（假设为N点）x(n)的N点离散傅里叶变换X(k)以N为周期的周期延拓序列。所以，既然是x(n)的傅里叶变换以2π／N为间隔的等间隔采样，那么X(k)作为

的主值序列，就是x(n)的傅里叶变换以2π／N为间隔在［0,2π］区间的等间隔采样。

因此，离散傅里叶变换(DFT)与周期序列的离散傅里叶级数(DFS)之间存在本质的联系。虽然DFT是针对有限长序列的，但根据傅里叶变换的时域、频域的对偶关系，频域的离散化必然对应着时域序列本身的周期性。因此，进行离散傅里叶变换的原始序列必然是周期的，即DFT隐含着序列的周期性，有限长序列都是作为周期序列的一个周期来表示的，都隐含有周期性意义。

5.3DFT的性质与定理

由于DFT与DTFT和Z变换之间以及DFS之间的本质联系，因此它们的许多性质具有相似性。但是，由于在DFT变换中x(n)和X(k)均为有限长序列，因此它们之间的性质还有一些重要差别。DFT的许多性质在数字信号处理中有广泛的应用。本节讨论DFT的一些主要性质。以下讨论中，假设序列x1(n)和x2(n)都是有限长序列，且设

X1(k)=DFT［x1(n)］,X2(k)=DFT［x2(n)］5.3.1DFT的隐含周期性

DFT的正反变换式定义了一对N点离散傅里叶变换，其中n和k的取值范围均定义为0~N－1，如果式（5.1-1）中k的取值域为［－∞,+∞］，就会发现X(k)是以N为周期的，即

X(k)=X(k+mN),m为任意整数(5.3-1)

我们把X(k)的这一性质称为DFT的隐含周期性。证明

DFT的隐含周期性可以从两种不同的角度得出：

（1）如前综述，X(k)是对X(ejω)的采样，由于X(ejω)是以2π为周期的周期函数，即X(k)是对X(ejω)在主值区［0,2π］上的N点等间隔采样。显然，当自变量k超出DFT变换区间时，必然得到［0,2π］之外区间上X(ejω)的采样，且以N为周期重复出现，得到=X((k))N。

（2）利用X(k)与x(n)的周期延拓序列x((n))N的离散傅里叶级数之间的关系，也可以得出DFT的隐含周期性。5.3.2线性性质

设有限长序列x1(n)和x2(n)的长度分别为N1和N2，若

x(n)=ax1(n)+bx2(n),a和b为常数

则有

X(k)=aX1(k)+bX2(k),k=0,1,…,N－1(5.3-2)

式中，N≥max［N1,N2］，X(k)、X1(k)、X2(k)分别是x(n)、x1(n)、x2(n)的N点DFT。注意，这里以N为变换长度，相对较短的序列需通过尾部补零使其长度增加为N。

该性质可以直接用DFT的定义式证明，请读者自己练习。5.3.3循环移位性质

1.有限长序列的循环移位

设序列x(n)的长度为M，对x(n)以N为周期进行周期延拓，得到则定义x(n)的循环移位序列为(5.3-3)式(5.3-3)表示，将序列x(n)以N为周期进行周期延拓得到的周期序列，再移位m个单位之后，得到的周期序列仍旧是一个周期为N的序列，并取的主值序列，就得到x(n)的循环移位序列y(n)。当m<0时，表示向右循环移位。x(n)及其循环移位的过程如图5-2所示，图(a)、(b)、(c)、(d)分别对应序x(n)、、和y(n)的波形，其中M=6,N=8,m=2。

2.循环移位性质

设序列x(n)的长度为M，x(n)的循环移位序列y(n)为

y(n)=x((n+m))NRN(n),N≥M

令X(k)=DFT［(x(n)］，则

Y(k)=DFT［y(n)］=WkmNX(k),0≤k≤N－1

(5.3-4)

图5-2序列的循环移位过程证明令n+m=l，得到上式中的x((l))N和WklN都是以N为周期的，对其在任意一个周期的求和结果相等。将上式的求和区间设在主值区间，就有因此5.3.4DFT的共轭对称性

我们知道，序列傅里叶变换的对称性是关于坐标原点的纵坐标的对称性。DFT也有类似的对称性，但在DFT中涉及的序列x(n)和X(k)均为有限长序列，且取值区间均为主值区间［0,N－1］，所以讨论对称性时，不能再以坐标原点作为对称点，而是以n=N/2点作为对称点。为了区别于无限长共轭对称序列，我们用xep(n)和xop(n)分别表示有限长（或圆周）共轭对称序列和共轭反对称序列，其定义分别为(5.3-5)(5.3-6)如同任何实数都可以分解为偶对称分量和奇对称分量一样，任何有限长序列x(n)都可以用它的共轭对称分量和共轭反对称分量之和来表示，即(5.3-7)其中(5.3-8)(5.3-9)同理可得(5.3-10)(5.3-11)式(5.3-10)和式(5.3-11)中的Xep(k)和Xop(k)分别表示有限长序列X(k)的共轭对称分量和共轭反对称分量。利用上述公式，可得出DFT的共轭对称性质：

（1）DFT［x*(n)］=X*(－k)=X*(N－k)(5.3-12)式中,x*(n)表示x(n)的共轭序列。

证明因为所以

（2）DFT［x*(－n)］=X*(k)(5.3-13)证明因为DFT的周期性，取x(n)的主值序列求和就有（3）如果序列x(n)可以表示为

x(n)=xr(n)+jxi(n),0≤n≤N－1其中则有(5.3-14)(5.3-15)式(5.3-14)和式(5.3-15)表明，如果序列x(n)的DFT为X(k)，则x(n)实部的DFT对应于X(k)的共轭对称分量Xep(k)，x(n)虚部的DFT对应于X(k)的共轭反对称分量Xop(k)。

证明利用式(5.3-8)和式(5.3-9)，可分别得到因此式(5.3-16)表明，序列x(n)的共轭对称分量的DFT为X(k)的实部；式(5.3-17)表明，序列x(n)的共轭反对称分量的DFT为X(k)的虚部乘以j。

上述公式(5.3-12)到公式(5.3-17)是序列DFT共轭对称的基本公式，下面给出x(n)为实序列时的DFT共轭对称公式，这些公式可以利用上述基本公式进行证明。

（5）实信号DFT的共轭对称性。设x(n)为长度为N的实序列，且X(k)=DFT［x(n)］。

①X(k)共轭对称，即

X(k)=X*(N－k),0≤k≤N－1(5.3-18)②如果x(n)实偶对称，即

x(n)=x(N－n)

则X(k)也实偶对称，即

X(k)=X(N－k)(5.3-19)

③如果x(n)实奇对称，即

x(n)=－x(N－n)

则X(k)纯虚奇对称，即

X(k)=－X(N－k)

(5.3-20)

实际中实序列的DFT应用非常广泛。利用上述性质可知，只要x(n)是实序列，则其离散傅里叶变换X(k)必然共轭对称。所以只要计算出X(k)的前一半N/2个值，后一半的X(k)值即可利用对称性求得。这样可以减少近一半的运算量，提高运算

效率。5.3.5循环卷积定理

对应于傅里叶变换中的线性卷积定理，DFT有循环卷积定理。下面首先介绍循环卷积的概念，然后介绍循环卷积定理和计算循环卷积的方法。

1.两个有限长序列的循环卷积

设序列h(n)和x(n)的长度分别为N和M。h(n)和x(n)的L点循环卷积定义为(5.3-21)式中，

L称为循环卷积的长度，且L≥max［M,N］。为了区别于第3章和第4章介绍的线性卷积，这里用来表示循环卷积，即yc(n)=x(n)

h(n)。

循环卷积可以采用图解法、列表法等多种方法来计算，在计算机中常采用矩阵相乘的方法来计算。下面介绍如何用矩阵计算来表示循环卷积的公式。在式(5.3-21)中，令n=0，则在求和区间内由x((n－m))L形成的关于x(n)的循环倒相序列为

{x(0),x(L－1),…,x(1)}

与由x((n))L的主值区间x(n)取值形成的序列{x(0),x(1),…,

x(L－1)}相比，该循环倒相序列相当于第一个序列值x(0)不动，将后面的序列反转180°再放在x(0)后面。在式(5.3-21)中，令n=1，则在求和区间内由x((n－m))L形成的x(n)循环倒相序列为

{x(1),x(0),x(L－1),…,x(2)}

观察可以发现，它相当于将上述n=0时形成的关于x(n)的循环倒相序列向右循环移一位而成。

继续增加n的取值，可以发现，所得到的序列都是在前一次移位的基础上再向右移一位形成的。因此，当m和n均从0变化到L－1时，我们就得到式(5.3-21)中关于

x((n－m))L

取值的矩阵为我们把上面的矩阵称为序列x(n)的L点循环卷积矩阵。其特点是：

（1）矩阵的第1行是x(n)的循环倒相序列，即由x((n))L的主值区间的x(n)取值形成的序列。注意，此时如果x(n)的长度M<L，则在x(n)后面补L－M个零。

（2）矩阵第1行以后的各行均由前一行向右循环移动一位所形成。

（3）矩阵各主对角线上的值均相等。

利用序列x(n)的L点循环卷积矩阵，式(5.3-21)可以写成矩阵形式(5.3-22)例5-3已知x(n)=R4(n)，h(n)={1,2,0,1}，试分别求x(n)和h(n)的4点与8点循环卷积。

解利用式(5.3-22)，计算x(n)和h(n)的4点循环卷积yc(n)=x(n)

h(n)如下同理，计算x(n)和h(n)的8点循环卷积yc(n)=x(n)h(n)如下

2.DFT的时域循环卷积定理

设有限长序列h(n)和x(n)，长度分别为N和M，L=max

［N,M］。h(n)和x(n)的L点DFT分别为

H(k)=DFT［h(n)］

X(k)=DFT［x(n)］

如果Y(k)=H(k)·X(k)

则(5.3-23)或证明假设式(5.3-23)成立，对式(5.3-23)的两边进行DFT，可得利用循环卷积公式(5.3-21)，有上式中，令n1=n－m，则有因为上式中x((n1))LWkn1L是以L为周期的，所以对其在任何一个周期上求和的结果都是相等的。因此，对x((n1))L的主值区间0~L－1进行求和，则有时域循环卷积定理表明，两序列时域的循环卷积，相当于两序列的DFT相乘。因此，利用该定理可以使两序列时域的循环卷积计算得到大大简化，计算过程如图5-3所示。图5-3用时域循环卷积定理计算两序列时域的循环卷积

3.DFT的频域循环卷积定理

利用时域与频域的对称性，可以得到DFT的频域循环卷积定理。

设有限长序列h(n)和x(n)，长度分别为N和M，L=max

［N,M］。h(n)和x(n)的L点DFT分别为

H(k)=DFT［h(n)］

X(k)=DFT［x(n)］

如果

y(n)=h(n)x(n)

则或关于DFT频域循环卷积定理的证明留作习题，请读者证明。

4.线性卷积与循环卷积的关系设x1(n)和x2(n)分别是N点和M点的有限长序列，对x1(n)和x2(n)作线性卷积(5.3-25)则yl(n)是长度为M+N－1点的有限长序列，其长度等于参与卷积的两个序列的长度和减1。再对x1(n)和x2(n)作循环卷积。首先对x1(n)和x2(n)分别补零，使之长度均为L，

L≥max［N,M］，然后进行L点循环卷积

(5.3-26)由于x2((n－m))L是周期延拓序列，因此有代入式(5.3-26)，得到(5.3-27)式(5.3-27)说明，有限长序列x1(n)和x2(n)的线性卷积yl(n)以L为周期的周期延拓序列的主值序列即为这两个有限长序列的循环卷积。yl(n)的长度为M+N－1，要想在周期延拓时不发生混叠，则要求延拓的周期L不小于M+N－1，如果满足了这个条件，延拓序列的主值序列就是yl(n)本身。此时，循环卷积和线性卷积相同。延拓周期L若小于M+N－1，则各延拓周期发生混叠，其循环卷积和线性卷积不相同。

一般在实际情况中常常要求线性卷积，但线性卷积计算复杂，计算量大，因而在满足一定条件下可以利用循环卷积代替线性卷积，简化计算，提高运算效率。5.3.6DFT形式下的帕斯维尔(Parseval)定理

设有限长序列h(n)和x(n)，长度分别为N和M，L=max

［N,M］。h(n)和x(n)的L点DFT分别为

H(k)=DFT［h(n)］

X(k)=DFT［x(n)］

则有

(5.3-28)证明式(5.3-28)中，令h(n)=x(n)，则有上式表明：一个序列在时域计算的能量与在频域计算的能量相等。

5.4频域采样

我们知道，对模拟信号进行时域等间隔采样，则时域采样信号的频谱是原模拟信号频谱的周期延拓。根据时域采样定理，针对带限信号，当采样频率大于等于奈奎斯特采样频率时，可以由时域离散采样信号恢复原来的连续信号，而不丢失任何信息。那么，根据时域和频域的对偶性质，对离散时间序列x(n)的连续频谱在频域等间隔采样，也应该有类似的情况出现。下面就对序列的频域采样加以讨论。

1.频域采样和频域采样定理

在5.2节中已经提到，周期为N的周期序列的离散傅里叶级数的系数与的一个周期x(n)的Z变换在z平面单位圆的N个均分点上的采样值相等。由于可以看做是x(n)的以N为周期的延拓序列，因此也可以这样认为：对长度为N的有限长序列x(n)的Z变换在单位圆上进行N点等间隔采样，则该频域采样序列的离散傅里叶级数的反变换就是周期延拓序列。那么，当x(n)的长度与采样点数不一样的时候又是什么情况呢？因此，需要一般性地讨论任意一个绝对可和的非周期序列x(n)的情况。对于绝对可和的非周期序列x(n)，它的Z变换为由于x(n)绝对可和，因此其傅里叶变换存在且连续，故Z变换的收敛域包括单位圆。如果在z平面的单位圆上，对X(z)以等间隔角2π／N进行采样，就可得到周期为N的频域周期序列(5.4-1)显然，这就是对x(n)的频谱函数X(ejω)在间隔为2π／N的N个均分频率点上的等间隔采样。对周期序列，令其离散傅里叶级数的反变换为，则将式(5.4-1)代入上式，得到由于所以(5.4-2)这说明频域采样序列所对应的时域周期序列是原序列x(n)的以N为周期的周期延拓，也就是说，频域采样会造成时域的周期延拓。根据DFT与DFS之间的关系知道，若分别截取和的主值序列xN(n)与

X(k)，则xN(n)和X(k)构成一对DFT，即由于X(k)是的主值序列，因而其频域的采样就是对X(ejω)在频率区间［0,2π］上的N点等间隔采样。而xN(n)是的主值序列，也就是原序列x(n)的以N为周期的周期延拓序列的主值序列。显然，只有在延拓周期大于等于x(n)的长度时，其周期延拓才不会发生混叠，才能使xN(n)=x(n)。综上所述，可以总结出频域采样定理:

如果原序列x(n)的长度为M，且对其频谱X(ejω)在频率区间［0,2π］上等间隔采样N点，得到的序列为X(k)，则仅当采样点数

N≥M

时，才能由频域采样X(k)不失真地恢复出x(n)=IDFT［X(k)］，否则将产生时域混叠失真，不能利用x(n)=IDFT［X(k)］恢复原序列x(n)。

频域采样定理告诉我们，如果时域序列x(n)是无限长的，则其频域采样必然会造成时域混叠，从而不可能无失真地恢复原序列；只有对有限长序列，才能以适当的采样间隔进行频域采样而不丢失信息。

2.z域内插公式与频域内插公式

我们已经知道，对于时域采样，其时域恢复过程可以用时域内插公式来表示。同样地，对于频域采样，也可由z域内插公式恢复出原序列的Z变换，由频域内插公式恢复出原序列的傅里叶变换。

设序列x(n)的长度为M，在z平面的单位圆上对X(z)进行N点等间隔采样，且满足频域采样定理（即N≥M），则有(5.4-3)(5.4-4)(5.4-5)将式(5.4-5)代入式(5.4-3)，得到由于，因此(5.4-6)令(5.4-7)则(5.4-8)式(5.4-8)就称为用X(k)即N个频率抽样来恢复X(z)的z域内插公式，其中φk(z)称为z域内插函数。

将z=ejω

代入式(5.4-6)，经化简求得其频率响应为(5.4-9)其中(5.4-10)利用时域和频域的对称特性，同样可以分析由式(5.4-9)如何得到序列x(n)的频率响应，请读者将其与时域采样信号的恢复重建相对比，采用类似的方法自行分析其恢复过程。5.5基2时间抽取FFT算法

快速傅里叶变换（FFT）并不是一种新的变换，而是计算离散傅里叶变换（DFT）的一种快速算法。

我们知道，有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列，这在数字信号处理中非常有用。但由于计算量太大，直接利用DFT计算很难实时地处理问题，因此引出了快速傅里叶变换(FFT)。FFT最初是在1965年，由Cooley和Tukey提出的计算离散傅里叶变换（DFT）的快速算法引出的，它将DFT的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT）算法的研究便不断深入，形成了一套高效的运算方法，使DFT的计算大大简化。数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。到目前为止，根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本上可以分为两大类：按时间抽选法（用DIT表示）和按频率抽选法（用DIF表示）。这里我们主要讨论两种基本算法：基2DITFFT（时间抽选FFT）和基2DIFFFT（频率抽选FFT）算法。

根据DFT的定义式在将所有WknN的值计算好的情况下，计算每一个X(k)值需要N次复数乘法和N－1次复数加法。因此计算出全部N点的X(k)即整个DFT运算共需N2次复数乘法和N(N－1)次复数加法。可见，直接计算DFT的计算量与N2成正比，当N很大时，运算量是很可观的。例如，当N=8时，DFT需64次复乘，而当N=1024时，DFT所需复乘为1048576次，即一百多万次复乘运算，这对实时性要求很高的信号处理来说，对计算速度的要求就太高了。因而需要对DFT的计算方法进行改进，以大大减少运算量。观察DFT的计算公式就可以发现，利用WknN因子的下列固有性质，可使DFT运算量减少：

（1）周期性：（2）对称性：利用这两个性质，可以使DFT运算中的有些项合并，以减少乘法次数。例如，当N=4时，求X(2)的值。通过合并，可以使上式的乘法运算次数由4次减少到1次。由于DFT运算中的运算量直接与N2成正比，因此N越小,运算量就越小。因而利用上述特性把长序列的DFT分解为短序列的DFT，可以大大减少运算量。FFT正是基于这种思路而发展起来的。5.5.1基2DITFFT算法原理

为了将大点数的DFT分解为小点数的DFT运算，要求序列的长度N为复合数，最常用的是N=2L的情况（L为正整数）。如果不满足这个条件，可以人为地加上若干零值点，使其满足这个条件。我们把这种N为2的整数幂的FFT称为基2FFT。下面我们首先讨论N/2点DFT的计算，在此基础上说明基

2DITFFT算法的计算原理和过程。

1.N/2点DFT

（1）将序列x(n)按n的奇偶分为两组子序列，这样有n为偶数时：n为奇数时:则序列x(n)的DFT就变为(5.5-1)由于因此，式(5.5-1)可表示为(5.5-2)其中,X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即(5.5-3)(5.5-4)式(5.5-2)表明，一个N点的DFT可以分解为两个N/2点的DFT。但由于x1(r)、x2(r)以及X1(k)、X2(k)都是N/2点的序列，因此利用式(5.5-2)只能得到X(k)的k=0,1,…,(N/２)－1个值，即X(k)的前半部分。要得到全部的X(k)值，就必须应用DFT中WmN因子的性质。（2）X(k)的后半部分的确定。利用WmN因子的周期性，有得到(5.5-5)(5.5-6)又由于所以根据式(5.5-2)，有(5.5-7)这样，只要求出0~(N/2－1)区间的所有X1(k)和X2(k)值，就可根据式(5.5-7)求出N/2到(N－1)区间内的所有X(k)值，也就是说后半部分的X(k)完全可由X1(k)和X2(k)确定。综上所述，X(k)的前后两部分均可由X1(k)和X2(k)表示，而X1(k)和X2(k)即为x1(r)和x2(r)的N/2点DFT。这样，一个N点的DFT完全可由两个N/2点的DFT来计算。为方便起见，以下的X1(k)和X2(k)均指x1(r)和x2(r)的N/2点DFT。（3）蝶形运算。由X1(k)、X2(k)表示X(k)的运算是一种特殊的运算，我们定义为蝶形运算:(5.5-8)式(5.5-8)的蝶形运算可用图5-4所示的蝶形信号流图符号表示。流图的表示方法将在第6章讨论。图5-4按时间抽取的蝶形运算流图符号（4）计算工作量分析。

通过将序列x(n)按奇、偶分组后，可将其N点DFT分解为两个N/2点DFT，然后经过N/2个蝶形运算来完成。从图5-4可以看出，每个蝶形运算需要一次复数乘法及两次复数加法。这样，一个N点DFT分解为两个N/2点DFT后，如果直接计算N/2点DFT，则每一个N/2点DFT只需要(N/2)2次复数乘法和(N/2)(N/2－1)次复数加法，因此两个N/2点DFT共需N2/2次复数乘法和N(N/2－1)次复数加法。再加上利用蝶形运算把两个N/2点DFT合成为N点DFT时的N/2个蝶形运算的计算量，通过这种分解方法计算N点DFT，总共需要复数乘法次数和复数加法次数都近似为N2/2，因此通过这种分解后，运算工作量差不多减少到原来的一半。

现以N=23=8点DFT为例，详细分析这一分解过程。

①n为偶数时的序列为x(0)，x(2),x(4),x(6),分别记作:

x1(0)=x(0),x1(1)=x（2）

x1(2)=x(4),x1(3)=x(6)对这组序列进行N/2=4点的DFT，得②n为奇数时的序列为x(1),x(3),x(5),x(7)，分别记作:

x2(0)=x(1),x2(1)=x(3),x2(2)=x(5),x2(3)=x(7)对这组序列进行N/2=4点的DFT，得因此③对X1（k)和X2(k)进行蝶形运算,这里前半部分项数为X(0)~X(3),后半部分项数为X(4)~X(7)。整个过程如图5-5所示。图5-5按时间抽取，将一个N点DFT分解为两个N/2点DFT（N=8）

2.继续分解

既然如此，由于N=2L，则N/2仍是偶数，可以采用同样的方法进一步将每个N/2点DFT的输入再按奇偶分组，进而将每个N/2点DFT分解为两个N/4点DFT。

首先，将原来n为偶数时的N/2点序列x1(r)分解为(偶中偶）(偶中奇）对这两个N/4点的序列进行DFT，得到从而可得到X1(k)前N/4项为(5.5-9)则X1(k)的后N/4项为(5.5-10)同理，对序列中原来n为奇数的N／2点序列x2(r)再按奇偶分为两个N／4点序列：(奇中偶）(奇中奇）则得到(5.5-11)(5.5-12)其中，X5(k)与X6(k)分别为序列x5(l)和x6(l)的点DFT。下面，我们以将N=8时的DFT继续分解为四个N/4点DFT的计算过程为例进行说明，具体如下。

(1)将原序列x(n)的“偶中偶”部分，即x3(l)=x1(r)=x(n)取出如下

x3(0)=x1(0)=x(0),x3(1)=x1(2)=x(4)

计算所构成序列的N/4点DFT，从而得到X3(0)和X3(1)。

(2)将原序列x(n)的“偶中奇”部分，即x4(l)=x1(r)=x(n)取出如下

x4(0)=x1(1)=x（2），x4(1)=x1(3)=x(6)

计算所构成序列的N/4点DFT，从而得到X4(0)和X4（1)。

(3)将原序列x(n)的“奇中偶”部分，即x5(l)=x2(r)=x(n)取出如下

x5(0)=x2(0)=x(1),x5(1)=x2(2)=x(5)

计算所构成序列的N/4点DFT，从而得到X5(0)和X5(1)。

(4)将原序列x(n)的“奇中奇”部分，即x6(l)=x2(r)=x(n)取出如下

x6(0)=x2(1)=x(3)

x6(1)=x2(3)=x(7)

计算所构成序列的N/4点DFT，从而得到X6(0)和X6(1)。

(5)由X3(0)、X3(1)、X4(0)、X4(1)进行蝶形运算,得到X1(0)~X1(3)。

(6)由X5(0)、X5(1)、X6(0)、X6(1)进行蝶形运算,得到X2(0)~X2(3)。

(7)由X1(0)~X1(3)和X2(0)~X2(3)进行蝶形运算,得到X(0)~X(7)。

整个过程如图5-6所示。图5-6按时间抽取，将一个N点DFT分解为四个N/4点DFT的组合（N=8）这样，经过又一次分解(得到四个N/4点DFT)和两级蝶形运算，其运算量又大约减少一半，即为N点DFT的1/4。

这样的分解可以一直进行到最后，也就是最后剩下的是两点DFT，它只有加减运算，但为了统一运算结构，我们仍然采用系数为W0N的蝶形运算来表示。比如，对于N=8时的DFT，N/4点即为两点DFT，因此亦即

因此，8点DFT的FFT的蝶形运算流图如图5-7所示。图5-78点FFT的时间抽选算法信号流图这种方法按照输入序列在时间上的次序是属于偶数还是奇数，将序列进一步分解为两个更短的子序列，所以称做“按时间抽取法”。5.5.2DITFFT运算量分析

由上述分析可知,当N=2L时共需L级蝶形运算，而且每级都由N/2个蝶形运算组成，每个蝶形运算都有一次复数乘法和两次复数加法。因此，最终计算一个序列的N点FFT，所需的总的复数乘法次数和复数加法次数分别为

复数乘复数加由于计算机的乘法运算比加法运算所需的时间多得多，故通常以乘法作为比较基准。

根据前面的分析我们知道，如果直接用DFT的计算公式计算N点DFT，所需的复数乘法次数是mF=N2。相比之下，在计算大点数的DFT时，FFT的优势体现得更为明显，如N=1024时，直接用定义式计算所需的复数乘法次数为N2=10242=126976，而FFT所需的复数乘法次数为，计算次数不到原来的1/20，因而计算时间大大减少。在N的数值继续增大时，与直接计算DFT相比，FFT的运算效率还将更加突出。5.5.3DITFFT算法的特点

针对任意的N=2L点的DITFFT算法，我们来讨论一下它的运算特点。

1.原位运算

由上述运算流图可知，一共有N个输入/输出行，log2N=L级蝶形运算(基本迭代运算)。设用m(m=1,2,…,L)表示第m级迭代，用k、j表示蝶形输入数据所在的（上/下）行数(k,j=0,1,…,N－1)，则任何一个蝶形运算都可用下面的通用式来

表示(5.5-13)令当m=1时，则有（只表示前两个蝶形）当m=2时，则有（只表示前两个蝶形）当m=3时，则有（只表示前两个蝶形）可见，通过在某级进行蝶形运算的两个节点变量k和j就完全可以确定该蝶形运算的结果，而与其他节点变量无关。这样，蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中，最终实现输入数据、中间运算结果和最后输出均使用同一组存储器的原位运算。因此，虽然DIFFFT算法共有L级且每级均有N/2个蝶形运算，但实际上只需

N个存储单元。可见，这种原位运算结构可以节省存储单元，降低设备成本。

2.倒位序规律

由图5-8可知,输出X(k)按正常顺序排列在存储单元,而输入则是按以下顺序排列的

x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)这种顺序看起来好像“混乱无序”，但实际上也是有规律的，我们将这种顺序称做倒位序,即二进制数倒位。

这种倒位序是由奇偶分组造成的，以N=8为例示于图5-8。图5-8倒位序的树状图

3.倒位序实现

一般实际运算中，总是先将输入序列按自然序存入存储单元，为了得到倒位序的排列，可以通过变址运算来完成。假设输入序列的序号为n，用二进制表示为(n2n1n0)2，则其倒位序的二进制数即为(n0n1n2)2。例如，N=8时的自然序与倒位序的关系如表5-1所示。

有了自然序号n与倒位序号，就可以通过变址处理将按自然顺序存放在存储单元中的数据换成所要求的倒位序。具体的变址功能如图5-9所示。图5-9码位倒置的变址处理当n=时，不必调换；当n≠时，必须将数据x()调入原来存放数据x(n)的存储单元内，而将数据x(n)调入存放x(

)的存储单元内。注意，在具体实现时必须避免把已调换过的数据再次调换（否则又回到原状），这可通过比较n

和的大小来实现：若比n小，则意味着此x(n)已经和

x()调换过，不必再调换了；只有当>n时才进行调换。

尽管变址运算所占运算量的比例很小，但对某些高要求的应用（尤其在实时信号处理中），也可设法用适当的电路结构直接实现之。例如，单片数字信号处理器TMS320C25就有专用于FFT的二进制码变址模式。

4.蝶形运算两节点的距离

以图5-7所示的8点DITFFT为例，其第一级每个蝶形运算的两节点间的距离为1，第二级每个蝶形的两节点间的距离为2，第三级每个蝶形的两节点间的距离为4。由此类推，对于N=2L点的DITFFT，其第m级运算中的每个蝶形的两节点间的距离为2m－1。因此，我们可以将蝶形运算表示为(5.5-14)

5.WrN的确定

在每个蝶形运算中都要乘以因子WrN，我们称其为旋转因子。由于N=2L为已知，因此只需确定r的值。通过观察图5-7不难发现，r的变化是有一定规律的，针对这些规律可以采用多种方法确定r的值。一种简便的方法是：

(1)将蝶形运算两节点中的第一个节点标号值k表示成L位(N=2L)二进制数：k=(nL－1…n1n0)2。

(2)将k=(nL－1…n1n0)2左移(L－m)位(m指m级的运算),右边位置补零,就可得到(r)2的值，即(r)2=(k)22L－m

。例如,N=8=23,有

(1)k=2,m=3时的r值:

k=2=(010)2，左移L－m=3－3=0，则r=(010)2=2

(2)k=3，m=3时的r值：

k=3=(011)2,左移L－m=3－3=0，则r=(011)2=3

(3)k=5，m=2的r值:

k=5=(101)2，左移L－m=3－2=1，则r=(010)2=25.6基2频率抽取FFT算法

同样是基于将长序列分解为短序列的思想，如果把输出序列X(k)（假设也是N点序列）按其顺序的奇偶性分解为越来越短的子序列，就形成了另一种FFT算法，称为按频率抽取（DIF）的FFT算法。

5.6.1DIFFFT算法原理

1.N点DFT的另一种表达形式

仍设序列点数N=2L，L为整数。我们把输入序列x(n)按n的顺序分成前后两部分（注意，这不是频率抽取），可得到N点DFT的另一种表达形式由于，故因此X(k)可进一步表示为(5.6-1)

2.X(k)按奇偶分组当k为偶数时,(－1)k=1；当k为奇数时,(－1)k=－1。因此，按k的奇偶可将X(k)分为两部分。令则(5.6-2)(5.6-3)

3.蝶形运算可以看出,式(5.6-2)为前一半输入与后一半输入之和的N/2点DFT，式(5.6-3)为前一半输入与后一半输入之差再与WnN乘积的N/2点DFT。如果令(5.6-4)则(5.6-5)式(5.6-4)所表示的运算关系可以用图5-10所示的蝶形运算来表示。图5-10按频率抽取蝶形运算流图符号这样，就可以把一个N点DFT按k的奇偶分解为两个N/2点DFT了。比如当N=8时，上述分解过程如图5-11所示。图5-11按频率抽取将一个N点DFT分解为两个N/2点DFT的组合（N=8）

4.继续分解

同时间抽取法的推导过程一样，由于N=2L，N/2仍是一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4点DFT。这两个N/4点DFT的输入也是先将N/2点DFT的输入（即前一级蝶形的输出）上下对半分开后通过蝶形运算之后形成的，这一步分解过程如图5-12所示。图5-12按频率抽取将一个N点DFT分解为四个N/4点DFT的组合（N=8）这样的分解可以一直进行到第L级（N=2L），第L级实际上就是作两点DFT，它只有加减运算，但为了统一运算结构，我们仍然采用系数为W0N的蝶形运算来表示，这N/2个两点DFT的N个输出即为x(n)的N点DFT的结果X(k)。图5-13表示了一个N=8的完整的按频率抽取的FFT运算流图。图5-13按频率抽取的FFT流图（N=8）5.6.2DIFFFT算法特点

1.原位运算

从图5-13可以看出，这种运算和按时间抽取法一样，是很有规律的。其每级(列)都是由N/2个蝶形运算构成的,每一个蝶形结构完成同样的基本迭代运算:(5.6-6)式中,m表示第m级迭代，k、j为数据所在行数。通过进行蝶形运算的两个节点变量k和j完全可以确定蝶形运算的结果，而与其他节点变量无关。因而蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中，实现原位运算。

2.蝶形运算两节点的距离

从图5-13中可以看出，当计算第一级蝶形(m=1)时，参与同一蝶形运算的两节点之间的距离为4；第二级每个蝶形的两节点之间的距离为2，第三级每个蝶形的两节点之间的距离为1。因此对于N=2L的一般情况，可推出其蝶形运算两节点之间的距离为2L－m=N／２m，其中，m是指第m级蝶形，m=1,2,…,L。

3.WrN的计算

由于第m级DIF蝶形运算的两节点间的距离为N／2m，因此该蝶形运算可表示为(5.6-7)现在的问题是在不同级（m）的情况下，如何求解r。一种简便的方法是：

(1)将蝶形运算中第一个节点标号值k表示成L位二进制数：k=(nL－1…n1n0)2。

(2)将k=(nL－1…n1n0)2左移m－1位，右边位置补零，就可得到(r)2的值，即(r)2=(nL－1…n1n0)22m－1。

例如，N=8时，有

(1)m=1，k=2时，k=(010)2，左移1－1=0位，则r=(010)2=2

(2)m=2，k=1时，k=(001)2，左移2－1=1位，则r=(010)2=2

(3)m=2，k=5时，k=(101)2，左移2－1=1位，则r=(010)2=25.6.3DIF法与DIT法的异同

通过对两种算法的原理及特点进行比较，我们可以看出，DIT法与DIF法在以下两个方面是相同的：

(1)两种算法都可以进行原位运算。

(2)两种算法运算量相同,即都有L级（N=2L时）运算，每级运算均需N/2个蝶形，每个蝶形均有一次复乘、两次复加，故总的运算量均为次复乘和Nlog2N次复加。两种算法的不同点在于:

(1)DIT输入为倒位序,输出为自然序，而DIF恰好相反。初看起来，这一点不同是显而易见的。但对于任何流图，只要保持各节点所连的支路以及传输系数不变，则不论节点位置怎么排列，所得流图都是等效的，因而不论是哪种算法，都可以将输入或输出进行重排，使其变为自然序或倒位序，所以这并不是二者的本质区别。

(2)基本蝶形不同。

DIT中的蝶形运算是先作复乘后再作加减法，而DIF的复数乘法则仅出现在减法之后，这一点才是二者本质上的不同。如果将DITFFT的基本蝶形运算写为矩阵形式，则有(5.6-8)同样，将DIFFFT的基本蝶形运算写为矩阵形式，则有(5.6-9)从式(5.6-8)和式(5.6-9)可以看出，两个蝶形运算的传输矩阵互为转置。另外，从图5-4和图5-10也可以看出，如果将DIT的基本蝶形的所有支路方向都反向，并且交换输入和输出，就得到了DIF的基本蝶形；反过来将DIF的基本蝶形的所有支路方向都反向，并且交换输入和输出，就得到了DIT的基本蝶形。因此，DIT与DIF的基本蝶形运算是互为转置的。5.7IDFT的快速实现——IFFT算法

DFT的FFT算法同样可以适用于离散傅里叶反变换的运算，即快速傅里叶反变换（IFFT）。IFFT可以通过以下两种方法来快速实现。一种方法需要稍微改动FFT程序和参数来实现。由于比较两式可知，只要把DFT运算中的每个系数WnkN换成

W－nkN，最后再乘以常数１／N，就可以得到IDFT的快速算法——IFFT。

此外,可以将常数１/N分配到每级运算中,由于

因此每级蝶形运算均乘以1/2即可。

另一种方法则完全不用改变FFT的程序而可以直接实现IFFT。考虑到将上述的IDFT公式两边同时取共轭，可得因此(5.7-1)式(5.7-1)说明，只要先将X(k)取共轭，然后将X*(k)作为输入直接调用FFT程序计算DFT，最后再取一次共轭，并乘以1/N，即可得到x(n)。所以，FFT和IFFT的运算完全可以调用同一个子程序，使用起来十分方便。

5.8MATLAB实现

5.8.1用MATLAB计算序列的DFT与IDFT

1.直接计算DFT与IDFT

根据DFT与IDFT的定义式可以编写MATLAB函数来完成DFT与IDFT的计算。考虑到每计算一个X(k)都需要一个for…end循环来计算N项的和式，为了计算所有的X(k)又需要另一个for…end循环，这就会得到一种两重循环嵌套的实现方式。为了提高计算效率，在MATLAB中可以采用矩阵向量乘法来实现上述计算。令行向量X表示X(k)的各点的值，行向量x表示

x(n)的各点的值，则上面两式可表示为

X=xXN（5.8-1）（5.8-2）式中矩阵WN为0≤n,k≤N－1（5.8-3）

矩阵WN是一个方阵，称为DFT矩阵。因此，可以编写下面的MATLAB函数来实现上述的计算过程。计算DFT的函数为：

function［Xk］=dft(xn,N)

%计算DFT

%［Xk］=dft(xn,N)

%Xk是频域的DFT序列

%xn是时域的有限长序列

%N是DFT的长度

n=［0:1:N－1］;

k=［0:1:N－1］;