第3章
CHAPTER 3


离散傅里叶变换












3.1〓DFT引入，DFT、IDFT的定义

微课视频



3.1引言

前面介绍了离散时间域信号(序列)的傅里叶变换(DTFT)和z变换。DTFT和z变换是与连续时间信号的傅里叶变换(CTFT)和拉普拉斯(Laplace)变换相对应的，是对绝对可和、任意长序列进行谱分析和变换域分析的数学工具。


回顾在时间域具有不同特征的信号的谱分析及其频域表述，习惯上称对连续时间域信号x(t)的频谱分析计算为傅里叶变换，如果时间连续信号又是周期的，其频域特征常用傅里叶级数(Fourier series,FS)描述； 而对于离散时间信号x(n)的频谱分析计算则称为序列的傅里叶变换，又称为离散时间傅里叶变换，当离散时间信号又具有周期性时，其频域特征则常以离散傅里叶级数(DFS)描述。

傅里叶变换处理的连续时间域信号及所得到的模拟域频谱函数的特点如图3.1.1所示。其中，信号时间参数t和角频率参数Ω均属于模拟域，是连续变量，在实数域取值，有对应的物理量纲。



图3.1.1模拟域连续时间信号及其傅里叶变换波形示意


连续时间周期信号的时域、频域波形特点如图3.1.2所示。设x(t)的周期为T0，则对应的基波角频率Ω0=2πT0(单位为rad/s)是确定的实数。图3.1.2中时间参数t可在实数域上连续取值； 频谱函数又称傅里叶级数的系数,是线谱,对应的频点为kΩ0，k为整数。这里时间和频率参数也可理解为模拟域参数。



图3.1.2连续时间周期信号的时域、频域波形特点


序列的傅里叶变换(DTFT)对应的信号序列时域、频域特征如图3.1.3所示。由图3.1.3可见，主要有如下特点。



图3.1.3非周期序列及其傅里叶变换

注： DTFT以数字频率为参数时，其周期为2π


(1) DTFT分析对欲处理的时域序列x(n)的长度无限制,一般针对的是非周期序列。参数T是时域序列样值点之间的时间间隔。如果x(n)是由连续函数采样得到的，采样频率为fs(单位为Hz)，则有T=1fs。

(2) 连续时间信号x(nT)的傅里叶变换(FT)结果X(jΩ)是模拟角频率Ω(单位为rad/s)的连续函数，频率变量Ω是连续的； 而序列x(n)的DTFT分析结果X(ejω)一般以数字频率ω=ΩT为变量，需注意频率变量ω(单位为rad)也是连续的。

由


e-j(ω+2π)n=e-jΩ+2πTTn=e-j2πne-jΩTn=e-jωn


得：


X(ej(ω+2π))= X(ejω)
XjΩ+2πT=X(jΩ)



因此，序列x(n)的DTFT结果X(ejω)是数字频率的周期函数，周期为2π； 而连续时间信号x(nT)的FT结果X(jΩ)是模拟角频率的周期函数，周期为2πT。


(3) 序列x(n)的DTFT计算在定义上要求做无限项的求和，实际数值计算总是取有限长的截断序列，这使计算结果相对理论结果会产生误差。

当离散时间序列又具有周期性时，其时域波形特点如图3.1.4(a)所示。因为周期序列的周期特性，其时域序列可展开成离散傅里叶级数(DFS)，从周期序列的DFS展开式可见，周期序列的频谱函数定义在无限多个频点上，其频谱结构是周期线谱； 各离散频点的幅度为|X～(k)|，如图3.1.4(b)所示。



图3.1.4离散时间周期序列及其离散傅里叶级数的系数


本章将介绍便于在计算机上实现的，可对有限长的离散时间信号进行谱分析的新的变换算法： 离散傅里叶变换(DFT)。离散傅里叶变换的特点是： 

(1) 对时域有限长的信号序列,用有限长的频谱序列(DFT的结果)，估计信号序列的频率成分。

(2) DFT的结果是离散数字频率的函数，隐含有周期性。

(3) DFT的快速实现算法统称为快速傅里叶变换(FFT)，是数字信号处理的基础。

离散傅里叶变换(DFT)及其快速算法(FFT)主要有三个用处: 

(1) 离散时间信号(序列)频谱的数值计算。计算机只能对有限的离散数据进行处理，处理结果也是离散数值。

(2) 通过FFT能够实现高效卷积。

(3) 波形编码，尤其是利用DFT的变形——离散余弦变换(discrete cosine transform,DCT)实现。

总结： 在信号处理过程中，如傅里叶变换，常涉及信号在时域和频域间的相互转换； 若一个域中的函数是周期的,则其相应的变换式在另一个域中是采样形式(离散的);反过来,若一个域中的函数是采样形式(离散的),其相应的变换式在另一个域中一定是周期的(可理解为周期延拓的结果)，如表3.1.1所示。


表3.1.1时域和频域的波形特征的对应关系



类型函 数 性 质

时域函数连续的非周期具有周期(T)采样形式(样值点间隔为Ts,对应的采样率为fs)
频域函数非周期连续的离散的线谱，角频率间隔Ω=2π/T具有周期性

(模拟频率域周期为fs； 数字频率域周期为2π)



本章还将介绍DFT和其他信号类型的傅里叶变换相关联时所应满足的条件，最后对用DFT分析现实中经常遇到的连续时间信号和离散时间信号的频谱时遇到的主要问题(混叠,频率泄漏,栅栏效应)及解决思路进行归纳总结。






3.2.1〓N点DFT的特点
微课视频


3.2离散傅里叶变换的定义
3.2.1离散傅里叶变换和逆离散傅里叶变换的定义


设离散时间信号x(n)，n=0,1,…,L-1,是长度为L的有限长序列,其N(N≥L)点离散傅里叶变换(DFT)定义为


X(k)≡DFT［x(n)］=∑N-1n=0x(n)e-j2πNkn(3.2.1a)


定义N点DFT的变换因子为WN=e-j2πN，则式(3.2.1a)可写为


X(k) = ∑N-1n=0x(n) WknN(3.2.1b)


一般情况下，可假设序列x(n)长度区间外的序列值为零，则式(3.2.1b)的求和具体为


X(k)=∑L-1n=0x(n)WknN(3.2.1c)
≡|X(k)|ejθ(k),k=0,1,2,…,N-1(3.2.1d)


其中，N称为离散傅里叶变换的区间长度或DFT窗长度，习惯上称为DFT的变换点数； |X(k)|、θ(k)分别称为序列x(n)的幅度谱和相位谱。显然DFT是把时域有限长序列x(n)变换成数字频率域的有限长序列X(k)。

对X(k)进行逆离散傅里叶变换(inverse discrete Fourier transform,IDFT)，可得到对应的时域有限长序列x(n)，IDFT定义为


x(n)=IDFT［X(k)］
=1N∑N-1k=0 X(k)ej2πNkn
=1N∑N-1k=0X(k)W-knN,n=0,1,2,…,L-1,…,N-1(3.2.2)


由式(3.2.2)可见，N点DFT的逆变换得到的时域有限长序列是N点的，而原时域序列长是L点的。须特别注意，当N>L时，从理论上讲，IDFT求出的x(n)在n=0,1,…,L-1的区间上等于原有限长时域序列； 而在L-1<n≤N-1的区间上，x(n)序列值等于0。实际因计算机的运算误差，N点IDFT算出的N点时间序列可能和原输入时域序列有误差。

另外，在保证从IDFT结果中能唯一地恢复原时域序列的前提下，为了减少计算量,对长为L的序列x(n)进行DFT分析时,通常取变换区间长度为N=L的DFT就可以。

实际应用中有趣的一点是： 从DFT、IDFT的数学定义上，根据WN=e-j2πN因子的特性WN*=ej2πN=WN-1，计算序列x(n)的N点DFT及计算N点IDFT可以共享同一个已有的DFT算法程序实现。说明如下。 

(1) 设通过已有算法程序，得到序列x(n)的N点DFT X(k)：

X(k)=∑N-1n=0x(n)WknN，k=0,1,…,N-1

(2) 根据X(k)，计算N点IDFT，求时域序列x(n)：

 将DFT的定义式X(k)取共轭，得

X*(k) =∑N-1n=0x*(n)W-knN，k=0,1,…,N-1

 对X*(k)计算N点DFT,得

x′(n)=∑N-1k=0X*(k)WknN，n=0,1,…,N-1

 对结果再取共轭: 

x′(n)*=∑N-1k=0X(k)W-knN，n=0,1,…,N-1

 对结果乘以1N，得到序列:

x(n)=1Nx′(n)*=1N∑N-1k=0X(k)W-knN，n=0,1,…,N-1

恰好是式(3.2.2)给出的IDFT的定义式。根据这个特点,编一个DFT子程序,完成DFT运算； 只要再配合两次共轭运算，就能通过DFT实现IDFT的运算。

从DFT的定义可知，序列x(n)的长度L和DFT变换区间长度N可以不同。N>L的DFT计算，数学上可理解为原序列x(n)后面补N-L个零的序列的N点DFT。显然，序列x(n)后补零不影响其N点DFT的结果。证明如下。 

证明: 设x1(n)长为L,n=0,1,…，L-1； x′1(n)是后面增补N-L个零后的序列,长N,n=0,1,…，L-1,…,N-1，如图3.2.1所示。它们的N点DFT分别为


X1(k)=DFT［x1(n)］
=∑L-1n=0x1(n)WknN,k=0,1,…，N-1
X′1(k)=DFT［x′1(n)］=∑N-1n=0x′1(n)WknN
=∑L-1n=0x1(n)WknN,k=0,1,…，N-1


显然: X′1(k)=X1(k)。

通过上面的证明，我们已了解，序列x(n)及其后补零序列进行相同变换区间长度的DFT的结果是一样的； 但序列x(n)及其前补零序列进行相同变换区间长度的DFT的结果则是不一样的。下面介绍这两个序列的DFT之间的差异。

对长为L的x(n)序列，前面补D个零后,得到长为L+D的序列，记为g(n)，n=0,1,…,L+D-1。为了比较x(n)和g(n)的DFT之间的关系，取DFT的点数N≥L+D，为计算简单，取N=L+D，如图3.2.2所示。根据定义，x(n)、g(n)的N点DFT分别为


X(k)=∑N-1n=0x(n)WknN
=∑L-1n=0x(n)e-j2πNkn,k=0,1,2,…,N-1(3.2.3)
G(k)=∑N-1n=0g(n)WknN=∑N-1n=Dg(n)e-j2πNkn(3.2.4)





图3.2.1L点序列及其后补零序列




图3.2.2L点序列及其前补零序列




式(3.2.4)取变量代换n=n′+D,则g(n)=x(n′)；  n=D,N-1时，对应n′=0,L-1，故：


G(k)=∑L-1n=0x(n)e-j2πNk(n+D)
=e-j2πNDk∑L-1n=0x(n)e-j2πNkn
=WDkNX(k)，k=0,1,…,N-1(3.2.5)


比较式(3.2.3)和式(3.2.5)可见： x(n)前补D个零的序列g(n)的N点DFT G(k)的相位相对于原x(n)的N点DFT X(k)的相位，叠加了一线性相移-2πNDk； 因为|WDkN|=1,G(k)和X(k)的幅度是相同的。

在同一个工程实践中，当有两个有限长序列x1(n),x2(n)均要进行DFT分析，如果x1(n)长为L,x2(n)长为N,L<N，为了模块化的计算程序的统一应用，DFT分析子程序常按长的序列取变换区间长度，如进行N点DFT。按式(3.2.1a)的定义分别计算：


X1(k)=DFT［x1(n)］
X2(k)=DFT［x2(n)］


这样做带来的方便之处： 各不同信号片断可以用同一个DFT子程序实现处理； 且当要对这些数据片段进行加权组合，产生新的处理对象时，对应的数字频率域的序列(即频谱序列)也遵循同样的加权组合，结果的长度就等于原离散傅里叶变换的点数，也是规格化的。

例3.2.1给定长度L=8的序列x(n)=［1，1，1，1，0，0，0，0］。

(1) 求其N=16点DFT X(k)；

(2) 前补3个零。求其N=16点DFT XD(k)。

解： 例3.2.1的MATLAB程序参考实现见3.8.1节。

运行3.8.1节给出的程序，结果如图3.2.3和图3.2.4所示。

也可根据定义计算：


X(k)=∑7n=0x(n)e-j2π16kn=∑3n=0x(n)e-j2π16kn,k=0,1,2,…,15
XD(k)=e-j2π163kX(k)=W3k16X(k)，k=0,1,…,15


图3.2.3给出了8点序列x(n)及其前补零序列xD(n)的时域波形及其16点DFT的幅度谱，显然对同样点数的DFT，前补零序列和原序列的幅度谱相同。图3.2.4是两个序列相应的DFT之相位谱； 为了清楚地观察两个序列的相位谱关系，图3.2.4还画出了两个序列相位谱的差，可见前补零序列的相位谱相对于原序列的相位谱随数字频率的变化而出现线性附加相位移。



图3.2.3例3.2.1中8点序列及其前补零序列时域及DFT幅度谱波形




图3.2.48点序列x(n)及其前补零序列的16点DFT相位谱及线性相移


在工程实践中，经常会遇到的问题是要对已获得的时域离散信号(序列)进行频谱分析时，先要确定用多少点数的DFT比较合适； 或者考虑，用现成的DSP芯片对给定的数据进行分析是否适当。回答这些问题要考虑多种因素。其中之一是要了解对某一有限长的时域序列，其不同变换区间长度DFT之间存在什么样的联系。

给定时域序列为x(n),0≤n≤N-1，根据定义，其N点DFT：


X(k)=∑N-1n=0x(n)WknN
=∑N-1n=0x(n)e-j2πNkn，k=0,1,…,N-1(3.2.6)


