第3章离散傅里叶变换和快速傅里叶变换
通过第1章的学习，我们已经知道了傅里叶变换。傅里叶变换的本质是建立以“时间”为变量的信号与以“频率”为变量的频谱函数之间的变换关系，换言之，傅里叶变换定义了时域和频域之间的一种变换，或者说映射。这里的“时间”和“频率”变量可以取连续值和离散值，从而形成几种形式的傅里叶变换对。离散傅里叶变换是有限长序列的傅里叶变换，它相当于把信号的傅里叶变换进行等间隔采样。离散傅里叶变换除了在理论上具有重要意义之外，由于存在快速算法，在数字信号处理中的应用也越来越广泛。下面首先回顾已学过的连续时间信号的傅里叶变换，然后详细介绍离散傅里叶变换和快速傅里叶变换。
3.1连续时间信号的傅里叶变换
3.1.1时间连续频率离散的傅里叶变换
周期为T的周期性连续时间函数x(t)，可以展开成傅里叶级数，级数的系数为X(kΩ0)。x(t)和X(kΩ0) 组成变换对，正变换为
X(kΩ0)=1T∫T/2-T/2x(t)e-jkΩ0tdt(3.1.1)
反变换为
x(t)=∑∞k=-∞X(kΩ0)ejkΩ0t(3.1.2)
式中，X(kΩ0)是离散频率的非周期函数； Ω0=2π/T 为离散频谱相邻两谱线的角频率间隔； k为谐波序号。
式(3.1.1)表明，周期性连续时间信号x(t)可分解为由无穷次谐波叠加而成，级数系数的绝对值X(kΩ0)代表谐波成分的大小。
3.1.2时间连续频率连续的傅里叶变换
时间连续的非周期信号x(t)的傅里叶变换的结果是连续的非周期的频谱密度X(jΩ)，变换公式为
X(jΩ)=∫∞-∞x(t)e-jΩtdt(3.1.3)
x(t)=12π∫∞-∞X(jΩ)ejΩtdΩ(3.1.4)
式中，Ω是模拟信号的角频率，
简称模拟角频率，单位是rad/s。由此可见，时域函数的连续性造成了频域函数的非周期性，而时域的非周期性造成了频谱的连续性。
3.2离散傅里叶变换及性质
离散傅里叶变换是有限长度序列的傅里叶表示式，同时它本身也是一个序列，而不是一个连续函数，它相当于把信号的傅里叶变换进行频率间隔采样。由于发明了计算离散傅里叶变换的快速算法，所以，离散傅里叶变换在离散时间信号分析与处理中应用非常广泛。
3.2.1序列的傅里叶变换
对于一般序列x(n)，定义傅里叶变换为
X(ejω)=∑∞n=-∞x(n)e-jωn(3.2.1)
式中，ω为离散信号的圆周频率
，简称圆周频率或数字频率，单位是rad。它和模拟角频率Ω的关系在下面说明。序列的傅里叶变换也称为离散时间傅里叶变换（DTFT）。
由式(3.2.1)得
X(ej(ω+2π))=∑∞n=-∞x(n)e-j(ω+2π)n=e-j2πn∑∞n=-∞x(n)e-jωn
=X(ejω)
所以，X(ejω)是ω的周期函数，周期为2π。
用ejωm乘以式(3.2.1)，并在-π~π内对ω积分，并考虑到虚指函数的正交性，有
∫π-πX(ejω)ejωmdω=∫π-π∑∞n=-∞x(n)e-jωnejωmdω
=∑∞n=-∞x(n)∫π-πejω(m-n)dω=2π∑∞n=-∞x(n)δ(m-n)=2πx(m)
所以
x(n)=12π∫π-πX(ejω)ejωndω(3.2.2)
这是DTFT的反变换公式。式(3.2.1)和式(3.2.2)为序列的傅里叶变换对，即离散时间信号的傅里叶变换对。
我们已经学习了模拟角频率Ω，圆周频率ω，它们之间的关系如下：
ω=ΩTs=2πf/fs(3.2.3)
其中，Ts为采样时间间隔； fs=1/Ts是采样频率。如果在式(3.2.3)中令
f′=f/fs(3.2.4)
则 ω=2πf′，f′称为归一化频率或相对频率。这样便可得到对离散序列作DTFT时频率轴定标的物理解释［1］，如图3.2.1所示。
图3.2.1四种频率
由于存在以下关系：
X(ejω)=X(z)z=ejω=∑∞n=-∞x(n)z-nz=ejω=∑∞n=-∞x(n)e-jωn
所以，序列的傅里叶变换就是单位圆上(z=1)的z变换，故序列的傅里叶变换的一切特性，皆可由z变换得到。表3.2.1是序列傅里叶变换的主要性质，其中，性质13~性质18是傅里叶变换的对称性质，它有助于简化运算和求解，这里做以下说明。
表3.2.1序列傅里叶变换的主要性质
序号序列傅里叶变换
1x(n) X(ejω)
2h(n) H(ejω)
3ax(n)+bh(n) aX(ejω)+bH(ejω)
4x(n-m) e-jωmX(ejω)
5 anx(n) X1aejω
6x(n)ejω0n Xej(ω-ω0)
7 nx(n) jdX(ejω)dω
续表
序号序列傅里叶变换
8x(n)*h(n) X(ejω)H(ejω)
9x(n)h(n) 12π∫π-πX(ejθ)H(ej(ω-θ))dθ
10 x*(n) X*(e-jω)
11x(-n) X(e-jω)
12 x*(-n) X*(ejω)
13 Re［x(n)］Xe(ejω)=X(ejω)+X*(e-jω)2
14jIm［x(n)］Xo(ejω)=X(ejω)-X*(e-jω)2
15 xe(n)=x(n)+x*(-n)2Re［X(ejω)］
16 xo(n)=x(n)-x*(-n)2jIm［X(ejω)］
17x(n)为实序列X(ejω)=X*(e-jω)
Re［X(ejω)］=Re［X(e-jω)］
Im［X(ejω)］=-Im［X(e-jω)］
X(ejω)=X(e-jω)
arg［X(ejω)］=-arg［X(e-jω)］
18 xe(n)=x(n)+x(-n)2(x(n)为实序列)Re［X(ejω)］
19xo(n)=x(n)-x(-n)2(x(n)为实序列)jIm［X(ejω)］
20 ∑∞n=-∞x(n)y*(n)=12π∫π-πX(ejω)Y*(ejω)dω(Parsval公式)
21 ∑∞n=-∞x(n)2=12π∫π-πX(ejω)2dω(Parsval公式)
共轭对称序列定义为满足
xe(n)=x*e(-n)(3.2.5)
的序列xe(n)。当共轭对称序列是实数时，条件变为xe(n)=xe(-n)，即xe(n)为偶对称序列。
共轭反对称序列定义为满足
xo(n)=-x*o(-n)(3.2.6)
的序列xo(n)。当共轭反对称序列是实数时，条件变为xo(n)=-xo(-n)，
即xo(n)为奇对称序列。
任意一个序列x(n)总可以表示成一个共轭对称序列和共轭反对称序列之和，即
x(n)=xe(n)+xo(n)(3.2.7)
其中，
xe(n)=12x(n)+x*(-n)(3.2.8a)
xo(n)=12x(n)-x*(-n)(3.2.8b)
同理，x(n)的傅里叶变换X(ejω)也可分解成共轭对称分量和共轭反对称分量之和，即
X(ejω)=Xe(ejω)+Xo(ejω)(3.2.9)
其中，
Xe(ejω)=12X(ejω)+X*(e-jω)(3.2.10a)
Xo(ejω)=12X(ejω)-X*(e-jω)(3.2.10b)
由以上的定义可以得到列于表3.2.1中的傅里叶变换的一些对称性质。
由第2章离散系统的基本概念可知，线性移不变系统单位采样响应h(n)和系统的频率响应H(ejω)也是一对DTFT变换对，即存在如下关系：
H(ejω)=∑∞n=-∞h(n)e-jωn(3.2.11)
h(n)=12π∫π-πH(ejω)ejωndω (3.2.12)
例 3.2.1设矩形窗
d(n)=1,0≤n≤N-1
0,n为其他值
若h(n)=d(n)，求N=5时系统的频率响应［2］。
解由式(3.2.11)得到
H(ejω)=∑∞n=-∞h(n)e-jωn=∑N-1n=0e-jωn=1-e-jωN1-e-jω=e-jωN/2(ejωN/2-e-jωN/2)e-jω/2(ejω/2-e-jω/2)
=e-jω(N-1)/2sin(ωN/2)sin(ω/2)(3.2.13)
H(ejω)=HR(ejω)+HI(ejω)=H(ejω)expφ(ω)(3.2.14)
可见，矩形窗幅频响应和相频响应分别为
H(ejω)=sin(ωN/2)/sin(ω/2)(3.2.15)
φ(ω)=arctanHI(ejω)HR(ejω)=-ω(N-1)/2(3.2.16)
当N=1时，φ(ω)=0
当N=2时，φ(ω)=-ω/2,0≤ω＜π
-ω/2+π,π≤ω＜2π
当N=3时，φ(ω)=
-ω,0≤ω＜23π
-ω+π,23π≤ω＜43π
-ω+2π,43π≤ω＜2π
当N=4时，……
所以，
φ(ω)=
-ω(N-1)/2,0≤ω＜2πN
-ω(N-1)/2+π,2πN≤ω＜4πN
︙
-ω(N-1)/2+(n-1)π,2(n-1)πN≤ω＜2nπN
︙
-ω(N-1)/2+(N-1)π,2(N-1)πN≤ω＜2π
可见，每隔2πN的整数倍相位翻转（频率响应由负变正或由正变负），所以，每隔2πN的整数倍相位要加上π，如图3.2.2所示。
图3.2.2矩形窗的频率响应
3.2.2离散傅里叶变换
有限长序列在数字信号处理中占有很重要的地位，计算机只能处理有限长序列，3.2.1节讨论的序列的傅里叶变换可以分析有限长序列，但是，无法利用计算机进行数值计算。本节以有限长度序列和周期序列之间的关系为出发点推导离散傅里叶变换(DFT)表示式，首先讨论周期序列的傅里叶级数表示式，它相当于有限长度序列的离散傅里叶变换，然后讨论可以作为周期函数一个周期的有限长序列的离散傅里叶变换。
1. 周期序列的离散傅里叶级数(DFS)
设 x～(n) 是周期为N的周期序列，即［3］
x～(n)=∑∞r=-∞x(n+rN)，r 为任意整数(3.2.17)
周期序列不是绝对可和的，换言之，对于z平面内的任意z值，其z变换都不收敛，即
∑∞n=-∞x～(n)z-n=∑∞n=-∞x～(n)z-n=∞
所以，周期序列不能用z变换表示。
但是，和连续时间周期信号一样，周期序列也可用离散傅里叶级数表示，也就是用周期为N的复指数序列表示。
周期为N的复指数序列的基频序列为
e1(n)=ej2πNn
k次谐波序列为
ek(n)=ej2πNkn
由于存在关系
ej2πN(k+rN)=ej2πNk，r为任意整数
即
ek+rN(n)=ek(n)
所以，离散傅里叶级数的所有谐波中只有N个独立分量，这是和连续傅里叶级数不同之处，后者有无穷多个谐波分量。
将x～(n)展成如下的离散傅里叶级数:
x～(n)=1N∑N-1k=0X～(k)ej2πNkn(3.2.18)
式中，N是常数，选取它是下面的X～(k)表达式成立的需要； X～(k)是k次谐波的系数。下面求解系数X～(k)，要利用下面的性质，即
1N∑N-1n=0ej2πNrn=1N1-ej2πNrN1-ej2πNr=1,r=mN,m为任意整数
0，r为其他值(3.2.19)
将式(3.2.18)两边同乘以e-j2πNrn，并对n=0~N-1的一个周期内求和，则得到
∑N-1n=0x～(n)e-j2πNrn=1N∑N-1n=0∑N-1k=0X～(k)ej2πN(k-r)n
=∑N-1k=0X～(k)1N∑N-1n=0ej2πN(k-r)n= X～(r)
式(3.2.18)的系数X～(k)为
X～(k)=∑N-1n=0x～(n)e-j2πNkn(3.2.20)
由于
X～(k+mN)=∑N-1n=0x～(n)e-j2πN(k+mN)n=∑N-1n=0x～(n)e-j2πNkn= X～(k)
所以，X～(k)也是一个以N为周期的周期序列。时域离散周期序列的离散傅里叶级数的系数仍然是一个周期序列，因而我们把式(3.2.18)和式(3.2.20)一起称为周期序列的离散傅里叶级数(DFS)对。
习惯上采用以下符号：
WN=e-j2πN，WknN=e-j2πNkn
这样，式(3.2.18)、式(3.2.20)又可表示为
X～(k)=DFSx～(n)=∑N-1n=0x～(n)e-j2πNnk=∑N-1n=0x～(n)WknN(3.2.21)
x～(n)=IDFSX～(k)=1N∑N-1k=0X～(k)ej2πNnk=1N∑N-1k=0X～(k)W-knN
(3.2.22)
其中，DFS·表示离散傅里叶级数的正变换； IDFS·表示离散傅里叶级数的反变换。
从上面的表示式可以看出，求和时都只取N点序列值。这一事实说明，一个周期序列虽然是无限长序列，但是，只要研究一个周期(有限长序列)的性质，其他周期的性质也就知道了，因而周期序列和有限长序列有着本质的联系。
周期序列X～(k)可以看成是对x～(n)的一个周期x(n)作z变换，然后将z变换在z平面单位圆上按等间隔角2π/N采样而得到的。令
x(n)=x～(n),0≤n≤N-1
0,
其他
则x(n)的z变换为
X(z)=∑∞n=-∞x(n)z-n=∑N-1n=0x(n)z-n(3.2.23)
将式(3.2.23)与式(3.2.21)比较得到
X～(k)= X(z)|z=ej2πNk
所以，X～(k)是在z平面单位圆上N个等间隔角点上对X(z)的采样值。
连同连续时间信号的傅里叶变换和傅里叶级数，到目前为止，我们已经学习了四种形式的傅里叶变换，如图3.2.3所示。总之，若信号在时域是周期的，那么其频谱一定是离散的，反之亦然； 同样，若信号在时域是非周期的，其频谱一定是连续的，反之也成立。第四种周期序列的离散傅里叶级数在时域和频域都是离散的，且都是周期的，我们可以利用它引出有限长序列的离散傅里叶变换。
图3.2.3四种形式的傅里叶变换
2. 离散傅里叶变换(DFT)的定义
由于周期序列只有有限个序列值有意义，因而它的离散傅里叶级数表达式也适用于有限长序列，这就是有限长序列的离散傅里叶变换。如果把长度为N的有限长序列x(n)看成周期为N的周期序列的一个周期，就可以利用离散傅里叶级数计算有限长序列。
对于一个周期序列x～(n)，定义它的第一个周期的序列值为此周期序列的主值序列(principal value sequence)，用x(n)表示，设周期为N，则有
x(n)=x～(n),0≤n≤N-1
0,n为其他值(3.2.24)
显然，x(n)是一个有限长序列，周期序列x～(n)可以看作将x(n)以N为周期进行周期延拓(periodic delay)的结果，如式(3.2.17)所示。该式可简写为
x～(n)=x((n))N及x(n)= x～(n)d(n)(3.2.25)
其中，d(n)是长度为N的矩形序列，即
d(n)= 1,0≤n≤N-1
0,n为其他值 (3.2.26)
((n))N 是余数运算，也可以记为n mod N，表示以 N 为模对n求余数。
例3.2.2x～(n)是周期为N=8的序列，求n=19和n=-5两数对N的余数。
解由于
n=19=3+2×8
故((19))8=3，又由于
n=-5=3+(-1)×8
故((-5))8=3，即
x～(19)=x(3)，x～(-5)=x(3)
■
同理，频域周期序列X～(k)也可看成是对有限长序列X(k)的周期延拓，而有限长序列X(k)看成周期序列X～(k)的主值序列，即
X～(k)=X((k))N，X(k)= X～(k)d(k)(3.2.27)
由DFS和IDFS定义可知，由于求和运算分别只限定在n=0~N-1和k=0~N-1的主值区间进行，故完全适用于主值序列x(n)与X(k)。由式(3.2.21)和式(3.2.22)可得有限长序列的离散傅里叶变换定义。
正变换：
X(k)=DFTx(n)=∑N-1n=0x(n)WnkN，0≤k≤N-1(3.2.28)
反变换：
x(n)=IDFTX(k)=1N∑N-1k=0X(k)W-nkN，0≤n≤N-1(3.2.29)
或分别表示成
X(k)=DFS［x～(n)］d(k)，0≤k≤N-1(3.2.30)
x(n)=IDFS［X～(k)］d(n)，0≤n≤N-1(3.2.31)
x(n)和X(k)是一个有限长序列的离散傅里叶变换对。由于x(n)和X(k)都是长度为N的序列，都有N个独立值，因而已知其中一个序列，就能唯一地确定另一序列。有限长序列的傅里叶变换是作为周期序列的一个周期表示的，也含有周期性的意义。
离散傅里叶变换和前面曾介绍的序列的傅里叶变换都是处理有限长序列的重要工具，它们之间有什么关系呢？由于序列x(n)在单位圆上的z变换等于序列的傅里叶变换X(ejω)，所以，X(k)是序列傅里叶变换X(ejω)的等间隔采样值，采样间隔为ω=2π/N，即
X(k)=X(ejω)|ω=2πNk
=X(ej2πNk)(3.2.32)
该式表明，序列x(n)的离散傅里叶变换结果X(k)是连续频谱X(ejω)的等间隔采样。那么，这种频域的采样需要满足什么条件才能保证由X(k)不失真地恢复X(ejω)呢？设X(k)的逆傅里叶变换xa(n)，则
xa(n)=IDFTX(k)=IDFS［X～(k)］d(n)
= 1N∑N-1k=0∑N-1m=0x(m)WkmNW-nkNd(n)= x～(n+rN)d(n)(3.2.33)
该式表明，xa(n)等于x(n)以N为周期进行延拓后再取主值序列，其中，N为绕单位圆一周采样的点数。设x(n)长度为M，若 N<M，则x(n)在进行周期延拓时产生混叠，从而不能使上式成立，即不能恢复X(ejω)。由此引出频率采样定理，即当N≥M时，才能由频率采样值X(k)不失真地恢复 X(ejω)。X(ejω)可以由内插式 (3.2.34)表示：
X(ejω)=∑N-1k=0X(k)ω-2πNk(3.2.34)
其中，
(ω)=1Nsin(ωN/2)sin(ω/2)e-jω(N-1)/2
称为内插函数［4］。证明从略。
3. 离散傅里叶变换推导图解
为了进一步理解离散傅里叶变换的实质，下面给出DFT推导过程的图解［1］，如图3.2.4所示。
设x(t) 是长度为T的连续时间信号，其傅里叶变换为X(jΩ)，理论上它是无限带宽的。如果x(t)为无限长，则可以用长度为T的矩形窗截短。
图3.2.4DFT推导图解
(1) 通过一个周期冲激串去乘待采样的连续时间信号x(t)，以实现时域采样，如图3.2.4(a)和图3.2.4(b)所示。
时域采样的周期冲激串p(t)为
p(t)=δT(t)=∑+∞n=-∞δ(t-nTs)
（3.2.35）
式中，Ts为采样周期。周期冲激串p(t)的傅里叶变换为
P(jΩ)=2πTs∑+∞k=-∞δ(Ω-kΩs)（3.2.36）
其中，Ωs=2πTs为采样频率。此处，对周期冲激串做傅里叶变换，而不是展成傅里叶级数，目的是为了下面求解采样后信号的傅里叶变换。
采样后的信号为
xs(t)=x(t)p(t)=x(t)∑+∞n=-∞δ(t-nTs)=∑+∞n=-∞x(nTs)δ(t-nTs)
（3.2.37）
根据时域上两个信号的乘积对应于频域上两个信号频谱的卷积，再乘以12π，故有
Xs(jΩ)=12π［X(jΩ)*P(jΩ)
］ （3.2.38）
式中，Xs(jΩ)为xs(t)的傅里叶变换； P(jΩ)为p(t)的傅里叶变换。将式（3.2.36）代入式（3.2.38），于是有
Xs(jΩ)=1Ts
X(jΩ)*∑+∞k=-∞δ(Ω-kΩs)（3.2.39）
因为信号与一个单位冲激函数的卷积就是该信号的移位，即
X(jΩ)*δ(Ω-Ωs)=X(j(Ω-Ωs))（3.2.40）
于是
Xs(jΩ)=1Ts∑+∞k=-∞X（j(Ω-kΩs)）(3.2.41)
这就是说，Xs(jΩ)是频率为Ωs的周期函数，它由一组移位的X(jΩ)的叠加所组成。但是，其幅度乘以1Ts，如图3.2.4(c)所示。设X(jΩ)最高截止频率为Ωc，只要Ωs≥2Ωc，采样就不会发生混叠。注意，图3.2.4(c)中左边依然是连续时间信号，右边是连续的频谱。
下面要将连续时间信号离散化，即实现
xd(n)=x(nTs)（3.2.42）
这里我们要强调一下用周期冲激串采样后的信号xs(t)与连续时间信号离散化后的信号xd(n)的重要区别［4］。
(1) xs(t)是一个连续信号，它除了在Ts的整数倍有冲激外，其他值全为零； 而xd(n)只在整数变量n上取值(实际上是引入了时间上的归一化)，它不再明确包含采样周期Ts的信息。
(2) xs(t)对原连续信号x(t)的采样是用冲激强度(面积)来表示的，而xd(n)是有限值。
那么，离散时间序列xd(n)的傅里叶变换和采样后信号xs(t)的傅里叶变换之间有什么关系呢？由于
xs(t)=∑+∞n=-∞x(nTs)δ(t-nTs)（3.2.43）
对式(3.2.43)做傅里叶变换，根据δ(t-nTs)的傅里叶变换是
e-jΩnTs,所以,得到采样后信号的频谱为
Xs(jΩ)=∑+∞n=-∞x(nTs)e-jΩnTs（3.2.44）
现在再对xd(n)做离散时间序列傅里叶变换，即
X(ejω)=∑+∞n=-∞xd(n)
e-jωn=∑+∞n=-∞x(nTs)
e-jωn（3.2.45）
比较式（3.2.44）和式（3.2.45）可见，X(ejω)和Xs(
jΩ)有如式(3.2.46)所示的关系：
X(ejω)=Xs(jΩ)Ω=
ωTs（3.2.46）
根据式(3.2.41)，有
X(ejω)=1Ts∑+∞k=-∞X（j(ω-2πk)/Ts）（3.2.47）
可以看出，X(ejω)就是Xs(jΩ)的重复，唯频率坐标有一个尺度变换，即X(ejω)变成了以2π为周期的函数。因此，xd(n)和x(t)之间的频谱关系是通过先把x(t)的频谱X(jΩ)按式（3.2.41）进行周期重复，然后，按式（3.2.46）的线性频率尺度变换联系起来的，如图3.2.4(d)所示。频谱周期性重复是冲激串采样转换过程中的第一步结果； 而按式（3.2.46）的线性频率尺度变换，可以不严格地看成是由冲激串xs(t)转换到离散时间序列xd(n)时，所引入时间归一化的结果。根据傅里叶变换的时域尺度变换性质，时间轴上有一个1Ts的变换，一定在频率轴上引入一个Ts倍的变化。因此，ω=ΩTs的关系就与从xs(t)到xd(n)的转换过程中，时间轴上有一个1Ts的尺度变换，在概念上是完全一致的［5］。
(2) 再进行频域采样。与时域采样相似，要在频率上乘一个冲激串，如图3.2.4（e）所示。
Q(ω)=∑+∞k=-∞δω-2πNk(3.2.48)
对应于Q(ω)的时域信号为
q(n)=N2π∑+∞m=-∞δ(n-mN)(3.2.49)
频域采样结果为
Xq(ejω)=X(ejω)Q(ω)=X(ejω)∑+∞k=-∞δω-2πNk=∑+∞k=-∞X(ej2πNk)δω-2πNk(3.2.50)
根据频域乘积对应于时域卷积的性质，频域采样后所对应的时域信号为
x～q(n)=xd(n)*q(n)(3.2.51)
因为信号与一个单位采样序列的卷积就是该信号的移位，所以，
x～q(n)是周期为N的周期函数，它由一组移位的xd(n)的叠加所组成，但是，在幅度上乘以N2π，如图3.2.4(f)所示。
数学表达由式(3.2.51a)给出
x～q(n)=xd(n)N2π∑+∞m=－∞δ(n－mN)=N2π∑+∞m=－∞xd(n－mN)(3.2.51a)
可见，如果要时域上不发生混叠，那么，N不能小于原始的实际数据长度M，此即为频率采样定理。
根据频率采样定理，N要大于原始的实际数据长度M。而实际上根据前面的讨论，xd(n)是无限长的（因为采样冲激串p(t)是无限长的），所以，这里有一个矛盾。这个矛盾只能采用将xd(n)截断来解决。设截断的xd(n)的长度为M,那么要求N≥M。
注意，到目前为止，频率还是连续的函数，因为采样后的频率函数为一系列冲激串的叠加。采样后的频率函数为
Xq(ejω)=∑+∞k=-∞X(ej2πNk)δ(ω-2πNk)（3.2.52）
记xd(n)的周期重复(周期为N)信号为x～(n)，根据以上分析可知
x～q(n)=N2πx～(n)（3.2.51b）
x～q(n)对应的离散傅里叶变换(DTFT)为Xq(ejω)，则x～(n)对应的DTFT为
X～(ejω)=2πNXq(ejω)=∑+∞k=－∞2πNX(ej2πNk)δω－2πNk（3.2.53a）
式(3.2.53a)中的冲激串是连续信号，令
X～(k)=X(ej2πNk)（3.2.53b）
就得到了离散频谱。
上面我们得到了周期的离散频谱式(3.2.53b)，其周期为N。可以看出，它应为某一离散周期时间序列的傅里叶级数(DFS)。那么，它对应的时域周期序列是什么呢？根据周期序列DFS和离散时间DTFT的关系，如果一个周期序列的DFS为X(ej2πNk)，则它对应的DTFT为∑+∞k=－∞2πNX(ej2πNk)δω－2πNk。这与式（3.2.53a）完全吻合。
再根据信号的时域与频域有一一对应的关系，可以得到，X～(k)所对应的时域信号为x～(n)，即xd(n)的周期重复，如图3.2.4(g)所示。
以上所得到的是一对离散傅里叶级数(DFS)，各取一个周期就得到离散傅里叶变换(DFT)。
3.2.3离散傅里叶变换的性质
本节讨论DFT的一些性质，它们本质上和周期序列的概念有关，可以由有限长序列及DFT隐含的周期性得到。设有限长序列x1(n)与x2(n)，且
X1(k)=DFT［x1(n)］，X2(k)=DFT［x2(n)］
1. 线性
x1(n)与x2(n)线性组合的离散傅里叶变换等于它们各自离散傅里叶变换的线性组合，即
X3(k)=DFT［ax1(n)±bx2(n)］=aDFT［x1(n)］±bDFT［x2(n)］
=aX1(k)±bX2(k)(3.2.54)
其中，a，b为任意常数。若x1(n)和x2(n)长度均为N，则所得时间序列的长度也为N； 若x1(n)和x2(n)长度分别为N1,N2，则所得时间序列的长度取两者中的最大者，即
N=max(N1,N2)
例如，若N1>N2，取N=N1，这时x2(n)需在尾部补上N1-N2个零值点，从而使其长度与x1(n)的长度相等，再作N=N1点的DFT。
2. 选频性［6］
设复序列x(n)是对复指数函数x(t)=ejrωt采样得到的，即
x(n)=ejrnωT(3.2.55)
其中，r是整数； T为采样周期，ωT=2π/N。对x(n)作傅里叶变换得
X(k)=∑N-1n=0ej2πn(r-k)/N，0≤k≤N-1(3.2.56)
即
X(k)=∑N-1n=0ej2πN(r-k)n=1-ej2πN(r-k)N1-ej2πN(r-k)=N，r=k
0， r为其他值
这说明复指数函数的采样序列的离散傅里叶变换具有正交性。也就是说，当输入频率为 rω 的正弦波时，傅里叶变换后的离散频谱中只有一条谱线取值为N，其余的都为零。如果输入信号是若干频率不同的正弦波的线性组合，经过离散傅里叶变换后，将在不同的谱线位置有对应的输出，因此，离散傅里叶变换算法实质上对频率具有选择性，它相当于频谱分析仪，对信号处理很有用处。
3. 循环移位
设有限长序列x(n)位于0≤n≤N-1区间内，将其左移m位，得到 x(n+m)，这是序列的线性移位，对这两个序列求DFT，前者的求和范围为0~N-1，后者则为-m~-m+N-1，当m不同时，DFT的求和范围也要改变，这给位移序列DFT的求解带来麻烦。为了解决这个问题，以方便DFT的运算，要重新定义序列的移位。首先将x(n)周期延拓成周期序列x～(n)，然后，将x～(n)向左移动m位，再取x～(n+m)的主值区间(0≤n≤N-1)上的序列值。
有限长序列x(n)向左移动m位的循环移位定义为
xm(n)= x～((n+m))Nd(n)(3.2.57)
如图3.2.5所示。
图3.2.5循环移位
从图3.2.5中可以看出，有限长序列循环移位始终局限于n=0~N-1主值区间内。当某些样本从另一端移出该区间时，需将这些样本从此区间的另一端循环移回来，如果我们想象将序列x(n)按反时针方向排列在一个N等分圆周上，向左移动m位的循环移位就是将该序列在圆周上顺时针旋转m位，如图3.2.5(b)所示［3］。
序列循环移位后的DFT为
Xm(k)=DFT［xm(n)］=W-kmNX(k)(3.2.58)
证明
Xm(k)=DFT［xm(n)］=DFS［x((n+m))N］d(k)
=W-kmNDFS［x((n))N］d(k)=W-kmNX～(k)d(k)=W-kmNX(k)
这表明，有限长序列的循环移位，在离散频域中只引入一个和频率成正比的线性相移 W-kmN=ej2πNkm，对频谱的幅度没有任何影响。
证毕。
同理可得频域的移位特性如下：
IDFT［X(k+l)］=WlnNx(n)(3.2.59)
时间函数乘以指数项Wln N，则DFT相当于左移l位，这就是调制信号的频谱平移原理，也称调制定理(modulated theorem)，它说明时域序列的调制等效于频域的循环移位。
4. 对称性
下面分别讨论三种序列的对称性［6］。
1) 实序列的对称性
设x(n)是长度为N的有限长序列，x～(n)是x(n)以N为周期进行周期延拓的结果，周期序列x～(n)可表示为偶对称序列x～e(n)和奇对称序列x～o(n)之和，即
x～(n)= x～e(n)+ x～o(n)(3.2.60)
利用偶函数和奇函数的对称性，有
x～e(n)=12［x～(n)+ x～(-n)］(3.2.61a)
x～o(n)=12［x～(n)- x～(-n)］(3.2.61b)
由于时域周期序列的频谱也具有周期性，它的实部是偶函数，虚部是奇函数，因此，X～(k)具有以下的对称性：
Re［X～(k)］=Re［X～(N-k)］(3.2.62a)
Im［X～(k)］=-Im［X～(N-k)］(3.2.62b)
|X～(k)| = |X～(N-k)|(3.2.62c)
argX～(k)=-argX～(N-k)(3.2.62d)
对于实序列，Re［X～(k)］就是偶对称序列x～e(n)的DFS，Im［X～(k)］就是奇对称序列x～o(n) 的DFS。可以推断，x(n)的DFT X(k)也具有类似性质。因此，利用式（3.2.62）的对称性关系，很容易从一个序列的DFT得到两个有关序列的DFT。
2) 复序列的对称性
周期复序列x～(n)= x～r(n)+jx～i(n)，它的DFT是
X～(k)=∑N-1n=0［x～r(n)+jx～i(n)］e-j2πNnk= X～r(k)+jX～i(k)(3.2.63)
其中，X～r(k)=DFT［x～r(n)］，X～i(k)=DFT［x～i(n)］都是复数。
式(3.2.63)两边的实部和虚部分别相等，则得到
Re ［X～(k)］=Re ［X～r(k)］-Im ［X～i(k)］(3.2.64a)
Im［X～(k)］=Im［X～r(k)］+Re［X～i(k)］(3.2.64b)
因为实部Re［X～(k)］具有对称性，即
Re ［X～(N-k)］=Re ［X～r(N-k)］-Im ［X～i(N-k)］
=Re［X～r(k)］+Im［X～i(k)］(3.2.65a)
虚部具有反对称性，即
Im［X～(N-k)］=Im［X～r(N-k)］+Re［X～i(N-k)］
=-Im［X～r(k)］+Re ［X～i(k)］(3.2.65b)
与式(3.2.64)联立求解得
Re［X～r(k)］=Re ［X～(k)］+Re ［X～(N-k)］2(3.2.66a)
Im［X～r(k)］=Im ［X～(k)］-Im［X～(N-k)］2(3.2.66b)
Re［X～i(k)］=Im ［X～(k)］+Im［X～(N-k)］2(3.2.66c)
Im［X～i(k)］=Re［X～(N-k)］-Re［X～(k)］2(3.2.66d)
因此，一个复序列的DFT可以同时变换成两个序列的DFT，利用两个频域周期序列的实部的偶对称性和虚部的反对称性，只要计算N/2个样本点的值即可，这样运算次数可减少一半。
3) 复序列的共轭对称性
x(n)的共轭复序列x*(n)=xr(n)-jxi(n)，它的离散傅里叶变换为
DFT［x*(n)］=∑N-1n=0x*(n)WknN=∑N-1n=0x(n)W-knN*
=∑N-1n=0x(n)W(N-k)nN*
=X*((N-k))N，0≤k≤N-1(3.2.67)
考虑到主值区间的定义，共轭复序列x*(n)的DFT可以表示成
DFT［x*(n)］=X*(N-k)(3.2.68)
复序列x(n)或共轭复序列x*(n)的实部序列xr(n)的DFT为
DFT［xr(n)］=12DFT［x(n)+x*(n)］
=12［X(k)+X*(N-k)］=Xe(k)(3.2.69)
虚部序列jxi(n)的DFT为
DFT［jxi(n)］=12DFT［x(n)-x*(n)］
=12［X(k)-X*(N-k)］=Xo(k)(3.2.70)
由式(3.2.69)和式(3.2.70)还可以得到
X*e(N-k)=12［X(N-k)+X*(N-N+k)］*
=12［X*(N-k)+X(k)］
=Xe(k)(3.2.71)
-X*o(N-k)=Xo(k)(3.2.72)
这表明Xe(k)具有共轭对称特性，Xo(k)具有共轭反对称特性。如果把Xe(k)分布在N等分的圆周上，以k=0为原点，则左半圆周上的序列与右半圆周上的序列是共轭对称的，也就是模相等、幅角相反，即
|Xe(k)| = |Xe(N-k)|
arg［Xe(k)］=-arg ［Xe(N-k)］
对Xo(k)而言，以k=0 为原点，则左半圆周上的序列与右半圆周上的序列是共轭反对称的，也就是实部相反、虚部相等，即
Re ［Xo(k)］=-Re ［Xo(N-k)］
Im ［Xo(k)］=Im［Xo(N-k)］
利用共轭对称性，可以用一次DFT运算来计算两个实数序列的DFT，从而达到减少计算量的目的。
5. 循环卷积
前面曾介绍过线性卷积，采用翻转、移位、相乘及求和的计算过程，对于序列x(n)和h(n)，其线性卷积的表达式为
y(n)=x(n)*h(n)=∑∞m=-∞x(m)h(n-m)
两个周期为N的周期序列所进行的卷积，称为周期卷积(periodic convolution)，卷积的结果仍是周期为N的序列。对于两个周期为N的序列x～(n) 和h～(n)的周期卷积y～(n)，有
y～(n)= x～(n)h～(n)=∑N-1m=0x～(m)h～(n-m)=∑N-1m=0x((m))Nh((n-m))N
如果仅将周期卷积的结果截取主值序列，即
y(n)= y～(n)d(n)=∑N-1m=0x((m))Nh((n-m))Nd(n)
而x～(n)和h～(n)的主值序列为x(n)和h(n)，则y(n)就称为x(n)和h(n)的圆周卷积(circular convolution)，表示为
y(n)=x(n)h(n)=∑N-1m=0x((m))Nh((n-m))Nd(n)(3.2.73)
该式所表示的卷积过程可以这样理解： 把序列x(n)分布在N等分的圆周上，而序列h(n)经翻转后分布在另一个N等分的同心圆的圆周上，一个圆周相对于另一个圆周旋转移位，在不同的位置上两序列的对应点依次相乘求和，就得到全部卷积序列，因此圆周卷积又称为循环卷积。循环卷积与周期卷积的过程是一样的，前者仅取卷积结果的主值序列。
由以上的分析可以得知，循环卷积运算有两种方法： 一种方法是，首先计算周期卷积，然后取其主值区间 (0≤n≤N-1)内的值； 另一种方法是，先把x(n)的序列值逆时针方向分布在一个圆周(内圆)上，h(n) 按顺时针方向均匀分布在另一个同心圆(外圆)上，如图3.2.6(a)所示，然后，求两个圆上相应序列的乘积，并把N项乘积累加起来作为 n=0 时的循环卷积值 y(0)。若求 n=1 时的值 y(1)，可将外圆h(n)固定，把内圆上的序列x(n)顺时针旋转一个单位(或将x(n)固定，把外圆上的序列h(n)作逆时针旋转)，然后，把对应项的乘积累加起来，即为所求y(1) 值，如图3.2.6(b)所示。这样依次将内圆序列循环移位一周，便可求得所有y(n)值［2］。
图3.2.6序列的循环卷积
接下来讨论时域和频域循环卷积。令x(n)，h(n)，y(n)都是N点序列，其DFT分别是X(k)，H(k)，Y(k)。若y(n)=x(n)h(n)，则
Y(k)=X(k)H(k)(3.2.74)
证明由式(3.2.73)可得
1N∑N-1n=0Y(k)ej2πNkn=y(n)=∑N-1m=0x(m)h(n-m)
=∑N-1m=01N∑N-1k=0X(k)ej2πNkm1N∑N-1l=0H(l)ej2πNl(n-m)
=1N∑N-1k=0∑N-1l=0X(k)H(l)ej2πNln1N∑N-1i=0ej2πN(k-l)i
=1N∑N-1k=0X(k)H(k)ej2πNkn
式(3.2.74)得证，即两个时域序列的循环卷积y(n)的DFT等于它们的DFT的乘积。
证毕。
同理可以得到频域循环卷积。
若y(n)=x(n)h(n)，则
Y(k)=1NX(k)H(k)(3.2.75)
即两时域序列乘积的DFT等于它们的DFT的循环卷积乘以因子1/N。
例3.2.3设两个有限长序列相等，即
x1(n)=x2(n)=1，0≤n≤N-1
0，其他
求此两序列的循环卷积。
解序列的DFT为
X1(k)=X2(k)=∑N-1n=0WknN=N，k=0
0，其他
两序列的循环卷积y(n)的DFT为
Y(k)=X1(k)X2(k)=N2，k=0
0，其他
则卷积序列
y(n)=N，0≤n≤N-1
■
这个例子说明，循环卷积与线性卷积不同之处在于卷积后的序列长度不同，循环卷积的序列长度是N而不是2N-1。如果对序列x1(n),x2(n)各补N个零值点成为2N长序列，则2N点的循环卷积相当于两序列的线性卷积。
对于有限长序列，存在线性卷积和循环卷积两种形式的卷积。由于循环卷积与DFT相对应，因此，在以后的讨论中可以知道，它可以采用快速傅里叶变换算法(FFT)进行运算，且在运算速度上有很大的优越性。然而，我们在实际应用中遇到的均为线性卷积，例如，信号通过线性系统，系统的输出信号y(n)是输入信号x(n) 与系统单位抽样响应的线性卷积，即
y(n)=x(n)*h(n)
若x(n)，h(n)均为有限长序列，能否用循环卷积实现线性卷积是我们颇为关心的问题。设x(n)为N点有限长序列，h(n)为M点有限长序列，它们的线性卷积
y(n)=x(n)*h(n)=∑N-1m=0x(m)h(n-m)=∑M-1m=0h(m)x(n-m),
0≤n≤N+M-2
仍然是一个有限长序列，长度为N+M-1，即
y(n)=y(n),0≤n≤N+M-2
0, n为其他值
如果直接计算，则需 N×M 次乘法运算，(N-1)×(M-1)次加法运算，当N和M值较大时，运算量是比较大的。
现在讨论用循环卷积实现线性卷积的条件。设循环卷积的长度为L，要用循环卷积实现线性卷积，那么循环卷积的长度L必须大于或等于线性卷积的长度N+M-1，即
L≥N+M-1(3.2.76)
如图3.2.7所示，否则，循环卷积周期延拓时会产生混叠。因此，x(n)，h(n)必须扩展为L点序列，取 L=N+M-1，具体步骤如下。
图3.2.7有限长序列循环卷积和线性卷积(N=3，M=5)
(1) 将x(n)，h(n)分别补零增长到L点，即
x(n)=x(n),0≤n≤N-1
0,N≤n≤L-1(3.2.77)
h(n)=h(n),0≤n≤M-1
0,M≤n≤L-1(3.2.78)
(2) 计算x(n)，h(n)两序列L点的循环卷积即等于线性卷积：
y(n)=x(n)h(n)=x(n)*h(n)
或者用DFT求y(n)，则
y(n)=IDFT［X(k)H(k)］
其中，X(k)，H(k)分别为x(n)，h(n)的L点离散傅里叶变换。用循环卷积实现线性卷积过程如图3.2.8所示。
图3.2.8用循环卷积实现线性卷积
3.2.4离散傅里叶变换在应用中的问题
离散傅里叶变换在信号处理中得到广泛应用，利用离散傅里叶变换可以进行谱分析和时域卷积(实现滤波)。但是，由于离散傅里叶变换固有的局限性，当利用离散傅里叶变换作谱分析时可能会出现下面几个问题。
1. 混叠现象
通常待分析的信号x(t)是连续信号，为了能应用离散傅里叶变换需要对连续时间信号进行采样。若x(t)的频率范围为 0≤f≤fm，则当fm≥12fs时，采样信号的频谱中周期延拓分量互相重叠，这就是混叠现象，它会造成x(n)的频谱和原始信号的频谱不一致。
解决混叠问题的唯一方法是保证采样频率足够高，以防止频谱混叠，这意味着通常需要知道原信号的频率范围，以确定采样频率。若已知信号的最高频率fm，为防止混叠，选定采样频率fs≥2fm。但是，很多情况下可能无法估计信号频率，为确保无混叠现象，可在采样前利用模拟低通滤波器将原信号的上限频率fm限制为采样频率fs的一半，这种滤波器被称为抗混叠滤波器。
对于DFT的离散频谱来说，相邻谱线的频率间隔为Δf，通常称为频率分辨率，由它可以确定模拟信号x(t)的周期，也就是时间长度T为
T=1Δf=1fs/N=Nfs(3.2.79)
由频率分辨率还可以确定DFT所需的采样点数N=fs/Δf。我们希望Δf越小越好，但是，Δf越小，N越大，计算量和存储量也随之增大。一般取N为2的整数次幂，以便用FFT计算，若已给定N，可用补零方法使N为2的整数次幂（在实际应用中，应多采些数据，要尽量避免补零）。
2. 频谱泄漏
信号如果在频域上是带限的，则时域上信号为无限长，但是，离散傅里叶变换却是对有限长序列定义的，因此，为了作离散傅里叶变换，在时域上需要进行截断，使得采样后信号x(n)在区间 ［0,N-1］上，这相当于将x(n)和矩形窗函数d(n)相乘，即
x1(n)=x(n)d(n)(3.2.80)
式中，x1(n)表示截断后的序列。时域上两个序列的乘积等于频域上两个序列的傅里叶变换的卷积。已知d(n)的傅里叶变换
D(ejω)=DTFT［d(n)］=e-jN-12ωsin(Nω/2)sin(ω/2)(3.2.81)
显然，矩形窗函数的频谱为sinc函数。如果窗谱是 δ 函数，那么，时域窗宽应为无穷宽，实际上等于没有乘窗函数，则卷积结果仍为X(ejω)。现在窗谱是sinc函数，有一定宽度，X(ejω)和D(ejω)的卷积是x1(n)的DTFT X1(ejω)，从而产生了频谱变宽的泄漏现象。
泄漏现象是由截断造成的，改善泄漏可以增加采样点数N或采用其他形式的截断函数(通常称为窗函数)，这个问题将在功率谱估计中将详细讨论。
泄漏也会引起混叠。由于泄漏使信号的频谱展宽，如果它的高频成分超过了折叠频率(fs/2)就造成了混叠，这种可能性在矩形窗截断时尤为明显，因为矩形窗的频谱旁瓣收敛得比较慢。
3. 栅栏效应
用DFT计算信号频谱，结果是离散的，即只能给出频谱的采样值，而得不到连续的频谱函数，这就像隔一个“栅栏”观看景象一样，只能在一系列离散点上看到真实的景象，而其他点处却看不到，故称这种效应为“栅栏效应”。减小栅栏效应的一个方法就是要使频域采样更密，即增加频域采样点数N，在不改变时域数据的情况下，必然要在时域序列数据的末尾补一些零值点，使DFT计算周期内的点数增加，但不改变原有记录数据。频域采样为 2πNk，N增加，必然使采样点间隔更小，谱线更密，原来看不到的谱分量就有可能看到了。补零对原频谱 X(k)起到插值的作用，使谱的外观更加平滑。
4. 频率分辨率
分辨率是信号处理中的基本概念，频率分辨率是指所用的算法能将信号中两个靠得很近的谱峰保持分开的能力［4］，即为
Δf=fsN=1NTs=1T(3.2.82)
式中，T=NTs是模拟信号的长度。如果Δf不够小，可以通过增加信号的长度使其变小，该分辨率取决于数据窗的长度和形状。窗的长度指的是实际的信号长度，长度越长，分辨率越高。由DFT的定义可知，对一个N点序列进行DFT分析时，相邻谱线的频率间隔为
Δf=fs/N(3.2.83)
它也是频率分辨率的一种定义，如果Δf不够小，可以通过补零的方法改变谱线的间距。
但是，补零不能提高真实的频率分辨能力，由于没有增加有效数据长度，原数据的新的信息没有增加，因此，仍不能将信号中靠得很近的谱峰分开。
例3.2.4设x(t)的最高频率fm不超过3Hz，由采样定理可知，用fs=10Hz，即Ts=0.1s对其采样，不应发生混叠问题。设T=25.6s，即采样所得的x(n)的点数为256，那么，对x(n)作DFT时，能得到最大频率分辨率［1］，即有
Δf=fsN=10256=0.0390625Hz
如果信号x(t)由三个正弦信号组成，其频率分别是f1=2Hz，f2=2.02Hz，f3=2.07Hz，即
x(t)=sin(2πf1t)+sin(2πf2t)+sin(2πf3t )
用DFT求其频谱，其幅频特性如图3.2.9(a)所示。
图3.2.9幅频特性
显然，由于f2-f1=0.02<Δf，所以，不能分辨出由 f2 产生的正弦分量。由于f3-f1>Δf，所以，能分辨由f3产生的正弦分量。
如果增加点数N，即增加数据的长度N，如令N=1024，这时
T=1024×0.1s=102.4s
其幅频特性如图3.2.9(b)所示。
3.3快速傅里叶变换
3.3.1FFT的基本思想
FFT不是一种新的变换，而是DFT的快速算法。由于DFT的计算量很大，在应用上受到很大限制，FFT的出现使DFT的运算大大简化，从而使DFT在实际中得到广泛应用。
现在分析直接计算DFT的运算量。对N点有限长序列x(n)的DFT变换为
X(k)=∑N-1n=0x(n)WnkN，k=0,1,…,N-1(3.3.1)
假设x(n)是复序列，则 X(k)也是复数，X(k)共有N个点，所以，整个DFT运算需要N2次复数乘法和N(N-1)次复数加法。
快速傅里叶变换能减少运算量的根本原因在于，它不断地把长序列的离散傅里叶变换变为短序列的离散傅里叶变换，再利用系数WnkN的对称性和周期性，即
W(nk+N/2)N=-WnkN(3.3.2)
WknN=Wn(N+k)N=Wk(n+N)N(3.3.3)
将DFT运算中的有些项加以合并，达到减少运算工作量的效果。FFT的算法很多，且有专用的芯片和许多免费程序，在实际应用中，读者可以查阅有关数字信号处理的教材及有关文献。本章只介绍时间抽取(decimationintime，DIT)基2FFT算法和频率抽取(decimationinfrequency,DIF)基2FFT算法。
3.3.2时间抽取基2FFT算法
对式(3.3.1)，令N=2M，M为正整数，则可将x(n)按奇偶分成两组，即令n=2r及n=2r+1，r=0,1,…,N/2-1，于是
X(k)=∑N/2-1r=0x(2r)W2rkN+∑N/2-1r=0x(2r+1)W(2r+1)kN
=∑N/2-1r=0x(2r)(W2N)rk+WkN∑N/2-1r=0x(2r+1)(W2N)rk(3.3.4)
式中，W2N=e-j2πN2=e-j2πN/2=WN/2。于是，有
X(k)=∑N/2-1r=0x(2r)WrkN/2+WkN∑N/2-1r=0x(2r+1)WrkN/2=A(k)+WkNB(k)(3.3.5)
式中，A(k)=∑N/2-1r=0x(2r)WrkN/2，B(k)=∑N/2-1r=0x(2r+1)WrkN/2，k=0,1,…，N/2-1。
这样一个N点的DFT已被分解成两个N/2点的DFT。如果利用A(k)，B(k)表达全部的X(k)，必须要利用系数WnkN的周期性，即式(3.3.2)及Ak+N2=A(k)，Bk+N2=B(k)，从而得到
Xk+N2=Ak+N2+Wk+N2NBk+N2=A(k)-WkNB(k)，
k=0,1,…,N/2-1(3.3.6)
当N=8时，A(k),B(k)及X(k)的关系如图3.3.1所示。
图3.3.1N=8 时 A(k)，B(k)及 X(k)之间关系
当A(k)，B(k)仍是高复合数(N/2)的DFT时，可按上述方法继续加以分解，N点DFT可分成M级。令r=4l及r=4l+2，l=0,1，…，N/4-1，则A(k)表示为
A(k)=∑N/4-1l=0x(4l)W2lkN/2+∑N/4-1l=0x(4l+2)W(2l+1)kN/2
=∑N/4-1l=0x(4l)WlkN/4+WkN/2∑N/4-1l=0x(4l+2)WlkN/4
=C(k)+WkN/2D(k)(3.3.7)
Ak+N4=Ck+N4+Wk+N4N/2Dk+N4=C(k)-WkN/2D(k)，
k=0,1,…,N/4-1(3.3.8)
同理，B(k)表示为
B(k)=∑N/4-1l=0x(4l+1)W2lkN/2+∑N/4-1l=0x(4l+3)W(2l+1)kN/2
=∑N/4-1l=0x(4l+1)WlkN/4+∑N/4-1l=0x(4l+3)WlkN/4
=E(k)+WkN/2F(k)(3.3.9)
Bk+N4=E(k)-WkN/2F(k)k=0,1,…,N/4-1(3.3.10)
若N=8，这时C(k)，D(k)，E(k)，F(k)都是2点的DFT，无须再分，即
C(0)=x(0)+x(4)，E(0)=x(1)+x(5)；
C(1)=x(0)-x(4)，E(1)=x(1)-x(5)
D(0)=x(2)+x(6)，F(0)=x(3)+x(7)；
D(1)=x(2)-x(6)，F(1)=x(3)-x(7)
以上算法是将时间下标n按奇、偶不断进行分组，故称时间抽取算法。上述过程如图3.3.2所示，其基本运算单元如图3.3.3所示，由于运算单元呈蝴蝶形，又称蝶形运算单元(butterfly computation unit)图，图中p,q为第m级蝶形运算单元上、下节点序号，且q-p=2m-1。一个蝶形单元可以将运算量减少至一次复数乘法和两次复数加法，即输入端先与WrN相乘，再与另一输入端分别作加减。第m-1级运算(m=1,2,…，M)中，序号为p，q两点只参与这一个蝶形单元的运算，其输出在第m级，且这一蝶形单元也不再涉及别的点，这一特点称为“同址运算”。由图3.3.2可见，按FFT同址运算的特点，FFT输出的X(k)按自然顺序排列在存储单元中，即按X(0),X(1),…，X(7)的顺序排列，而输入x(n)却不是按正常顺序存储的，而是按 x(0),x(4),x(2),…，x(7)次序排序，服从所谓的“码位倒置”的规律，也叫倒位序。分级运算、同址运算和倒位序也是时间抽取基2FFT算法，简称DIT基2FFT算法的主要特点。
图3.3.28点FFT时间抽取流图
图3.3.3第m级蝶形单元
从上面的分析中可以看出，DIT基2FFT共有M级运算，每一级都含有N/2个蝶形单元，每一个蝶形单元又需要一次复数乘、两次复数加运算，那么，完成 M=log2N级共需要 N2log2N次复数乘法和Nlog2N次复数加法。DIT算法所需的运算量与Nlog2N 成正比，而直接运算DFT的运算量与N2成正比，显然，DIT基2FFT算法大大减少了运算量。
3.3.3频率抽取基2FFT算法
3.3.2节讨论的DIT基2FFT算法将输入序列x(n)按时间下标n的奇偶分解为短序列，还有一种FFT算法是将代表频域输出序列的X(k)按频率下标k的奇偶分解成短序列，称为按频率抽取基2FFT算法。下面给出算法的简单推导，仍讨论长度为N=2M的序列x(n)。首先将x(n)按序号分成前后两部分，得
X(k)=∑N-1n=0x(n)WnkN
=∑N/2-1n=0x(n)WnkN+∑N-1n=N/2x(n)WnkN
=∑N/2-1n=0x(n)WnkN+∑N/2-1n=0xn+N2WnkNWNk/2N
=∑N/2-1n=0［x(n)+WNk/2Nx(n+N/2)］WnkN(3.3.11)
式中，WNk/2N=(-1)k。令 k=2r,k=2r+1，r=0,1,…,N/2-1，则
X(2r)=∑N/2-1n=0x(n)+xn+N2WnrN/2(3.3.12)
X(2r+1)=∑N/2-1n=0x(n)-xn+N2WnrN/2WnN(3.3.13)
令g(n)=x(n)+xn+N2
，h(n)=x(n)-xn+N2WnN，则
X(2r)=∑N/2-1n=0g(n)WnrN/2(3.3.14)
X(2r+1)=∑N/2-1n=0h(n)WnrN/2(3.3.15)
由于g(n)和h(n)是两个N/2点序列，所以，式(3.3.14)和式(3.3.15)表示的是N/2点DFT运算。N点DFT按频率k的奇偶分解为两个N/2点的DFT。频率抽取法和时间抽取法一样，由于N=2M，N/2仍是一个偶数，所以，可以把每个N/2点的DFT的输出再进一步分解为奇数组与偶数组，这样把一个N/2点DFT分解为两个N/4点的DFT。这两个N/4点DFT输入也是将N/2点DFT输入的前一半和后一半分开，再通过蝶形运算而形成，类似的分解可以一直进行下去，直到第M次(M=log2N )分解。第M次分解实际上作2点DFT，2点DFT运算包含一次乘法和二次加法运算。N=8时的频率抽取FFT如图3.3.4所示，图中每一蝶形运算如图3.3.5所示，其中，p,q为第m级蝶形运算单元的上、下节点序号，且q-p=N/2m。
图3.3.48点FFT频率抽取结构图
图3.3.5第m级蝶形单元
由以上的讨论可知，频率抽取FFT算法的运算量与时间抽取FFT算法的运算量相同。频率抽取FFT算法也具有同址运算的优点，不过其输入x(n)是正序排列，而输出 X(k)是倒位序。
前面介绍了两种基本的快速傅里叶变换算法，下面对它作一些讲解和点评。
FFT是实现离散傅里叶变换的快速算法，是一种非常重要的工具。仔细分析FFT算法提出的过程，我们可以发现，这是一个发现问题、分析问题和解决问题的过程，是一个不断创新的过程。傅里叶变换可以把时域的信号变换到频域，也可以把频域的信号变换到时域，以便从不同的角度更好地观察和分析信号。但是，必须知道被分析信号的解析表达式。由于在实际应用中，一般不知道被分析信号的解析表达式，所以，只能对实际信号进行采样，得到离散的信号； 同时，分析信号的工具是数字计算机，且只能处理有限长的离散信号，为此，人们提出了离散傅里叶变换(DFT)。DFT架起了时域与频域之间的桥梁，是信号分析的有力工具。但是，DFT涉及复数的乘法和加法，其运算量比实数的大； 当离散信号的点数多时，运算量将更大，不能满足在实际应用对信号处理实时性的要求。对DFT运算量大的问题进行分析发现，运算量与点数的平方成正比。为此，将长度为N的离散信号分成两个N/2长度的离散信号，再分成四个N/4长度的离散信号……以减少计算量； 再利用系数的固有特性，对运算过程中的有些项进行合并，从而实现了DFT的快速运算，即FFT。可以从时域对离散信号进行分解，这就是按照时间抽取(DIT)的FFT算法； 也可以从频域对长的离散信号进行分解，这就是按频率抽取(DIF)的FFT算法。可以将离散信号的长度定为2的N次方，这就是基2FFT，也可以将离散信号的长度定为4的N次方，这就是基4FFT。对FFT的讲解和点评如表3.3.1所示［13］。
表3.3.1FFT的点评要点
讲 解 过 程进 行 点 评
DFT重要、有用； 但是，运算量大发现问题
与N2有关分析问题
将N分解，再利用系数的固有特性进行合并解决问题
输入乱序、输出顺序,进行倒序寻找规律
从DIT到DIF； 从基2FFT到基4FFT不断突破
3.3.4快速傅里叶逆变换
上面讨论的FFT算法同样可以用于IDFT运算，简称IFFT，即快速傅里叶逆变换。将IDFT的定义公式
x(n)=1N∑N-1k=0X(k)W-knN
与DFT公式
X(k)=∑N-1n=0x(n)WknN
比较可以看出，只要把DFT运算中的每一系数WknN改为W-knN，并在最后再乘以常数1/N，那么前面所讨论的FFT算法就可用来计算IDFT。当把频率抽取FFT算法用于计算IDFT时，由于输出变量变成x(n)，相当于按x(n)的下标的奇偶来分组，因而改称为按时间抽取IFFT算法。同样，当把时间抽取FFT算法应用于IDFT时，输入变为X(k)，是按X(k)的奇偶分组的，故改称为按频率抽取IFFT算法。例如，将频率抽取FFT算法用于计算IDFT，作如下修改： 把WnkN换成W-nkN； 每级运算中都乘以因子1/2(将常量1/N分解成M个1/2连乘，可以防止溢出，有利于减小量化误差)； 输入序列改为自然顺序的X(k)，输出序列改为倒位序的x(n)，就可以得到图3.3.6所示的IFFT结构图。
图3.3.6按时间抽取的IFFT结构(N=8)
另外，还有一种完全不改变FFT计算程序就可以计算IFFT的方法。对IDFT公式取共轭运算，有
x*(n)=1N∑N-1k=0X(k)W-nkN*=1N∑N-1k=0X*(k)WnkN(3.3.16)
式(3.3.16)两边再取共轭，有
x(n)=1N∑N-1k=0X*(k)WnkN*=1NDFTX*(k)*(3.3.17)
式(3.3.17)说明，只要先将输入序列X(k)取共轭，然后，直接利用FFT子程序，最后，把运算结果取一次共轭，并乘以系数1/N，就可以得到x(n)的结果。
例3.3.1已知有限长序列
x(n)=δ(n)+2δ(n-1)-δ(n-2)+3δ(n-3)
按FFT运算流图求X(k)，再用所得的X(k)按IFFT反求x(n)［7］。
解求DFT的过程如图3.3.7所示，求IDFT的过程如图3.3.8所示。
图3.3.7FFT运算
图3.3.8IFFT运算
X(k)的结果为
X(0)=5，X(1)=2+j，X(2)=-5，X(3)=2-j
3.3.5快速傅里叶变换的应用
1. 快速卷积
若长度为N1的序列x(n)和长度为N2的序列h(n)作线性卷积，得
y(n)=x(n)*h(n)=∑∞m=-∞x(m)h(n-m)
y(n)也是有限长序列，长度为N1+N2-1，此卷积运算需要N1N2次乘法，当N1=N2=N时，需要N2次乘法。
如果用循环卷积实现线性卷积，需要将两序列补零加长至N1+N2-1，这时，利用FFT技术可以大大减少求卷积所需要的运算工作量，这种快速卷积运算如图3.3.9所示。由图可见，在快速卷积中需要两次FFT、一次IFFT运算。在一般的数字滤波器中，由h(n)求 H(k)是预先计算好的，故实际只需要两次FFT运算。若N1=N2=N，所需乘法运算次数为
2×(N/2)×M+N=N(M+1)
显然，N 值越大，式中乘法运算次数比N2越小。
图3.3.9快速卷积
快速卷积可以实现信号的实时处理。但是，在工程实际中，有时遇到要处理的信号很长，对于这类信号只能采用分段卷积的方法。
一般代表滤波器特性的h(n)是有限长序列，其长度为 N，信号x(n)的长度 N1很大，且N1N，将N1分成若干小段，每段长M，以xi(n)表示第i小段。为完成xi(n)和h(n)之间的循环卷积，将xi(n)补零，使其长度达到N+M-1，输入序列为
x(n)=∑mi=0xi(n)
其中，
xi(n)=x(n)，iM≤n≤(i+1)M-1
0，其他，m=N1M
输出序列为
y(n)=x(n)*h(n)=∑mi=0xi(n)*h(n)=∑mi=0［xi(n)*h(n)］=∑mi=0yi(n)
其中，yi(n)=xi(n)*h(n)。
由于yi(n)的长度为N+M-1，xi(n)的非零值长度为M，故相邻的 xi(n) 必有 N-1 长度的重叠。将 yi(n) 求和得 y(n)，其重叠部分必然相加，这种分段卷积再相加求和的方法称为重叠相加法［7］，如图3.3.10所示。
图3.3.10重叠相加法
2. 快速相关
快速相关的原理和快速卷积类似，也是借助于FFT技术实现。相关运算通常用来确定隐含在可加性噪声中的信号，在时域分析中将进一步讨论。利用相关计算还可以求序列的功率谱。快速相关的实现如图3.3.11所示。
图3.3.11快速相关
3.4与本章内容有关的MATLAB函数
本节介绍与傅里叶变换有关的MATLAB函数［8］。
1. fftfilt
该函数用DFT实现长序列的卷积，采用重叠相加法，其调用格式为
y=fftfilt(h,x)
其中，x(n)的长度为N； h(n)的长度为M； 将x(n)分成L段，程序自动确定对x(n) 分段的长度。
2. fft
该函数用来实现快速傅里叶变换，其调用格式为
X=fft(x)或X=fft(x,N)
对前者，若x的长度是2的整数次幂，则按长度实现x的快速变换； 对后者，N应为2的整数次幂，若x的长度小于N，则补零，若超过N，则舍弃N以后的数据。
3. ifft
该函数用来实现快速傅里叶逆变换，调用格式同fft。
例3.4.1已知信号x(t)=2sin(
5πt)+5cos(18πt)，求N点DFT的幅值谱和相位谱(N=64)。
解MATLAB程序如下：
%examp3_1.m
N=64;%设置采样点数
fs=100;%设置采样频率
dt=1/fs;
n=0:N-1;
x=2*sin(5*pi*n*dt)+5*cos(18*pi*n*dt);
y=fft(x,N);%N点傅里叶变换
mag=2*abs(y)/N;%计算信号幅值，abs函数用来求信号的模
pha=angle(y);% angle函数用来求相角
f=n*fs/N;%计算频率，fs/N为频率间隔
subplot(121);%设置绘图窗口，在一幅图中产生两个窗口1×2，在第一个窗口画图
plot(f,mag);
title('Magnitude')
subplot(122);%在第二个窗口画图
plot(f,pha);
title('Phase')
N点DFT的幅值谱和相位谱如图3.4.1所示。
图3.4.1幅值谱与相位谱
小结
本章主要介绍离散时间信号分析的重要工具——傅里叶变换。首先回顾了连续时间信号的傅里叶级数和傅里叶变换，然后，从序列的傅里叶变换入手，重点阐述离散傅里叶变换的产生及物理意义，详细介绍了离散傅里叶变换的一些性质，如线性、对称性、循环移位、循环卷积等，利用循环卷积性质可简化线性卷积的运算； 进一步分析了DFT在实际应用中存在的问题，给出了解决的办法。由于离散傅里叶变换计算量很大，难以实时处理信号，其快速算法——FFT能大大提高运算速度，因此，本章又介绍了时间抽取基2FFT算法和频率抽取基2FFT算法的原理、基本实现及FFT算法的应用； 最后介绍相关的MATLAB函数，给出相应的例程。
习题和上机练习
3.1求序列
h(n)=αn，0≤n<N
0，其他，x(n)=β n-n0，n0≤n
0，n0>n
的卷积。
3.2序列x(n)的傅里叶变换为 X(ejω)，求下列各序列的傅里叶变换：
(1) ejω0nx(n)； (2) nx(n)； (3) x(-n)； (4) x*(n)；
(5) x(n-k)； (6) x2(n)； (7) jIm［x(n)］。
3.3计算下列信号的傅里叶变换：
(1) 2nu(-n)； (2) a|n|u(n)sin(ω0n),|a|<1；
(3) 12n［u(n+3)-u(n-2)］；
(4) cos(18πn/7)+sin(2n)； (5) ∑∞k=014nδ(n-3k)。
3.4已知周期序列
xp(n)=10，2≤n≤6
0，n=0,1,7,8,9
周期N=10，试求Xp(k)=DFS［xp(n)］，并画出Xp(k)的幅度和相位特性。
3.5已知序列
x(n)=an，0≤n≤9
0，其他
分别求其10点和20点离散傅里叶变换。
3.6已知有限长序列x(n)的DFT为X(k)，试用频移性质求序列x(n)sin(2πrn/N)的DFT。
3.7对于有限长序列x(n)，若X(k)=DFT［x(n)］，试证明：
(1) 若x(n)满足x(n)=-x(N-1-n)，则X(0)=0；
(2) 若N为偶数，且有x(n)=x(N-1-n)，则XN2=0。
3.8已知有限长序列x(n)和h(n)如图题3.8所示，试画出：
(1) x(n)和h(n)的线性卷积；
(2) x(n)和h(n)的5点循环卷积；
(3) x(n)和h(n)的8点循环卷积。
图题3.8
3.9设序列x(n)的DFT为X(k)，将它分解为实部和虚部，即
X(k)=XR(k)+jXI(k)
证明:
(1) 若序列x(n)是实序列，则XR(k)是偶函数，XI(k)是奇函数;
(2) 若序列x(n)是纯虚序列，则XR(k)是奇函数，XI(k)是偶函数。
3.10设N点序列x(n)的DFT为X(k)，再按k对X(k)作DFT运算，得到
x1(n)=∑N-1k=0X(k)WknN
试求x1(n)与x(n)的关系。
3.11已知x(n)是长度为N的有限长序列，X(k)=DFT［x(n)］，现将长度扩大r倍(补零增长)，得到长度为rN的有限长序列
y(n)=x(n),0≤n≤N-1
0,N≤n≤rN-1
求DFT［y(n)］与 X(k)的关系。
3.12已知序列x(n)=anu(n)，0<a<1，令对其z变换X(z)在单位圆上N等分点采样，采样值为 X(k)=X(z)|z=W-kN，求有限长序列 IDFT［X(k)］。
3.13已知有限长序列x(n)=δ(n-2)+3δ(n-4)。
(1) 求它的8点离散傅里叶变换X(k);
(2) 已知序列y(n)的8点离散傅里叶变换Y(k)=W4k8X(k)，求序列y(n);
(3) 已知序列m(n)的8点离散傅里叶变换M(k)=X(k)Y(k)，求序列m(n)。
3.14在离散傅里叶变换中产生泄漏和混叠效应的原因是什么？怎样才能减小这种效应？
3.15简略推导按时间抽取基2FFT算法的蝶形公式，并画出N=8时的算法流图，说明该算法的同址运算特点。
*3.16已知两序列
x(n)=(0.9)n,0≤n≤16
0，其他，
h(n)=1,0≤n≤8
0,其他
编写程序以实现序列的线性卷积和N点循环卷积。
*3.17对下面信号进行频谱分析，求幅度谱|X(k)| 和相位谱θ(k)。
(1) x1(t)=at，a=0.8，0≤t≤4ms，fmax=400Hz；
(2) x2(t)=sin t/t，T=0.125s，N=16。
*3.18给定信号x(t)=sin(2πf1t)+2sin(2πf2t)，f1=15Hz，f2=18Hz，现在对x(t)采样，采样点数N=16，采样频率fs=50Hz，设采样序列为x(n)。
(1) 编写程序计算x(n)的频谱，并绘图；
(2) 改变采样频率，得到序列x1(n)，计算x1(n)的频谱，并绘图；
(3) 增大采样点数，得到序列x2(n)，计算x2(n)的频谱，并绘图；
(4) 采样点数N=64，采样频率fs=300Hz，在采样点后补零得到新序列x3(n)，计算x3(n)的频谱，并绘图。
参考文献
［1］胡广书.数字信号处理理论、算法与实现［M］.2版.北京： 清华大学出版社，2003.
［2］徐科军，全书海，王建华.信号处理技术［M］.武汉： 武汉理工大学出版社，2001.
［3］程佩青.数字信号处理教程［M］.2版.北京： 清华大学出版社，2002.
［4］ Oppenheim A V,Schafer R W. DiscreteTime Signal Processing［M］.3rd ed.Publishing House of Electronics Industry,2011.
［5］Oppenheim A V,Willsky A S.Signal and Systems［M］.2nd ed.Prentice Hall,1997.
［6］何振亚. 数字信号处理的理论与应用(上)［M］. 北京： 人民邮电出版社，1983.
［7］靳希，杨尔滨，赵玲.信号处理原理与应用［M］.北京： 清华大学出版社，2004.
［8］楼顺天，李博菡.基于MATLAB的系统分析与设计——信号处理［M］.西安： 西安电子科技大学出版社,1998.
［9］姚天任.数字信号处理学习指导与题解［M］.武汉： 华中科技大学出版社，2002.
［10］谢红梅，赵健. 数字信号处理常见题型解析及模拟试题［M］.西安： 西北工业大学出版社，2001.
［11］程佩青.数字信号处理教程习题分析与解答［M］.2版.北京： 清华大学出版社，2002.
［12］高西全,丁玉美.数字信号处理学习指导［M］.2版.西安： 西安电子科技大学出版社，2001.
［13］徐科军.以能力为导向讲授研究生信号处理课［J］.电气电子教学学报，2019，41(2)： 5659，76.