式(3.2.6)表明，N点DFT离散数字频率间隔ωI=2πN。

(1) 考虑x(n)的L=rN点DFT G(k)，r为整数，即L是N的整倍数。

先对x(n)后补零成为L点序列，再根据定义求rN点DFT：


G(k)=∑rN-1n=0x(n)WknrN，k=0,1,2,…,N-1,N,…,rN-1
=∑rN-1n=0x(n) e-j2πrNkn(因为x(n)只有前N点，所以该求和具体如下)
=∑N-1n=0x(n)e-j2πNkrn


当k是r的整数倍时：


G（k）=Xkr(3.2.7)


G(k)的离散数字频率间隔ω′I=2πrN,比原N点DFT离散数字频率间隔ωI=2πN小。所以rN点DFT G(k)谱线数目在2π的数字频率范围内是X(k)的r倍，更加密集，如图3.2.5中z平面单位圆上等间隔相邻点之间的角度差形象地表示了相应点数的DFT的离散数字频率间隔。


(2) 考虑x(n)的L点DFT G(k)，其中L≠rN，不是原DFT变换区间长度的整数倍。

先对x(n)后补零成为L点序列，再根据DFT的定义，当L>N时,则x(n)的L点DFT G(k)的数字频率离散间隔ω′I=2πL。G(k)谱线数目在2π的数字频率范围内比x(n)原来的N点DFT X(k)多L-N根谱线，也更加密集。设z平面单位圆上的点的幅角为数字频率值，如图3.2.6形象地给出了10点DFT的数字频率值在z平面单位圆上的位置。

例3.2.2设长为8的时域序列x(n)=［1，1，1，1，0，0，0，0］。

(1) 求其N=8点DFT；




图3.2.58点及16点DFT中的数字频率值在z平面单位圆上的分布




图3.2.610点DFT中的10个数字频率值在z平面单位圆上的分布





(2) 后补8个零，即r=2，求其N=16点DFT；

(3) 后补3个零，即L=11，求其N=11点DFT。

解： 例3.2.2的MATLAB程序参考实现见3.8.1节。

运行3.8.1节给出的程序，结果如图3.2.7和图3.2.8所示。



图3.2.7例3.2.2中8点序列x(n)的时域波形





图3.2.8例3.2.2序列x(n)的8点、16点和11点DFT谱线的幅度波形比较



也可根据定义式(3.2.1a)计算出x(n)的N=8、16、11点离散傅里叶变换。

(1) 8点DFT： 

X(k)=∑7n=0x(n)e-j2π8kn=∑3n=0x(n)e-jπ4kn,k=0,1,2,…,7

(2) 16点DFT： 

X(k)=∑15n=0x(n)e-j2π16kn= ∑3n=0x(n)e-jπ8kn，k=0,1,…,15

(3) 11点DFT： 

X(k)=∑10n=0x(n)e-j2π11kn=∑3n=0x(n)e-j2π11kn，k=0,1,…,10

图3.2.7是给定的x(n)序列时域波形。图3.2.8画出了N=8、16、11点离散傅里叶变换的幅度分布。由图3.2.8可见，DFT的结果仅出现在数字频率2πNk，k=0,1,…,N-1处； 因16点DFT与8点DFT涉及8个相同的数字频率值，正如式(3.2.7)描述的，这些频率处的序列频谱值相等。






3.2.2〓DFT和z变换的关系
微课视频


3.2.2离散傅里叶变换和z变换的关系

设序列x(n)长M=N时，其z变换和N点DFT根据定义分别为


X(z)=∑N-1n=0x(n) z-n（且设其ROC含z平面上的单位圆）(3.2.8)
X(k)≡DFT［x(n)］=∑N-1n=0x(n)WknN,k=0,1,2,…,N-1(3.2.9)


显然，当z=zk=W-kN=ej2πNk,k=0,1,2,…,N-1,则z变换和DFT的函数值是相等的； 而zk是z平面单位圆上的点，相邻zk点的角度间隔ωI=2πN。另外根据z变换和DTFT的关系可知： 序列x(n)的z变换X(z)在单位圆上(即z=ejω时)的函数值就是序列(系统)的频谱特性(DTFT)。因此，zk处X(z)的函数值就是序列的DTFT，X(ejω)在［0,2π］范围的N点等间隔采样，可记为X(ejkωI)。

显然相邻zk的角度间隔ωI=2πN也就是序列x(n)的DTF X(ejω)在［0,2π］范围上N点等间隔采样的采样频率间隔。

综上所述，当给定时域序列x(n)，可按下列步骤从其z变换求出其N点DFT： 


x(n)→求z变换，X(z)，且其收敛域包含z平面上的单位圆

→令z=zk = W-kN=ej2πNk≡ejkωI，k=0,1,2,…,N-1

→代入X(z)，得到离散频点处的函数值X(ejkωI),其中： ωI=2πN

→ X(k)≡X(ejkωI),k=0,1,2,…,N-1,即为x(n)的N点DFT序列




数字频率间隔ωI=2πN,和对应的模拟频率间距fI的关系如下：


ωI=ΩIT=2πfIT(3.2.10a)
fI=ωI2πT= 1NT=fsN(3.2.10b)


其中，T为采样序列的时间间隔。如果给定系统的采样频率fs，则fs=1T。

通过上面的分析可以推知序列(离散时间信号)的DFT具有如下的物理意义： 

信号x(n)的离散傅里叶变换X(k)是对其z变换X(z)在z平面单位圆上的N点等间隔采样，k=0,1,…,N-1。等效于对x(n)的DTFT X(ejω)在［0,2π］频率范围内的N点等间隔采样X(ejkωI),序号k采样点处的离散数字频率值为： ωk=kωI=k2πN； 
将数字频率ωk隐去其频率单位，简记为频率序号k，对应的就是信号x(n)的DFT函数X(k)。显然当DFT变换区间长度N不同，在0~2π范围内对X(ejω)的采样频率间隔ωI和点数不同，得到的DFT结果自然不同。

总结DFT的主要特点： 

(1) DFT是关于有限长序列的离散傅里叶变换；

(2) DFT只有N个独立复值，是数字频率域上的有限长序列；

(3) DFT隐含有周期性。

这是因为在DFT、IDFT中，x(n)、X(k)均为有限长序列，且WkN=W(k+mN)N,所以序列x(n)的DFT结果X(k)在数字频率域“隐含有”以N为周期的特点，即X(k+mN)=X(k); X(k)的IDFT结果在离散时间域也“隐含有”以N为周期的特点，即x(n+mN)=x(n)。

下面讨论的周期序列与有限长序列的内在联系，可以从另一个角度，理解离散傅里叶变换隐含有周期性的特点。






3.2.3〓周期序列与有限长序列的关系

微课视频


3.2.3周期序列与有限长序列的内在联系

因为周期序列x～(n)及其频谱(DFS)X～(k)都是周期离散序列，所以只要知道它们一个周期内的内容,其他的所有情况就完全确定了。因此把前面介绍的离散傅里叶级数(DFS)和逆离散傅里叶级数(IDFS)的无限长序列的序号取值范围都限定在一个周期内,这时离散时间域周期序列x～(n)可由主值区间上的主值序列表示。

1. 主值区间

周期为N的序列x～(n)的序号取值为n=0,1,2,…,N-1的范围称为主值区间。

2. 主值序列

周期为N的序列x～(n)在主值区间上的函数值称为周期序列的主值序列，常记为x(n)。这时，x～(n)又可称为x(n)的周期延拓。

周期序列x～(n)的主值区间和主值序列x(n)如图3.2.9所示。



图3.2.9周期序列x～(n)的主值区间和主值序列x(n)


3. 周期序列与主值序列的关系

(1) 时域关系： x(n)是周期序列x～(n)的主值序列; x～(n)是x(n)的周期延拓。两者可以相互表示：


x(n)=x～(n),0≤ n≤N-1,N为周期序列的周期
0,其他(3.2.11)
x～(n)=∑∞m=-∞x(mN+n1)(3.2.12a)



其中，n=mN+n1,在(-∞,+∞)的整数区间任意取值，0≤ n1≤N-1。式(3.2.12a)就是周期延拓的数学表示。习惯上常省略等号右边的下标，周期延拓简记为


x～(n)=∑∞m=-∞x(mN+n)，n在(-∞，+∞)整数区间上任意取值(3.2.12b)


(2) 频域关系： 设X(k)是x(n)的N点DFT,X～(k)是周期为N的x～(n)的DFS。则在数字频率域X(k)是周期序列X～(k)的主值序列； X～(k)是X(k)的周期延拓。


4. 余数运算(求N模数)

定义： 余数运算的结果为整数域上任意数n除以N所得的余数n1。

由欧几里得除法可得，任给整数n> N，总存在整数m、n1，满足：


n=mN+n1，m为整数，0≤ n1≤N-1(3.2.13)


反过来，如果式(3.2.13)成立，则n1称为n除以N的余数。该余数运算可用式(3.2.14)表示：


((n))N=n1(3.2.14)


这里((n))N看成余数运算的运算符号。例如，序列x～(n)的周期为9，n为任意整数。其主值序列记为x(n1)，0≤n1≤8：


x～(25)=x((25))9=x(7)
x～(-5)=x((-5))9=x(4)


周期序列x～(n)的函数值可通过模N的余数运算，用其主值序列的相应函数值x(n1)求出。

利用列长为N的单位矩形序列符号RN(n)和余数运算符号，周期为N的周期序列和其主值序列的时域关系可记为紧凑形式：


x(n)= x～(n)RN(n)(3.2.15)


其中，RN(n)是单位矩形序列,0≤n≤N-1,列长为N。


x～(n)=x((n))N=x(n1)，0≤ n1≤N-1,(((n))N是余数运算）(3.2.16)


式(3.2.16)表明： 周期延拓后任意序号n对应的序列值x～(n)可用主值序列的相应值x(n1)求出。类似地，主值序列的DFT X(k)和周期序列的DFS X～(k)在数字频率域有如下的相互关系： 


X(k)=X～(k)RN(k)(3.2.17)


其中，RN(k)是长为N的单位矩形序列，k=0,1,…，N-1。


X～(k)=X((k))N，-∞<k<+∞（((k))N是余数运算）(3.2.18)


式(3.2.17)表明： 有限长序列x(n)的N点DFT X(k)是其周期延拓序列x((n))N的DFS X～(k)的主值序列。

N点DFT从物理上可看成对离散时间傅里叶变换(DTFT)的(频域)采样,采样间隔为数字频率ωI。类似于对连续信号以时间间隔T的采样，会导致信号频谱的周期延拓，并且延拓周期为1T=fs； 在频率域对DTFT以间隔ωI，即进行频率间隔为fI的采样也会引起时域序列的周期化——周期延拓。原有限长序列x(n)经过此周期延拓成为周期序列x～(n)，且周期为1fI=NT，其中N为频率范围0~fs,即2π内对DTFT采样的采样点数； 隐去具体时间间隔T,该周期序列x～(n)的周期可计为N。图3.2.10给出相关变换的关系示意。



图3.2.10序列的DTFT、DFT、IDFT的物理特性示意



反过来,采样后的DTFT频谱函数就是周期离散序列x～(n)的DFS的系数，等于离散数字频率点ωk=kωI=k2πN处的DTFT函数值：


X～(ejω)|ω=(2π/N)k =∑N-1n=0x～(n) e-j2πNkn
≡X～(k),k为任意整数



图3.2.11z平面单位圆上6等分点
的极角





X～(k)是离散数字频率的周期序列。虽然k任意，因为NωI=N2πN= 2π，一个周期内有N点，这N个独立的函数值就是周期序列的离散傅氏级数(DFS)系数X～(k)的主值序列。这N个离散的数字频率点从几何上可看成z平面上单位圆的N个等分点所对应的极角。如图3.2.11所示，N=6,对应的数字频率分别为0，π/3，2π/3，π，4π/3，5π/3。






3.3.1〓DFT性质——线性性质、序列的和、序列的初值
微课视频



3.3离散傅里叶变换的性质
3.3.1线性性质


设x1(n)和x2(n)是两个长度分别为N1，N2的有限长序列，其线性组合为序列：


x3(n)= a x1(n)+b x2(n)


其中，a、b为任意常数，序列x3(n)的长为N=max［N1,N2］。若序列x1(n)、x2(n)和x3(n)的N点DFT分别记为X1(k)、X2(k)和X3(k)，则：


X3(k) = DFT［a x1(n)+b x2(n)］=aX1(k)+bX2(k)(3.3.1)


3.3.2序列的和

对长为N的序列x(n),其N点DFT记为X(k),则式(3.3.2)成立。


∑N-1n=0x(n)=X(k)|k=0(3.3.2)


证明: 由DFT的定义


X(k)=∑N-1n=0x(n) WknN|k=0，当k=0,W0N=1

故


X(0)=∑N-1n=0x(n)


有限长序列x(n)的各采样值和等于其DFT在k=0处的函数值。

3.3.3序列的初值


对长为N的序列x(n),其N点DFT记为X(k),则式(3.3.3)成立。


x(0)=1N∑N-1k=0X(k)(3.3.3)


由逆离散傅里叶变换的定义可以证明式(3.3.3)。






3.3.4〓DFT性质——序列的圆周移位
微课视频


3.3.4序列的圆周移位(循环移位)
1. 圆周移位(循环移位)


考虑在一个直角坐标系中，以原点为圆心画一个圆，将有限长序列x(n),0≤ n ≤ N-1,均匀排列在这个圆周的N等分点上,正水平轴方向与圆周相交的点记为n=0的序列值； 显然，圆周上相邻两序列值间隔2πN的角度。

如果将整个序列x(n)逆时针方向沿圆周旋转(移位) m2πN的角度,其中m>0(若m<0，则整个序列顺时针旋转)； 再从水平正轴方向开始作为旋转后所得序列xm(n),0≤ n ≤ N-1的起始点,则就是原序列移|m|个单位的圆周移位序列。这种在圆周上的序列圆周移位，用数学定义式表示如下：


xm(n)=x((n-m))NRN(n)(3.3.4)


其中，x((n-m))N是周期为N的周期序列，由原长为N的序列x(n)周期延拓并移位m个单位得到。RN(n)为长为N的单位矩形序列，在这里起到取出主值的作用。

从一维时间轴上，如图3.3.1所示，对式(3.3.4)表达的序列循环移位(又称圆周移位，本书两种叫法不再区分)的物理过程进行如下解释： 

(1) 给定有限长序列x(n)；

(2) 将x(n)周期延拓，得到无限长的周期序列x～(n)=x((n))N；

(3) 若右移位m个单位,则m>0，得到移位m个单位的序列：


x～(n-m)=x((n-m))N


图3.3.1中m=-2，表示左移位2个单位: 


x～(n-m)= x～(n+2)=x((n+2))N


(4) 取出主值序列: 


xm(n)=x((n-m))NRN(n)


x(n)循环移|m|单位的结果是与原x(n)同样长的有限长序列。

当把有限长序列均匀地排列到一个圆周上时，如图3.3.2所示，又可以从二维几何角度，对序列圆周移位(循环移位)过程进行如下解释:

(1) 将x(n),0≤ n ≤ N-1，均匀排列在N等分的圆周上，图3.3.2(a)中不同的序列值用不同的符号点区别。正水平轴方向排列n=0的序列值。圆周上相邻两序列值间隔角度2πN；



(2) 将整个x(n)序列值逆时针方向沿圆周旋转(移位) m2πN的角度，m>0。图3.3.2(b)中m=2，如果是顺时针旋转，则对应于m<0；

(3) 从正水平方向开始逆时针地赋予序列值的序号n,0≤ n ≤N-1,得到x(n)的圆周移位m个间隔的序列，如图3.3.2(c)所示。


xm(n)=x((n-m))NRN(n)


从这个过程可见，圆周移m个单位,所得序列相对原序列x(n)，只是顺序的重新定义。

2. 时域圆周移位定理

设序列x(n)的长为N。其N点DFT记为

X(k)= DFT［x(n)］=∑N-1n=0x(n)WknN

x(n)在圆周上移位m个单位的时间移位序列记为


xm(n)=x((n-m))NRN(n)





图3.3.1循环移位过程一维示意图




图3.3.2圆周移位过程二维示意图




则xm(n)的DFT为


Xm(k)= DFT［xm(n)］=WkmNX(k),k=0,1,…,N-1(3.3.5)


其中，WN=e-j2πN。式(3.3.5)表明，时域移m个单位的序列之DFT相对于原序列的DFT有斜率为m2πN的附加线性相位移，附加相位移Δφ(k)=k·m2πN。

证明： 


Xm(k)=DFT［xm(n)］
=∑N-1n=0x((n-m))NRN(n)WknN,作变量代换i = n-m
=∑-m+N-1i=-mx(i)Wk(i+m)N


因为x(n)的圆周移位均匀地分布在圆周上,所以对x(n)在［-m,-m+N-1］范围上求和等同于对x(n)在［0,N-1］范围上的求和。所以


Xm(k)=WkmN∑N-1i=0x(i)WkiN
=WkmNX(k)
=e-j2πNkmX(k)


3. 频域圆周移位定理——调制特性

如果给定长为N的序列x(n)的离散傅里叶变换X(k)，则其IDFT为时域序列x(n)：


x(n)=IDFT［X(k)］ =1N∑N-1k=0X(k)W-knN


设X(k)在数字频率域的圆周移位为Xm(k)=X((k-m))NRN(k),则其IDFT为时域序列：


x′(n)=IDFT［Xm(k)］ =W-nmNx(n)
=ej2πNnmx(n)(3.3.6)


式(3.3.6)的证明类似于时域圆周移位定理，可用IDFT的定义直接证明。

根据欧拉公式： ej2πNnm=cos2πNm n+jsin2πNm n。当数字系统的采样频率为fs时，数字频率2πNm=2πf/fs,即对应于模拟频率f=fsN m。式(3.3.6)右边表达式物理上可解释为用时域序列x(n)对频率为fsN m的载波进行幅度调制。

根据幅度调制的特点，序列x(n)的频谱将被搬移到“载波频率”Ω=2πfsN m附近。故频域圆周移位定理具有“调制特性”的物理意义是： 把序列x(n)的DFT X(k)在数字频率域上进行圆周移位m个单位得到Xm(k)，对应的是原序列x(n)在离散时间域上的“幅度调制序列ej2πNnmx(n)”的谱。






3.3.5〓序列的圆周卷积

微课视频


3.3.5圆周卷积(循环卷积)及其与有限长序列线性卷积的关系
1. 时域圆周卷积定理

设x1(n),x2(n)长分别为N1,N2;它们的N（即max［N1,N2］）点DFT分别记为X1(k)、X2(k)。若X3(k)=X1(k)X2(k)，k=0,1,…,N-1，则


x3(n)=IDFT［X3(k)］,n=0,1,…,N-1
= ∑N-1m=0x1(m)x2((n-m))NRN(n)(3.3.7a)
≡x1(n)○N x2(n)

或

x3(n)=∑N-1m=0x2(m)x1((n-m))NRN(n)(3.3.7b)
≡x2(n)○N x1(n)



其中，x2((n-m))NRN(n)、x1((n-m))NRN(n)分别是x2(n)、x1(n)的N点圆周移位序列，x1(n)、x2(n)序列中较短的序列要先后补零成N长，所以圆周移位序列的长为N。符号○N 表示N点圆周卷积运算，也可以用符号○* 表示，以区别线性卷积常用运算符*。

式(3.3.7)表示的运算定义为两个序列x1(n)、x2(n)的N点圆周卷积，也称为N点循环卷积(简称圆周卷积或循环卷积)。

显然，因式(3.3.7)中序列x3(n)长为N，故两个序列x1(n)、x2(n)的N点圆周卷积的结果是长为N的序列； 一般地，设两个序列长分别为N1、N2，则两者的圆周卷积结果至少是长为N（即max［N1,N2］）的序列。圆周卷积的具体计算可由式(3.3.7)给出的定义进行，更形象的图解方法将在下面通过例题给出。

时域圆周卷积定理表明： 两时域序列的DFT之积，是离散时间域上这两个序列的圆周卷积的DFT。

2. 频域循环卷积定理

设x1(n)、x2(n)的N点DFT为X1(k)、X2(k)，k=0,1,…,N-1。如果x3(n) = x1(n) x2(n)，其N点DFT为


X3(k)=DFT［x3(n)］ ≡1NX1(k)○N X2(k)
=1N∑N-1j=0X1(j) X2((k-j))NRN(k)(3.3.8a)
=1N∑N-1j=0X2(j)X1((k-j))NRN(k)(3.3.8b)
≡1NX2(k)○N X1(k)


其中，X1((k-j))NRN(k)、X2((k-j))NRN(k)是X1(k)、X2(k)的圆周移位。式(3.3.8)就是频域循环(圆周)卷积定理。

频域圆周卷积定理表明： 离散时间域上两序列之积的DFT，对应着离散数字频率域上这两个序列的DFT之圆周卷积。

3. 循环卷积的计算

不管是时域序列的圆周卷积还是频域序列的圆周卷积，计算方法是一样的。下面根据式(3.3.7)表示的圆周卷积的定义，从表达式的物理意义出发，通过例子说明N点圆周卷积的图解计算方法。

例3.3.1设x1(n)=［1,1,1,0,0,0,0,0］; x2(n)=［1,1,0,0,0,0,0,1］; 求其8点圆周卷积，x3(n)= x1(n)○N x2(n)。

解: 因为x1(n)、x2(n)长度均为8，恰好与要求的圆周卷积的最少点数相等，均不再需要进行后补零处理。

图解法计算圆周卷积的第一步： 将x1(n)、x2(n)按逆时针方向均匀地排列到单位圆的8个等分点上。序号n换成m，正水平方向为m=0的序列值，以便区别最终圆周卷积结果的顺序号n，如图3.3.3(d)和图3.3.3(b)所示。针对式(3.3.7a)和式(3.3.7b)两个定义式，给出下面两个图解方法流程，两个方法的计算结果是完全等效的。

1） 方法1

(1) 将x1(m)循环反序： 等效于有限长序列在圆周上反序，即原序列相对水平对称轴对折，简记为x1(-m)，如图3.3.3(a)所示，原序列相对水平轴的对称点序列值互换。

(2) 循环反序序列x1(-m)逆时针移位n个单位，图3.3.3(a)表示n=0的情况。

(3) 步骤(2)得到的序列与如图3.3.3(b)所示序列x2(m)的对应点函数值两两相乘，并将8个积求和，得到x3(n)：


x3(n)=x2(n)○N x1(n)
=∑7m=0x2(m)x1((n-m))8R8(n)



(4) 令n=0,1,…,7,重复步骤(2)和步骤(3)。得到完整的圆周卷积结果，如图3.3.3(c)所示，图3.3.3(c)括号中的数值为该顺序点的序列值，正水平方向为n=0的序列值。表3.3.1也给出了x1(n)和x2(n)的8点循环卷积的结果。



图3.3.3序列(圆周卷积)循环卷积过程图解




表3.3.1例3.3.1的循环卷积结果



n01234567
x3(n)23210001



2） 方法2

(1) 由于序列的循环卷积具有可交换性，所以也可以把x2(m)循环反序，如图3.3.3(e)所示。

(2) 将x2(-m)逆时针移位n个单位，图3.3.3(e)表示n=0的情况。

(3) 将步骤(2)得到的序列与图3.3.3(d)所示的x1(m)序列对应顺序点序列值两两相乘并将8个积相加，得x3(n)：


x3(n)= x1(n)○N x2(n)=∑7m=0x1(m)x2((n-m))8R8(n)


(4) 令n=0,1,…,7,重复步骤(2)和步骤(3)得到最终的圆周卷积序列。结果与方法1的完全一样。


4. 循环卷积和线性卷积的关系


在实际数字信号处理中，(数字信号)序列x(n)通过单位采样脉冲响应为h(n)的数字系统，则输出y(n)是输入序列和系统单位采样脉冲响应的线性卷积，如图3.3.4所示。



图3.3.4序列通过单位采样脉冲
响应为h(n)的数字系统




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




线性卷积能否用运算更快的圆周卷积实现?回答是肯定的，但需要在一定的条件下才能实现。

首先，考虑到两列长为N的序列，它们的循环卷积结果仍是长为N的序列，而它们的线性卷积结果是长为2N-1的序列。更一般地，给定两列长分别为M、N的序列,它们的线性卷积结果是长为N+M-1的序列,它们之间圆周卷积结果的长为max(N,M)。

为了能用循环卷积代替线性卷积，就要保证循环卷积的结果恰好是期望的线性卷积结果； 显然，圆周卷积结果至少应该和线性卷积结果一样长。为此，按下列方法计算循环卷积的点数L，就可以通过L点的圆周卷积实现两序列的线性卷积。

设序列x(n)、h(n)的长分别为M、N，则知其线性卷积y(n)是长为N+M-1点的序列。

一般地，取L≥N+M-1。将x(n)和h(n)后面分别补L-M、L-N个零,使其长均为L点； 对这两个后补零后得到的新的长为L的序列进行L点的圆周卷积，则结果是长L点的序列。可以证明，这L点圆周卷积结果的前N+M-1点序列值即等于原两序列的线性卷积结果。

若取圆周卷积的点数L小于线性卷积结果长N+M-1，则由于圆周卷积隐含的周期延拓特性,两序列的L点圆周卷积结果可看成是周期为L(小于M+N-1)的周期序列的主值序列，圆周卷积结果不能正确地表示线性卷积结果。

例3.3.2设x1(n)和x2(n)是两个4点序列： 


x1(n)={1,2,2,1},x2(n)={1,-1,1,-1},n=0,1,2,3


(1) 确定它们的线性卷积x3(n)；

(2) 计算补零后的7点循环卷积x4(n)；

(3) 讨论不补零的4点循环卷积和其他补零数目的循环卷积结果之间的区别。

解： 例3.3.2的MATLAB程序参考实现见3.8.2节。

根据给定的两个序列，可画出x1(n)和x2(n)的时域波形，如图3.3.5所示。



图3.3.5例3.3.2的两时域序列


(1) x1(n)和x2(n)的线性卷积结果应长： 4+4-1=7点。运行3.8.2节的有关MATLAB参考程序，得到x1(n)和x2(n)的线性卷积结果。图3.3.6中x3画出了线性卷积结果波形x3(n)：


x3(n)={1,1,1,0,-1,-1,-1},n=0,1,…，6


(2) 依题目给定条件知，两序列的循环卷积至少长max［4,4］=4点，但为保证循环卷积与线性卷积结果一致，应将两序列后补至少7-4=3个零，再进行延长序列的循环卷积，即循环卷积点数N≥7。

运行3.8.2节的有关MATLAB参考程序结果如图3.3.6所示，其中x4、x5分别画出了x1(n)和x2(n)的7点、6点圆周卷积结果波形x4(n)、x5(n)：


x4(n)={1,1,1,0,-1,-1,-1},n=0,1,…，6
x5(n)={0,1,1,0,-1,-1},n=0,1,…，5




图3.3.6例3.3.2线性卷积和圆周卷积的比较


(3) 讨论。
比较图3.3.6中x4(n)、x3(n)波形可见，7点圆周卷积完全等于给定序列的线性卷积。验证了我们所介绍的，即在给定条件下，可通过圆周卷积实现线性卷积的原则。

再比较图3.3.6中6点圆周卷积结果x5(n)和7点圆周卷积结果x4(n)波形可见，x5(n)的长度为6点，比x4(n)少了1点，且x5(0)=0，而x4(0)=1。

可以证明： 6点圆周卷积结果x5(n)可以通过将7点圆周卷积结果x4(n)以6为周期进行周期延拓，并取延拓后的6点主值序列得到。显然，长为7的序列x4(n)在这个延拓过程中，其首尾序列值将会“混叠”，结果导致x5(0)=x4(0)+ x4(6)=1+(-1)=0，正如图3.3.6所示的波形x5(n)。

根据上述不同点数的圆周卷积结果的关系特点，请读者思考： 本例中，如果直接进行x1(n)和x2(n)的4点圆周卷积，将其结果和这两个序列的线性卷积相比较，4点圆周卷积结果哪些点与线性卷积的结果不同？提示： 把线性卷积结果序列以4为周期进行周期延拓，再取4点的主值序列，就是4点圆周卷积的结果。

5. 圆周卷积的应用

实现快速卷积: 因为有限长圆周卷积的快速计算可以通过DFT的快速实现(FFT)进行,从而，用圆周卷积代替线性卷积，则等效于快速线性卷积。

如图3.3.4所示，设输入序列x(n)长为L,系统单位采样脉冲响应h(n)长为M，系统的时域输出y(n)为线性卷积：


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


可通过计算x(n)、h(n)的N=L+M-1点圆周卷积，快速实现。具体步骤如下： 

(1) 计算x(n)、h(n)的N点DFT(FFT)，X(k),H(k)；

(2) 计算X(k)·H(k)；

(3) 根据时域圆周卷积定理，对步骤(2)的结果进行IDFT即为x(n)与h(n)的圆周卷积x(n)○N h(n)；

(4) 根据圆周卷积与线性卷积的关系，步骤(3)的结果即为系统的时域输出y(n)。

下面讨论快速卷积和直接卷积的运算量差异。根据定义，直接进行x(n)与h(n)的线性卷积，其乘法运算量为


md=LM(3.3.9)


后面将介绍，通过快速傅里叶变换(FFT)实现快速卷积的乘法运算量为 


mF=32Nlog2N+N,N=L+M-1(3.3.10)


当L不太大，且x(n)与h(n)的长近似时，比较适合用式(3.3.10)表示的N点FFT的乘法运算量估计快速卷积的运算量； 当L≈M>64时，随着M的增加，有mFmd。对于L很大，且LM的情况，N=L+M-1将会非常大，这时更常用的快速卷积是分段圆周卷积(如重叠相加,重叠保留法)。这方面的内容将在第4章具体介绍。






3.3.63.3.7〓序列相关性、Parseval定理

微课视频


3.3.6序列的相关性

两N点序列x(n)、y(n)的相关主要有线性相关和圆周相关(也称循环相关)。

1） N点圆周相关定义


rxy(m)=∑N-1n=0x(n)y*((n-m))NRN(m)(3.3.11a)
=∑N-1n=0x((n+m))NRN(m)y*(n)(3.3.11b)


其中，相关延迟m等于x(n)的时间序号减去y*(n-m)的时间序号； x((n+m))NRN(m)、y((n-m))NRN(m)分别是x(n)、y(n)的m单位圆周移位，但序列不需要反褶。

2) 两序列的线性相关定义


rxy(m)=∑∞n=-∞x(n)y*(n-m)(3.3.12a)
=∑∞n=-∞x(n+m)y*(n)(3.3.12b)
= r*yx(-m)(3.3.12c)


其中，m是相关延迟，x(n+m)、y(n-m)分别是x(n)、y(n)的m单位线性移位，也不需要反褶。

3) N点圆周相关与线性相关的关系

式(3.3.12)是计算两个序列线性相关的理论模型，求和限在{-∞,+∞}之间，表示线性相关计算应包括全部序列值的影响。当两个序列x(n)、y(n)为有限长，分别为N1点、N2点时，将序列x(n)、y(n)后补零成N≥N1+N2-1点长后，求出的N点圆周相关序列的前N1+N2-1点等于其线性相关序列。

4) 讨论

(1) 比较式(3.3.7)和式(3.3.11)可知： 圆周卷积要将一个序列循环反褶后再移位计算； 而圆周相关不需循环反褶,只要移位，就可以求出两序列的相关序列。

(2) 定义式(3.3.11)通常称为序列x(n)、y(n)的N点互圆周相关。如果序列的相关是对同一个序列进行的，则习惯上称之为序列x(n) N点自圆周相关。


rxx(m)=∑N-1n=0x(n)x*((n-m))NRN(m)(3.3.13)


当一个较短的序列和一个很长的序列计算圆周相关时，类似于圆周卷积，也存在分段相关的方法。本书不做介绍。

3.3.7帕斯瓦尔定理

帕斯瓦尔定理： 设有限长序列x(n)、y(n)的N点DFT分别为X(k)、Y(k)，则


∑N-1n=0x(n)y*(n)= 1N∑N-1k=0X(k)Y*(k)(3.3.14)


如果y(n)=x(n)，则式(3.3.14)左边是时域有限序列的能量E，右边是序列的谱域能量和，即


E=∑N-1n=0|x(n)|2 =1N∑N-1k=0|X(k)|2(3.3.15)


显然，帕斯瓦尔定理把序列时域能量表达与其频域能量表达联系起来了，因此，有的文献称该定理为“能量定理”。式(3.3.15)也可理解为信号能量守恒的体现。

比较序列自相关定义式(3.3.13)和式(3.3.15)可见，相关间隔为零的自相关序列值：


rxx(m)|m=0=rxx(0)=序列的能量


对有限长序列，式(3.3.15)表示能量E，也即rxx(0)有限，因此有限长序列一定是能量有限信号。这与“通信原理”课程中介绍的零均值平稳随机过程的自相关函数的性质类似，即自相关函数值R(0)为功率有限信号的平均功率。






3.3.8〓有限长序列及其DFT的对称性

微课视频


3.3.8有限长序列及其离散傅里叶变换的奇偶性和对称性

类似于第2章中讨论的DTFT的对称性，有限长序列x(n)及其DFT(离散傅里叶变换)X(k)也可以定义对称性。但是，假设有限长序列长为N，定义区间为［0,N-1］，其N点DFT的列长也为N，则其奇偶性和对称性应理解为关于“中心(N/2)点”的对称性。下面讨论有限长序列及其DFT相关的共轭对称性和奇偶性的定义和性质。

1. 有限长共轭对称序列和共轭反对称序列

设xep(n)，xop(n)分别表示有限长共轭对称序列和共轭反对称序列。则有如下定义式：


xep(n)=x*ep(N-n)，0≤n≤N-1(3.3.16a)
xop(n)=-x*op(N-n)，0≤n≤N-1(3.3.16b)


当N为偶数时，定义式(3.3.16a）和式（3.3.16b)中的N可换成N2-n，得到另一种定义式： 


xepN2-n=x*epN2+n，0≤n≤N2-1(3.3.17a)
xopN2-n=-x*opN2+n，0≤n≤N2-1(3.3.17b)


图3.3.7给出了N=7和N=8两种长度对应的共轭对称和共轭反对称序列示意图，可清楚地说明有限长序列的“中心”对称性。图3.3.7中*表示对应点为序列值取共轭后的值。



图3.3.7共轭对称和共轭反对称序列示意图


性质3.3.1任何有限长序列x(n)均可分解为共轭对称分量和共轭反对称分量之和，即


x(n)=xep(n)+xop(n)，0≤n≤N-1(3.3.18)


若给定序列x(n)，其共轭对称分量xep(n)和共轭反对称分量xop(n)可由如下计算获得。

将式(3.3.18)中的n换成N-n，并取复共轭，再将定义式(3.3.16a）,式（3.3.16b)代入得到


x*(N-n)=x*ep(N-n)+x*op(N-n)
=xep(n)-xop(n)，0≤n≤N-1(3.3.19)


式(3.3.18)分别加、减式(3.3.19)，整理可得


xep(n)=12［x(n)+x*(N-n)］，0≤n≤N-1(3.3.20a)
xop(n)=12［x(n)-x*(N-n)］，0≤n≤N-1(3.3.20b)


注意： 有限长序列的长度N不管是偶数还是奇数，其第1个样本x(0)，在［0,N-1］的定义区间上都没有其对称样本，如图3.3.7所示。

仔细观察定义式（3.3.16a）和式（3.3.16b）可见，当n=0时，有xep(0)＝x*ep(N)，xop(0)＝－x*op(N)； 而长为N的序列定义区间是［0,N－1］， 则xep(N)、xop(N)的取值是多少？这个问题的物理解释如下： 

首先，性质3.3.1中，有限长序列x(n)的共轭对称分量xep(n)和共轭反对称分量xop(n)在有的文献上分别称为“圆周共轭对称分量”和“圆周共轭反对称分量”。

其次，长为N的有限长序列x(n)若理解成某周期为N的周期序列x～(n)的主值序列， 则x(n)的xep(n)分量和xop(n)分量，对应地分别为周期序列x～(n)的共轭对称分量x～e(n)和共轭反对称分量x～o(n)的主值序列，参见第2章2.2.2节“对称性”的内容。因此，从周期为N的周期序列x～e(n)、x～o(n)看，xep(N)=x～e(N)，xop(N)=x～o(N)的取值也就得到了合理解释。

最后，考虑到周期序列x～(n)的共轭对称分量x～e(n)、共轭反对称分量x～o(n)亦是周期序列，故xep(N)、xop(N)数值与其主值序列xep(n)、xop(n)中的xep(0)、xop(0)呈对称关系。

类似地，式(3.3.20a)、式(3.3.20b) 中，n=0时，令序列值x(N)=x～(N)就是自然而然的了。

2. DFT的共轭对称性

根据已经掌握的基础知识，我们知道任何复序列都可以分解为纯实序列和纯虚序列分量之和，也可以分解为共轭对称和共轭反对称序列分量之和。

给定任意有限长复序列x(n)，其DFT记为X(k)。则x(n)有如下分解形式： 


x(n) = xep(n)+xop(n)


即，分解为共轭对称序列和共轭反对称序列之和，由式(3.3.20a）和式（3.3.20b)可计算各分量。或者：


x(n)= xr(n)+jxi(n)


即，分解为实序列xr(n)和纯虚序列jxi(n)之和。其中各分量可由式(3.3.21a）和式(3.3.21b)计算：


xr=Re［x(n)］=12［x(n)+x*(n)］（为x(n)的纯实分量序列）
jxi(n)=jIm［x(n)］=12［x(n)-x*(n)］（为x(n)的纯虚分量序列）(3.3.21a)(3.3.21b)



类似地，x(n)的DFT X(k)在数字频率域也可分解如下： 


X(k)=Xep(k)+Xop(k)（分解为共轭对称序列和共轭反对称序列之和）
=XR(k)+jXI(k)（分解为实序列XR(k)和纯虚序列jXI(k)之和）


实际上，由DFT的线性性和序列的共轭对称、共轭反对称的定义，可以证明：


DFT［xr(n)］ =Xep(k)(3.3.22a)
DFT［jxi(n)］=Xop(k)(3.3.22b)


其中，Xep(k)、Xop(k)分别为x(n)的DFT X(k)的共轭对称、共轭反对称分量序列。对应地：


DFT［xep(n)］ =Re［X(k)］= XR(k)(3.3.23a)
DFT［xop(n)］=jIm［X(k)］=jXI(k)(3.3.23b)


其中，XR(k)、jXI(k)分别为DFT X(k)的纯实、纯虚分量序列。



图3.3.8有限长序列的分解序列与其DFT
序列的分解序列之间的对应关系


图3.3.8给出了有限长复序列x(n)及其分解序列与相应的DFT之间的对应关系。即时域实序列的DFT是频域的共轭对称序列，共轭对称序列的DFT则是频域的实序列； 时域共轭反对称序列的DFT是频域的纯虚序列。性质3.3.2给出了更具体的时域序列的DFT特性。


性质3.3.2设x(n)是长度为N的实序列，则分解式x(n)=xep(n)+xop(n)，其中，0≤n≤N-1中的xep(n)、xop(n)均为实序列，且分别为“中心”偶对称和“中心”奇对称的序列。若记X(k)=DFT［x(n)］，则

(1) X(k)=X*(N-k),0≤k≤N-1（是共轭对称序列）；(3.3.24)

(2) 如果x(n)=x(N-m)，是“中心”实偶对称的，则


X(k)=X(N-k),0≤k≤N-1（是实偶对称序列）(3.3.25)


(3) 如果x(n)=-x(N-n)，是“中心”实奇对称的，则


X(k)=-X(N-k)


且


X(k)=jXI(k)，0≤k≤N-1(3.3.26)


是纯虚奇对称的，其中XI(k)是实序列。

证明： 以性质3.3.2中式(3.3.26)为例。根据DFT的定义：

X(k)=∑N-1n=0x(n)e-j2πNkn

将x(n)=-x(N-n)代入，并进行N-n=m变量代换：

X(k)=-∑N-1n=0x(N-n)e-j2πNkn
=-∑1m=Nx(m)ej2πNkm(3.3.27)


若把x(n)的N个序列值均匀地排列到单位圆上，显然式(3.3.27)右边求和等效为m在［0,N-1］范围内的求和：


X(k)=-∑N-1m=0x(m)ej2πNkm(3.3.28)


又因：


X(N-k)=∑N-1n=0x(n)e-j2πN(N-k)n=∑N-1n=0x(n)ej2πNkn(3.3.29)


比较式(3.3.28)和式(3.3.29)的右边，有


X(k)=-X(N-k)，0≤k≤N-1


将式(3.3.28)两边取共轭，考虑到x(n)是实序列，有


X*(k)=-∑N-1m=0x(m)e-j2πNkm=-X(k)(3.3.30)


仅当X(k)=jXI(k)，XI(k)为实序列时，式(3.3.30)才会成立。

3. 应用

实际应用中，经常要对实序列进行N点DFT计算，利用上面所介绍的离散傅里叶变换的共轭对称性和奇偶对称性可减少运算量，提高计算效率。

(1) 计算N点DFT，当N为偶数时，只要计算前N2+1点； 当N为奇数时，只要计算前N+12点； 其他点按照式(3.3.24)可求,可减少近一半的计算量。

(2) 利用DFT的共轭对称性，计算一个N点DFT X(k)可得两个不同实序列的N点DFT。设xr(n)、xi(n)为两个实序列，下面介绍通过一次N点DFT，得到这两个序列的DFT Xr(k)、Xi(k)，k=0,1,…,N-1。构造新复序列：


x(n)=xr(n)+jxi(n)


对x(n)进行N点DFT，得X(k)=Xep(k)+Xop(k)。由式(3.3.20a）和式（3.3.20b)可求出Xep(k)、Xop(k)分量，再利用式(3.3.22a）和式（3.3.22b)得


12［X(k)+X*(N-k)］=DFT［xr(n)］=Xr(k)，k=0,1,…,N-1
-j12［X(k)-X*(N-k)］=DFT［xi(n)］=Xi(k),k=0,1,…,N-1







3.4〓频域采样定理

微课视频


3.4频域采样

设列长为M的序列x(n)，其z变换为X(z)，其N点DFT为X(k),k=0,1,…,N-1。其N点IDFT为


xN(n)=IDFT［X(k)］，n=0，1，…,N-1(3.4.1)


根据IDFT隐含有周期性的特点，xN(n)可看成周期为N的时域序列的主值序列。数字信号处理中会遇到的问题如下： 

(1) x(n)和xN(n)的关系如何？

(2) X(z)和X(k)的关系如何？

(3) 能否由X(k)重建X(z)?

本节介绍的频域采样及频域采样定理，可以在一定程度上回答和解释上述问题。

性质3.4.1频域采样。

长为M的序列x(n)的z变换X(z)：


X(z)=∑∞n=-∞x(n)z-n

因x(n)长为M，有

X(z)=∑M-1n=0x(n)z-n

当N≥M时，可等效写为



X（z）=∑N-1n=0x(n)z-n(3.4.2)


且X(z)的收敛域包含z平面上的单位圆，即x(n)存在傅里叶变换。

令z=zk=ej2πNk,k=0,1,…,N-1,在z平面单位圆上对X(z)等间隔采样N点，得到: 


X(k)=X(z)|z=zk，k=0,1,…,N-1(3.4.3a)


代入式(3.4.2),得： 


X（k）=∑N-1n=0x(n)e-j2πNkn(3.4.3b)


式(3.4.3a)中右边表示在数字频率区间［0,2π］上，对x(n)的离散时间傅里叶变换(DTFT)X(ejω)的N点等间隔采样，式(3.4.3b)正是x(n)的N点DFT定义式，即x(n)的z变换X(z)在z平面单元圆上的N点等间隔采样就是x(n)的N点DFT序列X(k)。

如果对x(n)的N点DFT序列X(k)进行IDFT，根据IDFT隐含的周期特性，由式(3.4.1)得到的时域序列xN(n)长为N，而原x(n)的长为M。频域采样定理给出了xN(n)和x(n)两者之间的关系。

频域采样定理如果序列x(n)的长度为M，则只有当对其DTFT X(ejω)进行的频域采样点数N≥M时，才有


xN(n)=IDFT［X(k)］=x(n),n=0,1,…，M-1(3.4.4)


式(3.4.4)表示可由序列的DTFT的适当的频域采样序列X(k)无失真地恢复原有限长时域序列x(n)。否则由DTFT的采样序列逆变换得到的时域序列，相对原有限长时域序列有时域混叠现象，会产生失真。

根据DFT与DFS的关系： X(k)是xN(n)以N为周期的周期延拓序列x～(n)的DFS的主值序列。可以证明xN(n)是原序列x(n)以N为周期的周期延拓序列的主值序列。当N<M时,x(n)以N为周期的周期延拓会产生时域混叠现象。

显然，对无限长的时域序列,无论对其DTFT在数字频域的采样点N多大，由频域采样序列X(k)的IDFT恢复的时域序列都不能完全消除时域混叠现象。加大数字频域的采样点N,只能使恢复的信号逼近x(n)，减小误差。

下面给出在满足频域采样定理的条件下，用频域采样序列X(k)恢复原有限长序列x(n)的z变换X(z)的公式。根据有限长序列的z变换和IDFT的定义，具体推导如下。

设序列x(n)的长度为M，在频域0~2π对其z变换等间隔采样N点，N≥M，则有


X(z)=∑N-1n=0x(n)z-n
X(k)=X(z)|z=ej2πNk，k=0,1,…,N-1


式中，


x(n)=IDFT［X(k)］=1N∑N-1k=0X(k)W-knN（这里WN=e-j2πN）


将上面的x(n)代入X(z)的表示式中得：


X(z)=∑N-1n=01N∑N-1k=0X(k)W-knNz-n
=1N∑N-1k=0X(k)∑N-1n=0W-knNz-n
=1N∑N-1k=0X(k)1-W-kNNz-N1-W-kNz-1


式中，W-kNN=1，故：


X(z)=1N∑N-1k=0X(k)1-z-N1-W-kNz-1(3.4.5a)
=1-z-NN∑N-1k=0X(k)1-W-kNz-1(3.4.5b)


令


φk(z)=1N1-z-N1-W-kNz-1(3.4.6)


则


X(z)=∑N-1k=0X(k)φk(z)(3.4.7)


式(3.4.7)称为用有限长时域序列的DFT X(k)表示X(z)的z域内插公式，式(3.4.6)中的φk(z)称为z域内插函数。当z=ejω时，式(3.4.6)和式(3.4.7)就成为在满足频域采样定理条件下，由有限长时域序列x(n)的离散傅里叶变换X(k)恢复x(n)的傅里叶变换X(ejω)的频率域内插函数和内插公式，即


φk(ω)= 1N1-e-jωN1-e-jω-2πNk
=1Nsinω2Nsinω2-πkNe-jN-12ω+πkN
=1Nsinω2-πkNNsinω2-πkNejπkN(N-1)e-jN-12ω
X(ejω)=∑N-1k=0X(k)φk(ω)


进一步化简，得：


X(ejω)=∑N-1k=0X(k)φω-2πNk(3.4.8)


其中，


φ(ω)=1Nsinω2Nsinω2e-jωN-12(3.4.9)


在第5章中，我们将看到，z域X(z)的插值公式(3.4.5b)正是有限长单位采样响应系统的可能实现结构之一——“频率采样型”的理论根据。

3.5用离散傅里叶变换对连续时间信号进行谱分析

第4章将要介绍的DFT的快速算法(FFT)的出现，使DFT在许多领域,如数字通信、音视频信息处理、功率谱估计、系统分析、雷达理论、光学、医学、地震等领域有着广泛应用。在众多的数值分析应用中，时域以卷积运算和相关运算的处理为基础，频域以傅里叶变换类谱分析为基础。本节主要介绍用DFT对连续时间信号进行谱分析方面的最基本的概念和应用。






3.5.1〓DFT对连续时间信号进行谱分析

微课视频


3.5.1离散傅里叶变换是对连续时间信号傅里叶变换的近似

对信号进行谱分析的基本目标是： 得到待分析信号频率特性，即分析信号主要能量(功率)所在的频率范围或可能包含的频率分量，以及所得信号频谱中各分量的相对电平关系。对连续时间信号进行谱分析的基本方法——傅里叶变换(FT)虽然很精确，但其不能充分发挥计算机数值计算的优势。实际中，常用离散傅里叶变换(DFT)代替，近似实现连续时间信号的傅里叶变换。因为离散傅里叶变换可以通过快速算法(FFT)，在计算机或数字信号处理器(DSP)平台上高效、实时地实现。

考虑到DFT是对离散时间有限长序列x(n)定义的，因此对连续时间信号x′(t)应用DFT进行谱分析，必然涉及两方面的处理： 一是要对x′(t)采样，使之成为时间离散的序列； 二是要截断，取得有限长信号。这两个处理使得用DFT对连续时间信号进行谱分析，仅是傅里叶变换分析的近似。两种分析结果的近似程度(误差大小)与待分析信号x′(t)的固有带宽、采样频率、截取长度等物理参数以及DFT的点数等计算参数有关。

设连续时间信号x′(t)是带限的，其频率范围为［0，fh］，信号持续时间t任意,对其截断，截取持续时间总长度为TL的一段，记为x(t)。根据奈奎斯特采样定理，选择采样频率fs≥2fh，对x(t)在时间t轴上等间隔间隔T=1fs采样得序列x(n)，准备进行DFT分析。

这里要注意的问题是： 截取的有限长(TL)信号片段x(t)的频谱结构X(jΩ)是确定的，假设X(jΩ)能反映原始信号x′(t)的频谱结构，对x(t)进行采样，进行DFT分析才是有意义的。对时长TL的信号片段x(t)采样，采样信号序列的频谱是X(jΩ)以周期fs延拓的周期谱，采样频率fs足够高，可保证在奈奎斯特频率间隔-fs2,fs2内采样信号序列x(n)的频谱和截断所得的连续信号x(t)的频谱X(jΩ)一样。

还要注意选择合适点数的DFT进行谱分析，结果才能更好地近似有限长信号x(t)的频谱X(jΩ)。这些问题可通过具体的“物理”频率分辨率和“计算”频率分辨率指标进行定量描述。“频率分辨率”是谱分析质量的基本参数之一，指可分辨信号中频率分量或者频率变化的能力。

3.5.2频率分辨率

数字信号处理中的分辨率可指时间分辨率和频率分辨率，其中频率分辨率又可分为物理频率分辨率(physical frequency resolution)和计算频率分辨率(computational frequency resolution)。形象地讲，频率分辨率是指离散频谱中相邻两点间的频率间隔，即可识别的最小频率变化； 当两相邻离散谱线之一为零频率(直流分量)时，另一谱线对应的频率则是信号中的最小非零频率分量，所以频率分辨率也可理解为“可识别的最小频率分量”。	

1. 物理频率分辨率Δf

设给定的有限长信号x(t)(模拟/数字)的总持续时间长度为TL，如图3.5.1所示。物理频率分辨率Δf定义为


图3.5.1总时长为TL的信号中
可能的最小频率分量fmin



Δf=1TL(3.5.1)



从图3.5.1可见，时长为TL的信号x(t)中可能的信号周期小于或等于TL，表示信号x(t)中可能包含的最小频率分量fmin=1TL。显然这段信号的TL越长，信号能包含的最低频率分量就越小。对应地，对该信号做频谱分析时，

可识别的信号频率分量大于或等于物理频率分辨率Δf，根据式(3.5.1)，物理频率分辨率Δf由拟处理信号的总时长确定； 反过来讲，物理频率分辨率Δf是频谱分析应能分辨的最小频率差。


给定系统采样频率fs后，有限时长TL内信号x(t)可包含的采样点数L=TLfs。因此在给定系统采样频率fs的条件下，物理频率分辨率Δf还可由记录数据的采样序列长L确定：


Δf =1TL =fsL(3.5.2)


2. 计算频率分辨率Δfbin

设已知系统的采样频率fs，对截断得到的有限长采样序列x(n)进行N点DFT，可得到频域序列X(k)。X(k)等于对采样信号x(n)的DTFT频谱X(ejω)在奈奎斯特频率范围［0,fs］内进行频域采样，其中fs2，fs范围上的频谱对应于原DTFT频谱的负频率范围-fs2,0的频谱水平右移fs。从谱分析的物理意义上讲，DFT处理有限长连续时间信号x(t)的本质之一是用数字频率的频谱序列值X(k)近似原信号x(t)的连续频谱函数X(jΩ)。

因为在时域对连续时间信号x(t)采样后，x(t)的连续频谱函数X(jΩ)在模拟频率域以fs为周期进行周期延拓就是采用序列x(n)的频谱函数； 对应到数字频率域，即为周期为2π的频谱函数X(ejω)，所以N点DFT计算能反映的信号最小频率变化(频域采样间隔)为fsN(单位为Hz)或2πN。因此，N点DFT对应的计算频率分辨率Δfbin定义为


Δfbin=fsN(3.5.3)


3. 物理频率分辨率Δf和计算频率分辨率Δfbin的比较

对任意长的连续时间信号x′(t)的TL长的时域记录x(t)以采样频率fs进行采样得x(n)。这种截断处理对应的物理频率分辨率Δf根据式(3.5.1)的定义，有如下特点。

(1) 取决于实际有限长信号(采样)序列x(n)的参数。时域序列一定，Δf就是确定的。

(2) x(n)的固有频谱(DTFT)只能区分出x′(t)的频率分量中频差大于或等于物理频率分辨率的各频率分量。

(3) 若x′(t)的频率分量中确实有频差小于物理频率分辨率的分量,则这些分量在x(n)的频谱中会显示出“混叠效应”，从而在DFT谱分析中无法辨识。

用N点DFT进行序列频谱分析,给定采样频率fs，计算频率分辨率Δfbin，根据式(3.5.3)的定义可知，Δfbin是由DFT的点数N限定的。

因此，在给定采样频率fs的条件下，若长为TL(或长为L点)的序列x(n)的频谱物理频率分辨率Δf=1TL=fsL足够高，就可以清楚地显示原x′(t)的最小频率变化Δfmin情况，则DFT点数N增加，会得到更多离散频率点ωk=k2πN（k=0,1,…，N-1）处的频谱值，即DFT各谱线对应角频率间隔ΩI=2πNT=2πfsN减小，使计算频率分辨率Δfbin提高，可分辨出序列x(n)中具有更小频率变化的有关频率分量，从而提高对原连续时间信号x′(t)的频率成分的估计精度。

若截断信号的物理频率分辨率Δf差，即x(t)的总长度(TL)短，则对应x(n)总点数L少，x(n)的频谱已不能反映原x′(t)的实际最小频率变化Δfmin，即物理频率分辨率Δf >Δfmin； 则通过增加DFT的点数N，虽然提高了计算频率分辨率，Δfbin很小，但也无法识别实际x′(t)中的具有最小频率变化Δfmin的相关频率分量。






3.5.3〓DFT对连续时间信号进行谱分析——实例

微课视频


3.5.3离散傅里叶变换分析连续时间信号x(t)实例

例3.5.1设模拟信号x(t)有3个频率分量f1=2kHz,f2=2.5kHz,f3=3kHz，具体如下：


x(t)=cos(2πf1t)+cos(2πf2t) +cos(2πf3t),t的单位是ms


以采样频率fs=10kHz对其采样。取采样序列的采样点数L=10,20两种情况，对每种情况分别进行N=32,64点的DFT。问所得频谱能否区分x(t)中的3个频率分量?

解: 例3.5.1的MATLAB程序参考实现见3.8.3节。

给定模拟信号的最高频率成分为3kHz，故选择10kHz的采样频率对x(t)进行时域采样，满足奈奎斯特定理，采样不会引起信号失真。运行3.8.3节中的相应程序，可绘出对连续时间信号x(t)以fs=10kHz进行时域采样后，根据两种持续时间长度截断所得到的序列x1(n)、x2(n)的时域波形，如图3.5.2所示。




图3.5.2连续时间信号以fs采样后，分别取L=10，20个样值的时域波形


根据频率分辨率的定义，采样频率fs=10kHz时，取L个样点的序列x(n)对应的物理频率分辨率Δf=fsL; N点DFT对应的计算频率分辨率Δfbin=fsN。下面就采样点数L、DFT的点数N取不同值时，对模拟信号x(t)进行谱分析的结果（特别是频谱分析误差）进行研究。

依题意可知，模拟信号x(t)的最低频率分量为2kHz，而最小频率间隔Δfmin=0.5kHz； 综合考虑，为了能正确检测x(t)的各频率成分，选择的频谱分析系统的频率分辨率值最好小于0.5kHz。

（1） 采样点数L=10时，截断序列x1(n)的物理频率分辨率为


Δf=fsL=1010kHz=1kHz


数值大于实际模拟信号x(t)的最小频率间隔0.5kHz，因此截断对应的10点样值序列x1(n)的频率分辨率较差。x1(n)的频谱(DTFT)相对于原模拟信号x(t)的频谱会出现混叠现象，如图3.5.3所示，平坦的频谱包络不能区分模拟信号x(t)中的2kHz、2.5kHz和3kHz的频率分量。



图3.5.3长分别为L=10,20的截断序列的DTFT幅度谱


这时若进行N点DFT，频谱分析对应的计算频率分辨率分别为 


N=32时，Δfbin=fsN=1032kHz=0.3125kHz
N=64时，Δfbin=fsN=1064kHz=0.15625kHz


数值都小于原模拟信号x(t)的最小频率间隔，计算频率分辨率够好。但两种点数的DFT结果并不能区分实际信号中的3个频率分量，如图3.5.4(a)所示。说明截断信号的物理频率分辨率差的条件下，DFT的计算频率分辨率精度如何已无关紧要。



图3.5.4长为L=10及
N=32,64
的截断序列的DFT幅度谱


（2） 采样点数L=20时，截断序列x2(n)的物理分辨率为


Δf=fsL =1020kHz=0.5kHz



数值刚好等于实际模拟信号的最小频率间隔0.5kHz。因此截断对应的20点样值序列x2(n)的频谱(DTFT)恰能区分x(t)的3个频率分量，如图3.5.3中0~fs2范围内的3个谱峰所示。但有两个谱峰偏离实际的2kHz、3kHz频率点(虚线指示的位置)。在这种情况下，选用大点数N的DFT，对应的计算频率分辨率将提高，可以在一定程度上减少频谱估计误差。如图3.5.5所示，其中符号×处表示实际信号频率点，可见当N=64时，存在着更近似于实际频点的DFT谱线。




图3.5.5长为L=20及N=32,64的截断序列的DFT幅度谱


(3) 误差分析。

在长L点的截断序列x(n)的DFTF可以反映原连续时间信号x(t)的各频率分量条件下，如图3.5.3(b)所示情况，图3.5.5所示的N点DFT的离散谱线ki（i=0,1,…,N-1）对应的模拟域频率值为fki：


fki=kifsN(3.5.4)


N点DFT X(k)的序列号ki为整数,对应的频率值fki与实际的频率分量fi可能有误差。如果对实际频率分量fi,存在整数kfi,满足：


fi=kfifsN,
kfi= fiNfs(3.5.5)


则整数kfi必为DFT的某一谱线序列号，可精确指示实际频率分量fi。例3.5.1中，fs=10kHz，f2=2.5kHz，在N=32点的DFT中，由式(3.5.5)计算序号： 


kf2=2.5kHz×3210kHz=8


结果是整数，表明32点DFT序号为8的谱线精确地指示出了信号中的2.5kHz频率分量。进一步，对N/2=32/2=16中心对称的X(k)序号N-kf2=32-8=24则给出了频率分量f2的镜像频率点，图3.5.5中显示该频率点处，|X(24)|是局部峰值。

反之，对DFT的整数数字频率点序号ki，根据式(3.5.4)计算fki,如果仅有fki≈实际信号频率分量fi，这时整数ki表示的频率分量是模拟域实际频率分量fi的近似值。例3.5.1中，当N=64，fs=10kHz时，由式(3.5.5)计算连续时间信号中f1=2kHz的频率分量对应的
“序号”：


kf1=2kHz×6410kHz=12.8


结果不是整数，则在64点DFT X(k)中与kf2=12.8近似的整数序号是k=12、13，所以DFT数字频率点序号k1=12或13处的谱线可用来近似实际频率分量f1=2kHz的谱线。

根据DFT对N/2的中心对称性： N-k1=64-13=51点或N-k1=64-12=52点的谱线也近似对应着实际谱分量f1的镜像。这种近似误差的原因是物理频率分辨率或计算频率分辨率不够好（主要取决于两者中指标更差的那个频率分辨率）。例3.5.1中，20点序列的物理频率分辨率仅为0.5kHz，远劣于32点DFT对应的计算频率分辨率0.3125kHz，因此是32点DFT频谱分析近似误差的主要原因，可通过增加数据长度L，减小此偏差。

例3.5.1中实际频率fi（i=1，2，3）与N=32,64点DFT序列中对应的数字频率序号ki如表3.5.1所示，表中kfi同时给出了实际频率的镜像频率点对应“序号”，但ki仅给出了近似实际频率谱线的序号，其中粗体DFT序号是最接近实际kfi的序号。


表3.5.1由式(3.5.5)计算出kfi，近似表示fi的DFT数字频率序号ki




对应的序号kfif1 f2f3
N=326.4,25.68,249.6,22.4

N=6412.8,51.216,4819.2,44.8取整的序号ki=［kfi］

DFT离散频点
(省略镜像)序号
f1f3
6或79或10
12或1319或20



由图3.5.5和表3.5.1可见，除2.5kHz频率分量外，f1=2kHz,f3=3kHz两个频率，对应非整数的kfi。说明序列号为整数的DFT的峰值谱线仅近似表示了x(t)信号的这两个频率分量。若截断序列的物理频率分辨率够好，定义最接近实际的归一化离散频率为［kfi］N，［·］表示取整数，它与实际归一化频率fifs之差，可定义为误差: 


E=fifs-［kfi］N(3.5.6)


该误差随DFT的长度N的增加而减小，这是因为N增加，计算频率分辨率就会提高。例如fs=10kHz,f1=2kHz,f2=2.5kHz,f3=3kHz，实际的归一化频率和最接近的32点、64点DFT的归一化离散频率间的误差见表3.5.2，与理论分析一致。


表3.5.2不同点数DFT峰值谱线表示实际x(t)信号的频率分量的误差



点数f1/fsf2/fsf3/fs取整f1/fsf3/fs误差E


N=32
N=640.20.250.3［kfi］N
0.18750.3125±0.0125
0.2031250.296875±0.0031

备注DFT近似离散频率对实际频率的误差,随DFT的长度增加而减小







3.6〓DFT对离散时间信号的频谱分析

微课视频


3.6用离散傅里叶变换对离散时间信号进行谱分析


前面已经介绍，离散时间信号序列x(n)的谱是其DTFT函数X(ejω)，也是序列的z变换X(z)在z平面单位圆上的情况，即z=ejω时X(z)的特例。根据DFT的定义，DFT只适合对有限长序列和周期序列进行谱估计，下面就对这两种情况分别讨论。

1. 有限长序列x(n)

考虑到有限长序列x(n)的DFT——X(k)就是其z变换在z平面的单位圆，数字频率区间［0,2π］上的等间隔有限采样。所以，有限长序列的DTFT函数X(ejω)可直接用DFT估计，给出有限数字频率点的信号谱特性。

2. 周期为N的序列x～(n)

因周期序列x～(n)不满足绝对可和条件，其严格意义的DTFT不存在。但借助冲击函数δ(ω)，其频谱X(ejω)可记为


X(ejω)=2πN∑∞k=-∞X～(k)δω-2πNk


其中，


X～(k)=∑N-1n=0x～(n)e-j2πNkn，k∈(-∞,+∞)的整数


N为离散时间信号(周期序列)的周期。显然周期序列的谱X(ejω)是以N为周期的离散谱，谱线相对大小与周期序列离散傅里叶级数(DFS)的系数X～(k)成正比。用DFT分析周期序列x～(n)的谱，要先对其截断。

(1) 如果截取的是周期序列x～(n)的主值序列x(n)。

根据DFT隐含的周期性： 由截取的主值序列计算出的N点DFT序列X(k)正是周期序列的DFS X～(k)的主值序列，故X(k)可以表示周期序列x～(n)的谱结构。

(2) 如果截取的是周期序列的m个周期，序列长mN。

由序列的m个周期样值计算出的M=mN点的DFT XM(k)与由主值序列算出的N点DFT X(k)之间可以证明有如下关系： 


XM（k）=mXkm，km=整数
0，km≠整数，k=0，1，…，mN-1(3.6.1)


可见，只要截取x～(n)的整数个周期进行DFT，XM(k)也能表示x～(n)的频谱结构，达到谱分析的目的。需要注意，式(3.6.1)表示的m个周期的有限长序列的mN点DFT XM(k)和以前介绍的后补零成为长mN点序列，再进行mN点的DFT得到的谱G(k)是不同的。主要差别如下： 

取周期序列的m个周期参与离散傅里叶变换运算，其信号总能量是一个周期长的序列总能量的m倍。因此，各谐波幅度相对于单周期的X(r)是成m倍数的； 如果只补零，信号总能量未变，相应各谐波幅度自然不变。

(3) 如果时域带限x～(n)序列周期未知。

这时要对周期序列的谱进行DFT估计，可先取较短的M点一段xM(n)，做M点的DFT：


XM(k)=DFT［xM(n)］,k=0,1,…，M-1


再取2M长的一段时域序列x2M(n)，计算2M点的DFT：


X2M(k)=DFT［x2M(n)］,k=0,1,…，2M-1


比较XM(k)、X2M(k)，如果两者的主谱线差别满足谱分析频率误差要求，则可用XM(k)或X2M(k)来近似x～(n)的频谱。否则，继续将截取长度加倍，直到连续两次的主谱线对应的频率
误差满足误差要求。设最后截取的长度为rM，则XrM(k)表示数字频率ω=2πkrM（k=0,1，…,rM-1）处的离散时间周期信号谱线强度。






3.7〓DFT频谱分析的典型现象

微课视频


3.7离散傅里叶变换应用中的问题与参数选择

在应用DFT解决实际问题时常遇到下列几个问题： ①混叠现象； ②栏栅效应； ③频率泄漏。这些问题与应用中信号和DFT的参数选择有关。下面分别讨论这些问题。

3.7.1混叠现象
1. 混叠产生的机理

数字信号处理中的混叠可以指时域的混叠，也可以指频域的混叠。混叠主要是由采样引起的。

因为对连续时间信号x(t)以采样频率fs进行采样，所得采样序列x(n)的频谱是原连续时间信号频谱以fs为周期的周期延拓。所以，处理连续时间信号前，常令其通过截止频率为fh的前置低通滤波器，保证滤波器输出信号的双边谱带宽为2fh，即得到最高频率≤ fh的模拟信号xa(t)。然后，根据奈奎斯特采样定理，令:


fs≥2fh


即采样时间间隔T=1fs≤12fh，对xa(t)进行时域采样。工程应用中，常取fs≈(3~4)fh，保证不会因为对xa(t)的时域采样，导致信号xa(t)的频谱在频域以周期fs延拓，得到序列的谱时产生频谱混叠，形成xa(t)的频谱失真，失去进一步进行数字信号处理的依据。因为频谱混叠现象得到的是xa(t)失真的频谱，将不能从失真的频谱中恢复出原来的信号。

另一方面，在满足奈奎斯特采样定理的条件下，时域采样序列x(n)的频谱(DTFT)主周期即为实际连续时间信号x(t)谱结构，用N点DFT处理x(n)，结果X(k)等效于在频域对信号x(n)的频谱主周期，频率范围［0,fs］上的频谱函数进行N点等间隔采样。如果通过IDFT把X(k)变换回时域序列时，所得x′(n)可看成是原序列x(n)以N为周期的周期延拓序列的主值序列，如图3.7.1所示。

显然，如果原序列x(n)是L点有限长，只要N≥L，就能保证由DFT X(k)恢复原采样序列； 这个特点可以换个角度来应用： 如果我们希望处理信号的物理频率分辨率为Δf，则最短的时域记录长度TL=1Δf。当采样频率fs已知，即时域采样间隔T=1fs给定时，应选择DFT的点数N≥TLT。

如果原序列x(n)是无限长的，则DFT近似处理，无论取多大点数N，IDFT对应的有限长序列与原无限长序列相比，总有时域混叠效应，特别是以N为周期延拓的时域序列的主值序列两端时域混叠效应会更严重。

2. 避免混叠现象、保证处理要求的主要措施

对连续时间信号进行前置低通滤波，限定信号的最高频率为fh。在满足奈奎斯特采样定理时，如果期望的物理频率分辨率参数Δf给定，对采样序列x(n)截采样点数N必须满足: 


N=TLT≥2fhΔf(3.7.1)


其中，T为采样时间间隔。T不变时，可改变记录长度TL，即增加截取的样点数N以满足式(3.7.1)的不等号。

如果采样点数N给定(工程上一般对应着DFT的固定点数N)，则式(3.7.1)涉及信号最高频率fh,物理频率分辨率参数Δf的选取在考虑频率分辨率指标时是矛盾的。因为当fh变大时，奈奎斯特采样定理要求采样时间间隔T减小； 式(3.7.1)等号表明信号总时长TL必须相应减小，以保证N不变； 而物理频率分辨率数值Δf=1TL将变大，使分辨率指标变差。

例3.7.1对实信号进行谱分析，要求谱分辨率Δf≤10Hz，信号最高频率fc=25kHz，试确定最小记录时间TLmin，最大的采样时间间隔Tmax，最少的采样点数Nmin。如果fc不变，要求谱分辨率提高一倍，最少的采样点N和最小的记录时间又是多少? 

解： 已知谱分辨率Δfmax=10Hz，由式(3.5.1)知TL≥1Δf。因此最小记录时间: 


TLmin=110Hz=0.1s


因为信号最高频率fc=2.5kHz，奈奎斯特采样定理要求采样频率fs≥2fc，所以最大的采样时间间隔：


Tmax=12fc=12×2.5kHz=0.2ms


对应地，最少的采样点数: 


Nmin=TLminTmax=0.1s0.2ms=500(点)


第4章中将介绍DFT的点数为2的整数次幂时的快速实现，可取Nmin=512点。

如果信号最高频率fc不变，即最大的采样时间间隔Tmax=0.2ms不变； 要求谱分辨率提高一倍，即谱分辨率数值减小为Δfmax=5Hz，则最小的记录时间和最少的采样点N分别为 


TLmin=15Hz=0.2s
Nmin=TLminTmax=0.2s0.2ms=1000点


工程上DFT的点数为2的整数次幂时，能够快速实现，可取Nmin=1024点。

3.7.2栏栅效应
1. 栏栅效应产生的机理


如图3.7.1所示，非周期信号xa(t)的频谱幅度函数|Xa(jf)|是连续的。把xa(t)的N点采样x(n)进行DFT得到的频谱幅度序列|Xa(k)|只能是连续频谱Xa(jf)的幅度函数的有限离散频点采样。



图3.7.1频域采样X(k)的IDFT，对应于时域序列x(n)以N为周期延拓的主值序列


2. 栏栅效应

N点DFT结果好像是在栏栅的一边，通过N个有限离散频点的“缝隙”观察另一边的频谱Xa(jf)。故这种现象被形象地称为“栏栅效应”，只有在离散频点的“缝隙”处才看得见。两根谱线(“缝隙”)之间如有重要的频谱分量,则将被错过而检测不出。

3. 克服的措施

为了把被“栏栅”挡住的频谱分量检测出来，可在原记录序列后面补零,增加DFT的长度，即增加频域Xa(jf)上的采样点数N,改变离散谱线(“栏栅上缝隙”)的分布，就可能检测出原来看不到的频谱分量。

但这样不能提高物理频率分辨率，即原先不能分开的频率分量，仅后补零仍然不能分开，例如图3.5.4中，三个频谱分量混叠成的平坦的频谱带，没有因为增加了DFT的点数N而分开成三个。但在图3.5.5中，序列长L=20点的情况，物理频率分辨率能保证分开主要频率分量，增加频域采样点为N=64时,才会比N=32看到的离散频谱X(k)更接近图3.5.3(b)所示的实际情况。

3.7.3频率泄漏
1. 频率泄漏现象

频率泄漏现象指连续时间信号x(t)经系统处理后，所得的频谱在原来没有频谱的频率区间出现了频谱。

2. 频率泄漏产生机理

实际处理信号序列x(n)时，根据需要常要对其截短。最简单的截短就是将该序列限定为有限的N点，令该序列N点以外的序列值为零，等效于长的数据序列乘以长为N时间窗函数ω(n)，如矩形窗RN(n)。

例3.7.2设矩形窗RN(n)=1,n=0，1，…,N-1。取N=11，画出其时域和频谱图。

解： 矩形窗的频谱为


RN(ejω)=DTFT［RN(n)］= e-jωN-12sin(ωN/2)sin(ω/2)=sin(ωN/2)sin(ω/2)ejφ(ω)


当N=11时，其时域波形和幅度频谱如图3.7.2所示。图3.7.2(b)中，|ω|<2π/N的部分称为窗谱主瓣，其余部分称为窗谱旁瓣。



图3.7.2矩形窗R11(n)的时域、频域波形



根据频域卷积定理，时域中x(n)与ω(n)相乘，则频域中X(ejω)与ω(n)的谱函数W(ejω)进行卷积，频域卷积结果分布在更宽的频率范围上。这表明x(n)截短后的频谱不同于它原未截短的频谱。产生了原信号中没有的频率分量，从频谱上看,截短序列的谱扩散到原来没有频谱的频率区间——这就称为“频率泄漏”。

这种频率泄漏，有时还会产生不同频率分量间干扰，即频率“混叠现象”。

例3.7.3序列x(n)=cos(ω0n),ω0=π/4。其频谱为线谱X(ejω):


X(ejω)=π∑∞k=-∞δω-π4-2πk+δω+π4-2πk


x(n)的幅度谱主值波形如图3.7.3(a)所示。

将x(n)用RN(n)截短后，得y(n)=x(n)RN(n)，其幅度谱主值波形如图3.7.3(b)所示，显然原来仅在π/4处的频谱“泄漏”到了整个［0,2π］的频率范围。例3.7.3截短余弦序列y(n)的谱可以形象地认为是由矩形序列的窗谱主瓣、旁瓣引起的。



图3.7.3cos(nπ/4)加矩形窗前后的频谱


3. 克服频率泄漏的措施

频率泄漏现象是由对时域序列x(n)加窗造成的，如图3.7.2给出矩形窗的时域与频域波形。由图3.7.2可见，加大时域窗口宽度N，频域窗谱主瓣和窗谱旁瓣将变窄，可使频率泄漏减少。但无限增加时域窗口宽度，等于对原序列不截短，不满足截短要求。因此，为尽量减少频率泄漏，应选用窗谱旁瓣小、主瓣窄，即频率“泄漏”小的窗函数。这个问题将在第7章有限长脉冲响应滤波器设计中详细讨论。

3.8MATLAB实现


本节给出的MATLAB实现参考程序，均在MATLAB 7.4.0环境下调试通过。

3.8.1离散傅里叶变换的MATLAB实现

例3.2.1的参考MATLAB实现如下： 

%the effect of zero padding to the DFT (II,padding in front)

x=［1 1 1 1 0 0 0 0］; y=［0 0 0 1 1 1 1 0 0 0 0］;

D=3;   %number of zeros padding in front

n1=16; f1=fft(x,n1);  f2=fft(y,n1);

delta=-2*pi*D*［0:n1-1］./n1;

subplot(2,2,1);stem(x);axis(［0 16 0 max(x)］);text(8,0.8,'series of x(n)');

ylabel('x(n)');xlabel('n');

subplot(2,2,2);stem(［0:n1-1］./n1,abs(f1));axis(［0 1 0 max(abs(f1))］);

text(.2,3.7,'N=16 DFT of x(n)');ylabel('|X(k)|');xlabel('unit in 2*pi');

subplot(2,2,3);stem(y);axis(［0 16 0 max(y)］);text(8,0.8,'series of xD(n)');

ylabel('xD(n)'); xlabel('n'); 

subplot(2,2,4);stem(［0:n1-1］./n1,abs(f2));axis(［0 1 0 max(abs(f2))］);

text(.2,3.7,'N=16 DFT of xD(n)');ylabel('|XD(k)|');xlabel('unit in 2*pi');

figure

subplot(3,1,1);line(［0:n1-1］./n1, zeros(n1));hold on

stem(［0:n1-1］./n1,angle(f1));ylabel('Phase X(k) in rad');

subplot(3,1,2);line(［0:n1-1］./n1, zeros(n1));hold on

stem(［0:n1-1］./n1,angle(f2)); ylabel('Phase  XD(k)in rad');

subplot(3,1,3); line(［0:n1-1］./n1, zeros(n1));hold on; stem(［0:n1-1］./n1,delta); 

ylabel('shifted Phase in rad');xlabel('unit in 2*pi'); gtext('shifted Phase in XD(k) with X(k)')

例3.2.2的参考MATLAB实现如下： 

%the effect of zero padding to the DFT

x=［1 1 1 1 0 0 0 0］; n1=8; n2=16; n3=11;

f1=fft(x,n1); f2=fft(x,n2); f3=fft(x,n3);

stem(x); axis(［0 16 0 max(x)］); text(10,0.8,'series of x(n)'); ylabel('x(n)'); xlabel('n');

figure;

subplot(1,3,1); stem(［0:n1-1］./n1,abs(f1)); axis(［0 1 0 max(abs(f1))］); 

text(.2,3.7,'N=8 DFT of x(n)'); ylabel('|X(k)|'); xlabel('unit in 2*pi');

subplot(1,3,2); stem(［0:n2-1］./n2,abs(f2)); axis(［0 1 0 max(abs(f2))］); 

text(.2,3.7,'N=16 DFT of x(n)');

subplot(1,3,3); stem(［0:n3-1］./n3,abs(f3)); axis(［0 1 0 max(abs(f3))］);

text(.2,3.7,'N=11 DFT of x(n)');

figure;

subplot(3,1,1); stem(［0:n1-1］./n1,abs(f1)); axis(［0 1 0 max(abs(f1))］);

text(.4,3.2,'N=8 DFT of x(n)'); ylabel('|X(k)|');

subplot(3,1,2); stem(［0:n2-1］./n2,abs(f2)); axis(［0 1 0 max(abs(f2))］);

text(.4,3.2,'N=16 DFT of x(n)');

subplot(3,1,3); stem(［0:n3-1］./n3,abs(f3)); axis(［0 1 0 max(abs(f3))］);

text(.4,3.2,'N=11 DFT of x(n)'); xlabel('unit in 2*pi');

3.8.2用离散傅里叶变换计算线性卷积和圆周卷积的MATLAB
实现

例3.3.2的参考MATLAB实现如下： 

% comparison of convolution and circular convolution

x1=［1,2,2,1］;x2=［1,-1,1,-1］; x3=conv(x1,x2); x4=circonvf(x1,x2,7);

n=6; x5=circonvf(x1,x2,n);

subplot(2,1,1); stem(［0:length(x1)-1］,x1);hold on; plot(［0 ,max(n,7)］,［0 0］,'k');hold off

text(4,0.5*max(x1),'x1(n)'); xlabel('n');ylabel('x1');

subplot(2,1,2); stem(［0:length(x2)-1］,x2);hold on; plot(［0 ,max(n,7)］,［0 0］,'k');hold off

text(4,0.5*max(x2),'x2(n)'); xlabel('n');ylabel('x2');

figure

subplot(3,1,1); stem(［0:length(x3)-1］,x3);hold on; plot(［0 ,max(n,7)］,［0 0］,'k');hold off

text(4,0.5*max(x3),'Linear Convolution'); xlabel('n');ylabel('x3')

subplot(3,1,2); stem(［0:length(x4)-1］,x4);hold on; plot(［0 ,max(n,7)］,［0 0］,'k');hold off

text(4,0.5*max(x4),'7-Point Circular Convolutioné '); xlabel('n');ylabel('x4')

subplot(3,1,3); stem(［0:length(x5)-1］,x5);hold on; plot(［0 ,max(n,7)］,［0 0］,'k');hold off

string=［num2str(n),'-Point Circular Convolution'］;

text(4,0.5*max(x5),string); xlabel('n');ylabel('x5')

3.8.3连续时间信号谱分析的MATLAB实现

例3.5.1的参考MATLAB实现如下： 

% frequency resolution, frequency interferences

clear all;

freq1=2;freq2=2.5;freq3=3.;  fs=10;%frequency unit in kHz;

omega1=2*pi*freq1; omega2=2*pi*freq2; omega3=2*pi*freq3;

N=1024;% number of DFT used for simulate DFTF

L1=10;L2=20;% the number of sample poits

N1=32; N2=64;% the numbers of DFT

dt=1./fs; % x(n) time interval, unit is in ms

ddt=dt./100;% time interval decrease to a hundredth, unit is in ms

f=［freq1 freq2 freq3］; z=［0 0 0］;zz=［1 1 1］;

t1=0:ddt:(L1-1)*dt; tt1=0:ddt:(L2-1)*dt;

xt1=cos(omega1*t1)+cos(omega2*t1)+cos(omega3*t1);%short analog signal

xt2=cos(omega1*tt1)+cos(omega2*tt1)+cos(omega3*tt1);%longer analog signal

 

k1=0:dt:(L1-1)*dt; kk1=0:dt:(L2-1)*dt;

x1=cos(omega1*k1)+cos(omega2*k1)+cos(omega3*k1);%short sample sequence

x2=cos(omega1*kk1)+cos(omega2*kk1)+cos(omega3*kk1); %longer sample sequence

ydtft1=fft(x1,N);  ydtft2=fft(x2,N); fdtft=fs*linspace(0,1,N);

 

y1=fft(x1,N1);i1=0:N1-1;  y12=fft(x1,N2);i2=0:N2-1;

y2=fft(x2,N1);  y22=fft(x2,N2);

 

subplot(2,1,1);% time domain signals and its sequences

plot(10*t1,xt1,'--');hold on; stem(0:L1-1,x1,'b');hold off

title(' sampled signal to x1(t)'); ylabel('x1(n)'); axis(［0 L2 min(x1) max(x1)］);

subplot(2,1,2)

plot(10*tt1,xt2,'--');hold on; stem(0:L2-1,x2,'b');hold off

title(' sampled signal to x2(t)'); xlabel('n');ylabel('x2(n)'); axis(［0 L2 min(x2) max(x2)］);

 

figure % frequency domain,DTFTs to L=10,20 sequences 

subplot(2,1,1); plot(fdtft,abs(ydtft1));hold on;

stem(f, zz*max(abs(ydtft1)) , '--xr');hold off; ylabel('DTFT Spectrum |X_1(e^j^＼omega)|')

axis(［min(fdtft) max(fdtft) min(abs(ydtft1)) max(abs(ydtft1))］);

title('DFTF to signal x_1(n), L=10'); 

subplot(2,1,2); plot(fdtft,abs(ydtft2));hold on;

stem(f, zz*max(abs(ydtft2)) , '--xr');hold off; ylabel('DTFT Spectrum |X_2(e^j^＼omega)|')

axis(［min(fdtft) max(fdtft)  min(abs(ydtft2)) max(abs(ydtft2))］);

title('DFTF to signal x_2(n), L=20'); xlabel('frequency f(kHz)')

 

figure %DFTs to L=10 sequence

subplot(2,1,1); stem(i1,abs(y1));hold on; stem(f*N1/fs,z , 'xr');hold off

ylabel('Magnitude Spectrum |X(k)|'); axis(［0 max(i1) min(abs(y1)) max(abs(y1))］);

title('signal x1, L=10,N=32')

subplot(2,1,2)

stem(i2,abs(y12));hold on; stem(f*N2/fs,z, 'xr');hold off; xlabel('k'); 

ylabel('Magnitude Spectrum |X(k)|'); axis(［0 max(i2) min(abs(y12)) max(abs(y12))］); 

title(' signal x1 , L=10, N=64')

 

figure

subplot(2,1,1); stem(i1,abs(y2));hold on;stem(f*N1/fs,z , 'xr');hold off

ylabel('Magnitude Spectrum |X(k)|');axis(［0 max(i1) min(abs(y2)) max(abs(y2))］);

title('signal x2, L=20,N=32')

subplot(2,1,2);

stem(i2,abs(y22));hold on; stem(f*N2/fs,z, 'xr');hold off; xlabel('k')

ylabel('Magnitude Spectrum |X(k)|'); 

axis(［0 max(i2) min(abs(y22)) max(abs(y22))］); title('signal x2, L=20,N=64')

3.9习题


3.1计算下列信号的傅里叶变换在［0,2π］上的N点均匀采样值。

（1） acos(ω0n)RN(n)；

（2） anRN(n)；

（3） n2RN(n)；

（4） a|n|sin2πNk0nu(n),|a|<1；

（5） 12n(u(n+3)-u(n-2))；

（6） n(u(n+N)-u(n-N-1))；

（7） cos18πn7+sin(2n)；

（8） ∑∞k=014nδ(n-3k)；

（9） x(n)=cosnπ3,-1≤n≤40,其他n。

3.2根据下列离散时间信号傅里叶变换的N点离散采样值，确定各相应的离散时间域信号。

（1） X(k)=0,0≤2πNk<ω
1,ω≤2πNk≤π；

（2） X(k)=1-2e-j32πNk+4ej22πNk+3e-j62πNk；

（3） X(k)=∑+∞m=-∞(-1)mδ2πNk-π2m；

（4） X(k)=cos22πkN；

（5） X(k)=cosπkN+jsin2πkN,-π≤2πNk≤π；

（6） X(k)=e-j2πNk1+16e-j2πNk-16e-j4πNk。

3.3已知下列N点DFT X(k)，求x(n)=IDFT［X(k)］。

（1） X(k)=N2ejθ,k=mN2e-jθ,k=N-m0，其他k

（2） X(k)=-jN2ejθ,k=mjN2e-jθ,k=N-m0，其他k

其中，m为正整数，0<m<N2，N为变换区间长度。

3.4序列x(n)=δ(n)+2δ(n-2)+δ(n-3)： 

(1) 求x(n)的4点DFT；

(2) 若y(n)是x(n)与它本身的4点循环卷积，求y(n)及其4点DFT Y(k)；

（3） 若h(n)=δ(n)+δ(n-1)+2δ(n-3)，求x(n)与h(n)的4点循环卷积。

3.5设X(k)表示N点序列x(n)的N点离散傅里叶变换，X(k)本身也是N点序列。如果计算X(k)的离散傅里叶变换得到一序列x1(n)，试用x(n)表示x1(n)。

3.6设有两个序列：


x(n)=x(n),0≤n≤50,其他n
y(n)=y(n),0≤n≤140,其他n


各作15点的DFT，然后将两个DFT相乘，再求乘积的IDFT，设所得的结果为f(n)，问f(n)的哪些点(用序号n表示)对应于x(n)*y(n)应该得到的点。

3.7已知两个有限长序列为


x(n)=n+1,0≤n≤30,4≤n≤6
y(n)=-1,0≤n≤41,5≤n≤6


试用作图表示x(n),y(n)以及x(n)⑦y(n)。

3.8已知x(n)是N点有限长序列，X(k)=DFT［x(n)］，现将x(n)的每两点之间补进r-1个零点值，得到一个rN点的有限长序列y(n)：


y(n)=xnr,n=ir,i=0,1,…,N-10,其他n


试求rN点的DFT［y(n)］=Y(k)与X(k)的关系。

3.9设有一谱分析用的信号处理器，采样点数必须为2的整数幂，假定没有采用任何特殊数据处理措施，要求频率分辨率小于或等于10Hz，如果采用的采样时间间隔为0.1ms，试确定： 

(1) 最小记录长度；

(2) 所允许处理的信号的最高频率；

(3) 在一个记录中的最少点数。

3.10给定长为N的圆周共轭对称序列xep(n),圆周共轭反对称序列xop(n)， 0≤n≤N-1。试证明： 当N为偶数时，xepN2为实数； xopN2为纯虚数。

3.11序列x(n)的共轭对称和共轭反对称分量分别为


xe(n)=12［x(n)+x*(-n)］,xo(n)=12［x(n)-x*(-n)］


长度为N的有限序列x(n)(0≤n≤N-1)的圆周共轭对称和圆周共轭反对称分量分别定义如下：


xep(n)=12［x((n))N+x*((-n))N］RN(n)
xop(n)=12［x((n))N-x*((-n))N］RN(n)


(1) 证明：


xep(n)=［xe(n)+xe(n-N)］RN(n)
xop(n)=［xo(n)+xo(n-N)］RN(n)


(2) 把x(n)看作是长度为N的序列，一般来说，不能从xep(n)恢复xe(n)，也不能从xop(n)恢复xo(n)。试证明若把x(n)看作是长度为N的序列，且n≥N2时，x(n)=0，则从xep(n)可恢复xe(n)，从xop(n)可恢复xo(n)。

3.12设X(k)表示N点序列x(n)的N点离散傅里叶变换：

(1) 证明如果当x(n)满足关系式x(n)=-x(N-1-n),则X(0)=0；

(2) 证明当N为偶数时，如果x(n)=x(N-1-n),则XN2=0。

3.13证明DFT对称定理，即假设:


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


证明DFT［X(n)］=Nx(n-k)。

3.14如果X(k)=DFT［x(n)］，证明DFT的初值定理：


x(0)=1N∑N-1k=0X(k)


3.15证明离散相关定理。

若
X(k)=X*1(k)X2(k)
则
x(n)=IDFT［X(k)］=∑N-1l=0x*1(l)x2（(l+n)）NRN(n)
3.16已知序列x(n)=anu(n)，0<a<1,对x(n)的z变换X(z)在z平面单位圆上等间隔采样N点，采样值为


X(k)=X(z)|z=ej2πNk，k=0,1,…,N-1


求有限长序列xN(n)=IDFT(X(k))。

3.17已知调幅信号的载波频率fc=1kHz，调幅信号的频率fm=100Hz，用FFT对其进行谱分析，试问： 

(1) 最小记录时间Tpmin。

(2) 最低采样频率fsmin。

(3) 最少采样点数Nmin。

(4) 在频带宽度不变的情况下，将频率分辨率提高一倍的采样点数N值。

3.18以20kHz的采样率对最高频率为10kHz的带限信号xa(t)采样，然后计算x(n)的N=1000个采样点的DFT，即


X(k)=∑N-1n=0x(n)e-j2πNnk,N=1000


（1） k=150对应的模拟频率是多少？k=800呢？

(2) 采样频率点之间的间隔是多少？(提示： 结果请注明单位)


思政小课堂


思政主题：  国之重器中的数字信号处理技术

思政方向：  培养专业认同感，民族自豪感，创新思维




数字信号处理技术可以带领我们探索星辰、大海，带领我们拥有自己的指路“北斗”和太空“天宫”。

1. 中国天眼——500m口径球面射电望远镜

（fivehundredmeter aperture spherical radio telescope，FAST）



“中国天眼”（FAST，图C.1）是目前全球最大单口径、最灵敏的射电望远镜。它的口径有500m，有近30个足球场大的接收面积。其综合性能相比于美国“阿雷西博”望远镜，提高了10倍。从提出构想到建设完成，我国三代科技工作者历经22年，先后攻克了十多项世界级难题。



图C.1中国天眼


FAST的工作原理是：  通过抛物面把采集到的电磁场信号反射到接收信号的馈源舱里，馈源舱内部部署了数字信号处理系统，对接收到的信号进行处理。

为什么FAST比当前其他射电天文望远镜对无线电信号的采集都要灵敏呢？正是因为我国科学家一个创新的设计理念：  通过组合反射面的运动与馈源舱的运动相互配合，FAST可以实现对来自太空平行光的点聚焦。无论反射面如何变化，通过馈源支撑塔悬挂的轻型索拖动机构的牵拉和并联机器人的二次精调，馈源舱都能随着反射面的变化需求，在规定时间准确无误地到达规定的位置，其误差小于10mm！

FAST项目是人类直接观测遥远星系行星，寻找类似太阳系或地球的宇宙环境，以及探测潜在智慧生命的重要设施，在航天工程及其他领域具有广泛用途。总而言之，FAST不仅极大地推动了我国在天文领域基础研究与科研攻关的步伐，还催生了一系列具有国际影响力的科学产出。

2. 中国“蛟龙号”载人潜水器

“蛟龙号”载人潜水器是我国第一艘自行设计、自主集成研制、拥有自主知识产权的深海载人潜水器（图C.2）。


图C.2“蛟龙号”载人潜水器


2012年6月27日，“蛟龙号”载人潜水器成功下潜到马里亚纳海沟7062m深度，创造了世界同类作业型潜水器的最大下潜深度纪录。



“蛟龙号”是全球首台应用高速数字水声通信技术的载人潜水器，它利用高速率水声通信技术，通过各类声学元器件传送声波信号给海面的控制船，使“蛟龙号”和控制船之间实现了语音、文字、图像的数字通信。“蛟龙号”在探索深海中发挥着重要作用。

3.  中国北斗卫星导航系统

北斗卫星导航系统（Beidou  navigation satellite system，BDS，以下简称北斗）是中国着眼于国家安全和经济社会发展需要，自主建设、独立运行的卫星导航系统。北斗不仅广泛进入交通、能源、农业、通信、气象、自然资源、生态环境、应急减灾等重点行业，而且全面服务于绿色出行、外卖送餐、健康养老、儿童关爱、医疗教育等人民生活的方方面面。2021年国内智能手机出货量中支持北斗的占总出货量的94.5%。北斗既造福中国人民，也造福世界各国人民。

我们一起来看看北斗建设中的“自主创新”密码：   “北斗一号”（图C.3（a））中我们提出了新颖的全数字信号处理技术，解决了当时众人难以解决的信号“快捕精跟”关键问题，超越国际同类技术的崭新方案；  “北斗二号”（图C.3（b））中，我们面临国外对原子钟技术的封锁，靠着不服输的骨气和艰苦奋斗的精神，自主研发了100%国产化的“中国心”原子钟；  在覆盖全球的“北斗三号”（图C.3（c））建设中，我们又一次创造性地提出了“星间链路”方案，前所未有地让太空兄弟在天上“手拉手”，大大提高了定位精度，创造了“北斗耀苍穹”的奇迹。而奇迹的背后是一群平均年龄30岁出头的年轻人持续26年的奋斗，他们用青春书写了“世界的北斗”。



图C.3北斗卫星导航系统





图C.3（续）


今天，当我们享受着这些国之重器为我们的社会和国防带来的安全、为我们的生活提供的便利时，我们应该向所有呕心沥血、无私奉献研发它们的科研工作者、建设者们说一声，谢谢！