第3章离散傅里叶变换 Z变换和离散时间傅里叶变换都是数字信号处理领域中重要的数学变换,但由于DTFT和单位圆上的Z变换在频域上都是连续的,因而适合于理论分析和数学推导,无法直接通过计算机进行数值运算。相对而言,离散的时间和频率适合于计算机处理。 对离散时间信号(序列)而言,若满足绝对可和条件,DTFT存在; 若不满足绝对可和条件,如周期序列,DTFT则不存在。此时,可以将周期序列看作周期连续时间信号经过采样得到的,与傅里叶级数(FS)展开类似,周期序列可以展开为离散傅里叶级数(DFS),获得离散的频谱。 对于有限长序列,一方面可视为周期序列的一个周期,通过DFS能够获得离散频谱; 另一方面对有限长序列的DTFT进行离散化,也同样可以得到离散频谱。这两种途径得到的离散频谱实际上就是离散傅里叶变换(DFT)。 DFT具有时域和频域离散化、有限长的特点,非常适合于计算机处理,在实际数字信号处理中发挥着重要作用。为了更好地理解和掌握DFT,下面首先讨论周期序列的离散傅里叶级数,然后介绍DFT的定义和性质、频域采样与内插、DFT典型应用以及MATLAB仿真等内容。 3.1周期序列的离散傅里叶级数 3.1.1周期序列 设x(n)是长度为N的有限长序列,0≤n≤N-1; x~(n)表示以N为周期,对x(n)进行周期延拓后形成的周期序列,那么它们的关系可表示为 x~(n)=∑∞m=-∞x(n+mN),m为整数(3.1.1) x(n)= x~(n)RN(n)(3.1.2) 上述关系如图3.1.1所示。x~(n)称为x(n)的周期延拓序列,从n=0到N-1的一个周期称为主值区间,x(n)作为x~(n)主值区间上的序列,称为主值序列。 图3.1.1有限长序列与周期序列的关系示意图 为了简便描述上述关系,定义x((n))N,表示x(n)以N为周期进行延拓后的周期延拓序列,即 x~(n)=x((n))N(3.1.3) 式中: ((n))N表示n对N求余。 若n=mN+n′,0≤n′≤N-1,m为整数,则((n))N=n′。 例如,N=6,x~(n)=x((n))6,则有 x~(7)=x((7))6=x((7-6))6=x(1) x~(6)=x((6))6=x((6-6))6=x(0) x~(-1)=x((-1))6=x((6-1))6=x(5) 由于周期序列x~(n)不满足绝对可和条件,根据DTFT存在条件可知,周期序列DTFT不存在。但从时域采样角度看,周期序列可视为周期连续时间信号的采样,则周期序列应该可以展开为类似傅里叶级数的形式,这就是周期序列的离散傅里叶级数。 3.1.2周期序列的DFS 下面首先从周期连续时间信号的傅里叶级数入手,得到周期序列的离散傅里叶级数展开式,然后推导展开式的系数表达式,最后定义离散傅里叶级数正变换(DFS)和离散傅里叶函数逆变换(IDFS)。 1. 周期序列的离散傅里叶级数展开形式 假设对周期为T的连续时间信号x(t)采样,采样间隔Ts=T/N,得到周期为N的周期序列x~(n),则 x~(n)=x(t)|t=nTs=x(nTs),-∞<n<∞(3.1.4) x(t)的傅里叶级数展开式为 x(t)=∑∞k=-∞akejkΩ0t (3.1.5) 式中,ak=1T∫T0x(t)e-jkΩ0tdt为系数,Ω0=2πT=2πNTs表示谐波的角频率间隔,k为谐波序号。将式(3.1.5)代入式(3.1.4),可得 x~(n)=∑∞k=-∞akejk2πNTsnTs=∑∞k=-∞akej2πNkn(3.1.6) 可以看出,周期序列可以展开为复指数序列的加权和形式,复指数序列为 ek(n)=ej2πNkn(3.1.7) 式中,k=1时表示基频序列ej2πNn,其他k值表示谐波序列,数字频率为ωk=2πNk。由于数字频率ω以2π为周期,有 ek+mN(n)=ej2πN(k+mN)n=ej2πNkn=ek(n),m为整数(3.1.8) 这说明,第k次谐波和第k+mN次谐波完全相同,谐波数目实际上只有N个。因此,式(3.1.6)的傅里叶级数展开项也只有N个,可以重新表述为 x~(n)=1N∑N-1k=0X~(k)ej2πNkn,-∞<n<∞(3.1.9) 式中,X~(k)为第k次谐波的系数,且X~(k)=N∑∞m=-∞ak+mN(k=0,1,…,N-1),引入常数1/N是为了X~(k)的计算方便和描述简洁。式(3.1.9)即为x~(n)的离散傅里叶级数展开式,表明周期序列可以表示为N个谐波的加权和形式,这些谐波成分表征了周期序列的频谱分布规律。 2. 离散傅里叶级数的系数X~(k)的表达式 为了从周期序列x~(n)的DFS中得到系数X~(k),式(3.1.9)两边同乘以e-j2πNrn(r=0,1,…,N-1),并对n=0到N-1的一个周期求和,可得 ∑N-1n=0x~(n)e-j2πNrn=∑N-1n=01N∑N-1k=0X~(k)ej2πNkne-j2πNrn =1N∑N-1k=0X~(k)∑N-1n=0ej2πN(k-r)n(3.1.10) 由于 ∑N-1n=0ej2πN(k-r)n=N,k=r 0,k≠r(3.1.11) 因此 ∑N-1n=0x~(n)e-j2πNrn=X~(r)(3.1.12) 利用变量k表示,即有 X~(k)=∑N-1n=0x~(n)e-j2πNkn(3.1.13) 式(3.1.13)是第k次谐波系数X~(k)的表达式, k=0,1,…,N-1。由于 X~(k+mN)=∑N-1n=0x~(n)e-j2πN(k+mN)n=∑N-1n=0x~(n)e-j2πNkn=X~(k),m为整数 (3.1.14) 因此,k取值可拓展为-∞<k<∞,且X~(k)具有周期性,周期为N。这表明周期序列及其离散傅里叶级数的周期是相同的。 3. 离散傅里叶级数变换对 基于式(3.1.9)和式(3.1.13),定义离散傅里叶级数变换对: 离散傅里叶级数正变换: X~(k)=DFS[x~(n)]=∑N-1n=0x~(n)e-j2πNkn,-∞<k<∞(3.1.15) 离散傅里叶级数逆变换: x~(n)=IDFS[X~(k)]=1N∑N-1k=0X~(k)ej2πNkn,-∞<n<∞(3.1.16) 其中: DFS[·]表示从时域到频域的正变换,IDFS[·]则表示从频域到时域的逆变换。DFS和IDFS具有相同的周期N,若取一个周期,如主值区间0≤n≤N-1和0≤k≤N-1,可以代表x~(n)和X~(k)的完整信息。 从上述讨论可以看出,周期序列和有限长序列有着紧密的联系,若只考虑周期序列主值区间内的DFS,就是有限长序列的离散傅里叶变换,这将在3.2节中进行介绍。 例3.1.1设x(n)=R4(n),将x(n)以周期N=8进行周期延拓后得到周期序列x~(n),如图3.1.2(a)所示, 求x~(n)的DFS。 解: 根据DFS计算式(3.1.15),可得 X~(k)=∑7n=0x~(n)e-j2π8kn=∑3n=0e-jπ4kn=1-e-jπk1-e-jπ4k =e-jπ2k(ejπ2k-e-jπ2k)e-jπ8k(ejπ8k-e-jπ8k)=e-j38πksinπ2ksinπ8k X~(k)的幅度|X~(k)|如图3.1.2(b)所示。 图3.1.2例3.1.1中周期序列及其离散傅里叶级数的幅频特性 从图3.1.2中可以看出: |X~(k)|的周期也为8,每个周期内有8根谱线; 对比例2.4.5中R4(n)的DTFT幅频特性图2.2.1,|X~(k)|可以视为|X(ejω)|在频域上以2π/N(N=8)进行等间隔抽样后得到的,|X(ejω)|在[0,2π)内主瓣和旁瓣各对应2根谱线。显然,若N=16,则|X~(k)|每个周期内将会有16根谱线,|X(ejω)|主瓣和旁瓣各对应4根谱线。 3.2离散傅里叶变换及性质 3.2.1DFT的定义 设序列x(n)长度为M,n=0,1,…,M-1,并以N为周期进行延拓得到周期序列x~(n)。采用DFS定义式(3.1.15)和IDFS定义式(3.1.16),将主值区间上的离散傅里叶级数变换对定义为离散傅里叶变换对。即x(n)的N点离散傅里叶变换定义为 X(k)=DFT[x(n)]=∑N-1n=0x(n)WknN,k=0,1,…,N-1(3.2.1) X(k)的离散傅里叶逆变换(IDFT)为 x(n)=IDFT[X(k)]=1N∑N-1k=0X(k)W-knN,n=0,1,…,N-1(3.2.2) 式中,WN=e-j2πN,N表示DFT变换区间的长度,M≤N。当M<N时,x(n)补0进行运算。 这里从DFS、IDFS出发,分别定义了DFT和IDFT,两者构成一对离散傅里叶变换。下面将从DFT和IDFT定义本身出发,通过推导来证明两者之间的变换关系。 将DFT定义式(3.2.1)代入IDFT定义式(3.2.2),可得 IDFT[X(k)]=1N∑N-1k=0∑N-1m=0x(m)WmkNW-knN =∑N-1m=0x(m)1N∑N-1k=0Wk(m-n)N(3.2.3) 由于 ∑N-1k=0Wk(m-n)N=N,n=m+rN 0,n≠m+rN,r为整数 在限定变换区间n=0,1,…,N-1的情况下,有r=0,因此,式(3.2.3)表示为 IDFT[X(k)]=x(n),n=0,1,…,N-1 上述推导过程表明,DFT和IDFT存在一一对应的时频域映射关系。 例3.2.1设有限长序列x(n)=R4(n),求x(n)的DTFT以及N为4点、8点、16点DFT。 解: (1) x(n)的离散时间傅里叶变换为 X(ejω)=∑∞n=-∞R4(n)e-jωn =∑3n=0e-jωn=1-e-j4ω1-e-jω =e-j2ω(ej2ω-e-j2ω)e-jω/2(ejω/2-e-jω/2)=e-j3ω/2sin(2ω)sin(ω/2) (2) x(n)的4点DFT为 X(k)=∑3n=0x(n)Wkn4=∑3n=0e-j2π4kn=1-e-j4·2π4k1-e-j2π4k=4,k=0 0,k=1,2,3 (3) x(n)的8点DFT为 X(k)=∑7n=0x(n)Wkn8=∑3n=0e-j2π8kn=e-j38πksinπ2ksinπ8k,k=0,1,…,7 (4) x(n)的16点DFT为 X(k)=∑15n=0x(n)Wkn16=∑3N=0e-j2π16kn=e-j316πksinπ4ksinπ16k,k=0,1,…,15 图3.2.1给出了x(n)的DTFT和DFT幅频特性曲线示意图。可以看出,N点DFT相当于DTFT在频域[0,2π)上的N点等间隔采样。N决定了谱线的根数,N越大,谱线越密,可观察到的细节越多,|X(k)|的包络也越接近于|X(ejω)|。 图3.2.1R4(n)的DTFT与DFT的关系 1. DFT的矩阵表示 离散傅里叶变换对体现了x(n)和X(k)的时频域映射关系,式(3.2.1)和式(3.2.2)表明,X(k)可由x(n)线性表示; 反之,x(n)也可由X(k)线性表示。许多DFT应用场合中,为了简洁表述信号处理过程,DFT常采用矩阵形式。 设N个x(n)和X(k)分别构成列向量x和X,即 x=x(0) x(1) ︙ x(N-1),X=X(0) X(1) ︙ X(N-1) 则DFT的矩阵形式表示为 X=WNx(3.2.4) 式中 WN=111…1 1W1NW2N…WN-1N 1W2NW4N…W2(N-1)N ︙︙︙WijN︙ 1WN-1NW2(N-1)N…W(N-1)(N-1)N (3.2.5) 为N×N矩阵,其i+1行j+1列的元素为WijN(i,j=0,1,…,N-1); 并且满足: WN=WTN,WNWHN=NIN,其中(·)T表示转置,(·)H表示共轭转置,IN为N×N单位矩阵。 类似地,DFT的矩阵形式表示为 x=1NWHNX(3.2.6) 联立式(3.2.4)和式(3.2.6),可得 x=1NWHNWNx =1NNINx=x(3.2.7) 上式表明,依次经过DFT和IDFT运算后,将恢复出原始序列。当DFT点数大于序列x(n)长度时,对列向量x补零增加其长度。 例3.2.2设长度M=4的有限长序列x(n)=R4(n),利用DFT的矩阵形式求x(n)的4点DFT。 解: DFT点数N=M=4,无须对列向量x补零。DFT的矩阵形式计算为 X=WNx 式中 x=1 1 1 1 ,X=X(0) X(1) X(2) X(3),WN=1111 1W14W24W34 1W24W44W64 1W34W64W94=1111 1-j-1j 1-11-1 1j-1-j 因此 X=X(0) X(1) X(2) X(3)=1111 1-j-1j 1-11-1 1j-1-j·1 1 1 1=4 0 0 0 对比例3.2.1中4点DFT计算结果,两者是一致的。 2. DFT和Z变换、DTFT的关系 设序列x(n)长度为N,n=0,1,…,N-1,其Z变换、DTFT和N点DFT分别为 X(z)=ZT[x(n)]=∑N-1n=0x(n)z-n X(ejω)=DTFT[x(n)]=∑N-1n=0x(n)e-jωn X(k)=DFT[x(n)]=∑N-1n=0x(n)WknN,k=0,1,…,N-1 可以看出,三种变换的关系: X(k)=X(z)z=ej2πNk,k=0,1,…,N-1(3.2.8) X(k)=X(ejω)ω=2πNk,k=0,1,…,N-1(3.2.9) 上述两个关系式说明: N点离散傅里叶变换是Z变换在单位圆上的N点等间隔采样,也是离散时间傅里叶变换在频域[0,2π)上的N点等间隔采样,采样间隔为2π/N。 3. DFT和DFS的关系 有限长序列x(n)和X(k)构成一个离散傅里叶变换对,周期序列x~(n)和X~(k)构成DFS变换对,它们在时域和频域上都是离散的。若将有限长序列x(n)以周期N进行周期延拓后,将会形成一个周期序列x~(n),x~(n)通过DFS可得到X~(k)。因此,DFT和DFS的关系可以表示为 X(k)=∑N-1n=0x(n)WknNRN(k)=X~(k)RN(k)(3.2.10) x(n)=1N∑N-1k=0X(k)W-knNRN(n)=x~(n)RN(n)(3.2.11) 也就是说,DFT可以看作DFS在主值区间上的变换,序列x(n)和X(k)分别对应x~(n)和X~(k)的主值区间,x(n)和X(k)都隐含着周期性,周期均为N,即X(k+lN)=X(k),x(n+lN)=x(n)。这就是DFT的隐含周期特性。 根据DFT和IDFT的定义,利用性质Wk+lNN=WkN,l为整数,也可以推导得到隐含周期特性。即 X(k+lN)=∑N-1n=0x(n)W(k+lN)nN =∑N-1n=0x(n)WknN=X(k),k=0,1,…,N-1 (3.2.12) x(n+lN)=1N∑N-1k=0X(k)W-k(n+lN)N =1N∑N-1k=0X(k)W-knN=x(n),n=0,1,…,N-1(3.2.13) 3.2.2DFT的性质与定理 1. 线性特性 设x1(n)和x2(n)是有限长序列,长度分别为N1和N2,取N=max[N1,N2],计算N点DFT。令X1(k)=DFT[x1(n)],X2(k)=DFT[x2(n)],若y(n)=ax1(n)+bx2(n),a、b为常数,则有 Y(k)=DFT[ax1(n)+bx2(n)]=aX1(k)+bX2(k),k=0,1,…,N-1 (3.2.14) 需要说明的是,当N1和N2不相等时,需要对x1(n)和x2(n)分别补N-N1、N-N2个零,使序列y(n)长度增加到N。 2. 循环移位特性 1) 序列的循环移位 图3.2.2循环移位过程示意图 (N=6,m=2) 设x(n)是长度为N的有限长序列,其循环移位定义为 y(n)=x((n+m))NRN(n)(3.2.15) 整个循环移位过程可以分为三步: (1) 周期延拓: 将x(n)以N为周期进行周期延拓得到x~(n)=x((n))N。 (2) 移位: 将x~(n)左移m位得到x~(n+m)。 (3) 取主值区间: 取x~(n+m)的主值序列,得到循环移位序列y(n)。 序列x(n)及其循环移位过程如图3.2.2所示。显然,y(n)仍是长度为N是有限长序列。由图可见,循环移位的实质是将x(n)左移m位,而移出主值区间[0,N-1]的序列值又依次从右侧进入主值区间,因此称为“循环移位”。 2) 时域循环移位定理 设x(n)是长度为N的有限长序列,X(k)=DFT[x(n)],y(n)为x(n)的循环移位,即y(n)=x((n+m))NRN(n),那么,y(n)的DFT为 Y(k)=DFT[y(n)]=W-kmNX(k),k=0,1,…,N-1(3.2.16) 证明: Y(k)=DFT[y(n)]=∑N-1n=0x((n+m))NRN(n)WknN=∑N-1n=0x((n+m))NWknN 令n+m=n′,则有n=n′-m,可得 Y(k)=∑N-1+mn′=mx((n′))NWk(n′-m)N=W-kmN∑N-1+mn′=mx((n′))NWkn′N 由于上式中求和项x((n′))NWkn′N以N为周期,所以在任一周期上的求和结果相同。不妨将求和区间改为主值区间[0,N-1],可得 Y(k)=W-kmN∑N-1n′=0x((n′))NWkn′N=W-kmN∑N-1n′=0x(n′)Wkn′N =W-kmNX(k),k=0,1,…,N-1 由于W-kmN幅度为1,因此,序列经过循环移位后,其DFT的幅度特性保持不变,仅影响相位特性。 3) 频域循环移位定理 设x(n)是长度为N的有限长序列,X(k)=DFT[x(n)],Y(k)为X(k)的循环移位,即 Y(k)=X((k+l))NRN(k) 则 y(n)=IDFT[Y(k)]=WnlNx(n),n=0,1,…,N-1(3.2.17) 证明: y(n)=IDFT[Y(k)]=1N∑N-1k=0X((k+l))NRN(k)W-knN=1N∑N-1k=0X((k+l))NW-knN 令k+l=k′,则有k=k′-l,可得 y(n)=1N∑N-1+lk′=lX((k′))NW-(k′-l)nN=WnlN1N∑N-1+lk′=lX((k′))NW-k′nN 同理,上式中求和项X((k′))NW-k′nN以N为周期,可将求和区间改为主值区间[0,N-1],可得 y(n)=WnlN1N∑N-1k′=0X((k′))NW-k′nN=WnlN1N∑N-1k′=0X(k′)W-k′nN =WnlNx(n),n=0,1,…,N-1 3. 循环卷积定理 1) 循环卷积的定义和计算 设有限长序列x1(n)和x2(n)长度分别为N1和N2,N=max[N1,N2],则x1(n)和x2(n)循环卷积定义为 x(n)=x1(n)*○x2(n)=∑N-1m=0x1(m)x2((n-m))NRN(n)(3.2.18) 或 x(n)=∑N-1m=0x2(m)x1((n-m))NRN(n)(3.2.19) 与循环卷积相比较,线性卷积表达式为 x1(n)*x2(n)=∑∞m=-∞x1(m)x2(n-m)(3.2.20) 两者的区别: 一是卷积对象不同,循环卷积针对有限长序列,线性卷积无明确要求,可以是有限长序列或者无限长序列; 二是求和区间不同,循环卷积只在主值区间[0,N-1]求和,线性卷积的求和区间则是(-∞,+∞)。 观察循环卷积的定义式(3.2.18)或式(3.2.19),可以看出,循环卷积的计算步骤分为循环翻转、循环移位和乘累加。不妨以式(3.2.18)为参考,整个循环卷积过程描述如下: (1) 循环翻转: 将序列x2(m)进行周期为N的周期延拓,得到x2((m))N,再以纵轴为中心,左右翻转得到x2((-m))N,取主值序列得到x2((-m))NRN(m),该序列称为x2(m)的循环翻转序列。 (2) 循环移位: 将x2((-m))NRN(m)向右循环移位n,形成x2((-(m-n)))NRN(m),即为x2((n-m))NRN(m)。 (3) 乘累加: 将x1(m)和x2((n-m))NRN(m)相乘,并在区间[0,N-1]上对m求和,得到x(n)值。当n=0,1,2,…,N-1时,即可获得x1(n)和x2(n)的循环卷积x(n)。 例3.2.3设序列x1(n)=[1,1,1,1,0,0,0,0],x2(n)=[0,0,1,1,1,1,0,0],n=0~7,试画出N=8时两个序列的循环卷积示意图。 解: 根据循环卷积定义式(3.2.18),序列x1(n)和x2(n)的循环卷积过程如图3.2.3所示。其中图3.2.3(b)为x2(m)的周期延拓序列,图3.2.3(c)为x2(m)的循环翻转序列,图3.2.3 (d)和(e)为循环移位序列,图3.2.3 (f)为最终的循环卷积结果。 图3.2.3循环卷积过程示意图 2) 时域循环卷积定理 设有限长序列x1(n)和x2(n)长度分别为N1、N2,N=max[N1,N2],x1(n)和x2(n)的N点DFT为 X1(k)=DFT[x1(n)],X2(k)=DFT[x2(n)],k=0,1,…,N-1 若X(k)=X1(k)X2(k),则有 x(n)=IDFT[X(k)]=∑N-1m=0x1(m)x2((n-m))NRN(n),n=0,1,…,N-1(3.2.21) 证明: 对式(3.2.21)两边进行DFT,可得 X(k)=DFT[x(n)]=∑N-1n=0∑N-1m=0x1(m)x2((n-m))NRN(n)WknN =∑N-1m=0x1(m)∑N-1n=0x2((n-m))NWknN 令n-m=n′,则有n=n′+m,可得 X(k)=∑N-1m=0x1(m)∑N-1-mn′=-mx2((n′))NWk(n′+m)N =∑N-1m=0x1(m)WkmN∑N-1-mn′=-mx2((n′))NWkn′N 由于求和项x2((n′))NWkn′N以N为周期,可将求和区间改为主值区间[0,N-1],得到 X(k)=∑N-1m=0x1(m)WkmN∑N-1n′=0x2((n′))NWkn′N =X1(k)X2(k),k=0,1,…,N-1 由于式(3.2.21)表示序列x1(n)和x2(n)的时域循环卷积,因此有 x(n)=IDFT[X(k)]=x1(n)*○x2(n)=x2(n)*○x1(n)(3.2.22) 式(3.2.21)或式(3.2.22)表明,两个有限长序列的循环卷积的DFT就是两个序列DFT的乘积。该特性类似于Z变换和DTFT性质中的时域卷积定理,不同的是,时域卷积定理针对一般序列的线性卷积,这里是对有限长序列的循环卷积,因此该性质称为时域循环卷积定理。 3) 频域循环卷积定理 若有限长序列x1(n)和x2(n)的长度均为N,序列x(n)=x1(n)x2(n),则x(n)的DFT为 X(k)=DFT[x1(n)x2(n)]=1NX1(k)*○X2(k) =1N∑N-1l=0X1(l)X2((k-l))NRN(k),k=0,1,…,N-1(3.2.23) 或 X(k)=1NX2(k)*○X1(k)=1N∑N-1l=0X2(l)X1((k-l))NRN(k),k=0,1,…,N-1(3.2.24) 上述两式可利用循环卷积定义加以证明,具体过程略。 4. 复共轭序列的DFT 设有限长序列x(n)的长度为N,其复共轭序列表示为x*(n),若X(k)=DFT[x(n)],则有 DFT[x*(n)]=X*(N-k),k=0,1,…,N-1(3.2.25) 并且X(N)=X(0)。 证明: 根据DFT的定义可得 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),k=1,2,…,N-1 X*(0),k=0 (3.2.26) 上式中利用WNnN=1。由于DFT具有隐含周期特性且周期为N,X(0)=X(N)。为了简洁表示,上式可以统一表示为DFT[x*(n)]=X*(N-k),k=0,1,…,N-1。 类似地,复共轭对称序列的DFT可表示为 DFT[x*(N-n)]=X*(k),k=0,1,…,N-1(3.2.27) 并且x(N)=x(0)。 证明: DFT[x*(N-n)]=∑N-1n=0x*(N-n)WknN 令N-n=n′,则n=N-n′,有 DFT[x*(N-n)]=∑Nn′=1x(n′)W-k(N-n′)N*=∑Nn′=1x(n′)Wkn′N* 由于DFT具有隐含周期特性,周期为N,x(N)=x(0),且WkNN=Wk·0N=1,因此 DFT[x*(N-n)]=∑N-1n′=0x(n′)Wkn′N*=X*(k),k=0,1,…,N-1 5. 帕塞瓦尔定理 若有限长序列x(n)和y(n)的长度均为N,X(k)=DFT[x(n)],Y(k)=DFT[y(n)],则 ∑N-1n=0x(n)y*(n)=1N∑N-1k=0X(k)Y*(k)(3.2.28) 证明: ∑N-1n=0x(n)y*(n)=∑N-1n=0x(n)1N∑N-1k=0Y(k)W-knN* =1N∑N-1k=0Y*(k)∑N-1n=0x(n)WknN =1N∑N-1k=0X(k)Y*(k) 若y(n)=x(n),则有 ∑N-1n=0x(n)x*(n)=1N∑N-1k=0X(k)X*(k) 即 ∑N-1n=0|x(n)|2=1N∑N-1k=0|X(k)|2(3.2.29) 上式就是有限长序列DFT的帕塞瓦尔定理。对比第2章中序列DTFT的帕塞瓦尔定理可以看出,有限长序列在频域的能量有连续积分、离散求和两种形式,两者均等于时域的能量。 6. 共轭对称性 在第2章中已经详细讨论了离散时间傅里叶变换的对称性,即序列关于n=0对称,DTFT关于ω=0对称。离散傅里叶变换也具有类似的对称性,只不过DFT涉及的序列x(n)和X(k)均是有限长,且定义区间[0,N-1]。因此,DFT对称性是指关于N/2对称。下面讨论DFT共轭对称特性。 1) 有限长序列和DFT的共轭对称性 令xep(n)表示有限长共轭对称序列,xop(n)表示有限长共轭反对称序列,则xep(n)和xop(n)满足如下关系式: xep(n)=x*ep(N-n) xop(n)=-x*op(N-n),0≤n≤N-1(3.2.30) 式中,n=1与n=N-1对称,n=2与n=N-2对称; 由于DFT隐含周期特性,n=0与n=N序列值相同,因而xep(0)为实数,xop(0)为纯虚数。为区别第2章中共轭对称序列xe(n)和共轭反对称序列xo(n),这里针对有限长特别引入下标(·)p,蕴含隐含周期性之义。 与DTFT的对称性类似,有限长序列的共轭对称可以理解为在偶对称的基础上引入共轭,共轭反对称可以理解为奇对称的基础上引入共轭。对于有限长序列,共轭对称、共轭反对称、偶对称、奇对称均关于N/2对称。 图3.2.4(a)和(b)分别给出了N为偶数和奇数条件下共轭对称序列和共轭反对称序列的示意图,图中“*”表示序列取复共轭。 图3.2.4N取8和7时共轭对称、共轭反对称序列示意图 与DTFT的对称性相同,有限长序列可以分解为共轭对称序列和共轭反对称序列两部分,即 x(n)=xep(n)+xop(n),0≤n≤N-1(3.2.31) 为了得到xep(n)和xop(n)的表达式,将上式中n替换为N-n,取复共轭,可得 x*(N-n)=x*ep(N-n)+x*op(N-n)=xep(n)-xop(n)(3.2.32) 将式(3.2.31)、式(3.2.32)分别相加和相减,可得 xep(n)=12[x(n)+x*(N-n)] xop(n)=12[x(n)-x*(N-n)](3.2.33) 同理,X(k)也可以分解为共轭对称和共轭反对称两部分,分别记为Xep(k)、Xop(k),那么 X(k)=Xep(k)+Xop(k)(3.2.34) Xep(k)=12[X(k)+X*(N-k)] Xop(k)=12[X(k)-X*(N-k)](3.2.35) 2) 有限长序列和DFT的共轭对称性的对应关系 (1) 序列实虚部DFT的共轭对称性。 设有限长序列x(n)是一个复序列,其实部和虚部分别为xR(n)、xI(n),即 x(n)=xR(n)+jxI(n)(3.2.36) 式中 xR(n)=Re[x(n)]=12[x(n)+x*(n)],jxI(n)=jIm[x(n)]=12[x(n)-x*(n)] 利用复共轭序列的DFT特性,即式(3.2.25),对上述两式进行DFT,可得 DFT[xR(n)]=12DFT[x(n)+x*(n)]=12[X(k)+X*(N-k)]=Xep(k)(3.2.37) DFT[jxI(n)]=12DFT[x(n)-x*(n)]=12[X(k)-X*(N-k)]=Xop(k)(3.2.38) 根据DFT的线性特性,上述两式相加,正好得到式(3.2.34)。因此,序列实部的DFT对应于序列DFT的共轭对称部分,j和虚部的DFT对应于序列DFT的共轭反对称部分。 (2) 序列DFT实虚部的共轭对称性。 若将有限长序列x(n)分解为共轭对称序列xep(n)和共轭反对称序列xop(n),即式(3.2.31),分别对xep(n)和xop(n)求DFT,利用式(3.2.33)和式(3.2.27),可得 DFT[xep(n)]=12DFT[x(n)+x*(N-n)]=12[X(k)+X*(k)]=Re[X(k)](3.2.39) DFT[xop(n)]=12DFT[x(n)-x*(N-n)]=12[X(k)-X*(k)]=jIm[X(k)](3.2.40) 上述两式相加,可得 X(k)=DFT[x(n)]=DFT[xep(n)]+DFT[xop(n)]=XR(k)+jXI(k) (3.2.41) 式中,XR(k)和XI(k)分别表示X(k)的实部和虚部。 可见,序列共轭对称部分的DFT对应于序列DFT的实部,而共轭反对称部分的DFT对应于序列DFT的虚部和j相乘。 综上所述,可以归纳DFT的共轭对称性: 不管是序列还是DFT,实部始终对应于变换(DFT或IDFT)的共轭对称部分,j乘虚部对应于变换的共轭反对称部分。表3.2.1给出了有限长序列和DFT的共轭对称性的对应关系。 表3.2.1有限长序列和DFT的共轭对称性的对应关系 共轭对称性对应关系 x(n)DFTIDFTX(k) X(k)IDFTDFTx(n) 实部共轭对称 xR(n)DFTIDFTXep(k) XR(k)IDFTDFTxep(n) j乘虚部共轭反对称 jxI(n)DFTIDFTXop(k) jXI(k)IDFTDFTxop(n) 3) 共轭对称性的应用 在实际数字信号处理系统中,常通过A/D采样得到实数序列,进行DFT分析处理; 或者构造X(k)后利用IDFT生成实数信号,并通过D/A进行发送。下面以实序列为例,讨论如何利用共轭对称性来分析实序列DFT的特点,以及如何计算实序列DFT,从而达到减少DFT计算量、提高计算效率的目的。 设x(n)是长度为N的实数序列,X(k)=DFT[x(n)] (1) 实序列无虚部,DFT无共轭反对称部分,即X(k)共轭对称: X(k)=X*(N-k)(3.2.42) (2) 若x(n)=x(N-n),x(n)为共轭对称序列,那么X(k)只有实部,无虚部,并且关于N/2共轭对称,因此X(k)为实偶对称: X(k)=X(N-k) (3.2.43) (3) 若x(n)=-x(N-n),x(n)为共轭反对称序列,那么X(k)只有j和虚部,无实部,并且关于N/2共轭对称,因此X(k)为纯虚奇对称: X(k)=-X(N-k)(3.2.44) 当N为偶数时,只需要计算前面N/2+1点DFT值; 当N为奇数时,计算前面(N+1)/2点DFT值,其他值按照式(3.2.42)即可得到。如X(N-1)=X*(1),X(N-2)=X*(2),这样可减少近一半运算量。 利用共轭对称性,通过一个DFT可计算两个实序列的DFT。由于序列实部、j乘虚部的DFT分别对应于序列DFT的共轭对称和共轭反对称部分,如果将两个实序列合成一个复序列,通过计算复序列的N点DFT,就可同时得到两个实序列的N点DFT。 设x1(n)和x2(n)表示长度为N的两个实序列,将x1(n)和x2(n)分别作为实部和虚部,构造新的复序列x(n)如下: x(n)=x1(n)+jx2(n)(3.2.45) 令X(k)=DFT[x(n)],X1(k)=DFT[x1(n)],X2(k)=DFT[x2(n)],k=0,1,…,N-1,根据序列实虚部DFT的对称性,即式(3.2.37)和式(3.2.38),可得 X1(k)=Xep(k)=12[X(k)+X*(N-k)](3.2.46) X2(k)=1jXop(k)=12j[X(k)-X*(N-k)](3.2.47) 有限长序列的DFT的基本性质见表3.2.2。 表3.2.2有限长序列的DFT的基本性质 序号 序列 离散傅里叶变换 性质 1 x(n) X(k)=X(k+lN),l为整数 隐含周期性 2 ax1(n)+bx2(n) aX1(k)+bX2(k),a,b为常数 线性 3 x((n+m))NRN(n) W-kmNX(k) 4 WnlNx(n) X((k+l))NRN(k) 循环移位定理 5 ∑N-1m=0x1(m)x2((n-m))NRN(n) X1(k)X2(k) 6 x1(n)x2(n) 1N∑N-1l=0X1(l)X2((k-l))NRN(k) 循环卷积定理 7 x*(n) X*(N-k) 8 x*(N-n) X*(k) 复共轭序列 9 Re[x(n)] 12[X(k)+X*(N-k)]=Xep(k) 10 jIm[x(n)] 12[X(k)-X*(N-k)]=Xop(k) 11 12[x(n)+x*(N-n)]=xep(n) Re[X(k)] 12 12[x(n)-x*(N-n)]=xop(n) jIm[X(k)] 共轭对称性 13 x(n)实序列 X(k)共轭对称,X(k)=X*(N-k) 14 x(n)实偶对称,x(n)=x(N-n) X(k)实偶对称,X(k)=X(N-k) 15 x(n)实奇对称,x(n)=-x(N-n) X(k)虚奇对称,X(k)=-X(N-k) 实序列 共轭对称性 16 ∑N-1n=0x(n)y*(n)=1N∑N-1k=0X(k)Y*(k) 17 ∑N-1n=0|x(n)|2=1N∑N-1k=0|X(k)|2 帕塞瓦尔定理 3.3频域采样与内插 由时域采样定理可知,若采样频率大于或等于连续时间信号最高频率的2倍,那么可由离散时间信号恢复原来的连续时间信号,这个恢复过程称为内插; 对应地,也存在频域采样定理。根据DFT与Z变换、DTFT的关系,DFT可看作Z变换在单位圆上或者DTFT在[0,2π)上的等间隔采样,即DFT实现了频域采样。那么,能否由离散的频域采样值恢复出连续的频谱呢?如果可以,条件是什么?恢复连续频谱的内插公式又是什么形式?下面围绕这些问题进行讨论。 3.3.1频域采样定理 假设任意长序列x(n)满足绝对可和条件,则DTFT和Z变换存在,分别表示为X(ejω)和X(z)。由于X(z)收敛域中包含单位圆,若在单位圆上对X(z)进行N点等间隔采样,则可得 X(k)=X(z)|z=ej2πNk=∑∞n=-∞x(n)e-j2πNkn,0≤k≤N-1(3.3.1) 用DTFT表示,有 X(k)=X(ejω)|ω=2πNk,0≤k≤N-1(3.3.2) 式(3.3.2)表明: X(k)相当于在频率区间[0,2π)上对DTFT进行N点等间隔采样。 由于X(k)为离散的频域采样值,可以看作是一个长度为N的有限长序列xN(n)的DFT,即 xN(n)=IDFT[X(k)],0≤n≤N-1(3.3.3) 上述能否由离散X(k)恢复出连续频谱的问题,相当于能否由X(k)恢复出原始序列x(n); 而X(k)与xN(n)存在DFT变换关系。因此,频域上的问题完全可以转化时域上的问题,即xN(n)能否恢复出x(n)?这需要探讨有限长序列xN(n)和原始序列x(n)的关系来确定。下面推导xN(n)和x(n)的关系,并进一步导出频域采样定理。 由于原始序列x(n)长度没有指定为有限长或无限长,而xN(n)为有限长序列,无法直接建立两者之间的关系。但是,从DFT和DFS的关系可知,xN(n)和X(k)具有隐含周期性,若分别以周期N进行周期延拓得到x~N(n)和X~(k),则X(k)可以视为周期延拓序列x~N(n)的离散傅里叶级数X~(k)的主值序列。此时,从无限长周期延拓序列的主值序列角度来看待xN(n),就便于建立与x(n)的关系。即有 xN(n)=x~N(n)RN(n),X(k)=X~(k)RN(k) 根据IDFS,有 x~N(n)=xN((n))N=IDFS[X~(k)]=1N∑N-1k=0X~(k)W-knN =1N∑N-1k=0X(k)W-knN(3.3.4) 将式(3.3.1)代入式(3.3.4)可得 x~N(n)=1N∑N-1k=0∑∞m=-∞x(m)WkmNW-knN =∑∞m=-∞x(m)1N∑N-1k=0Wk(m-n)N(3.3.5) 式中 ∑N-1k=0Wk(m-n)N=N,m=n+rN 0,m≠n+rN,r为整数 由于m,n均没有限定,r=-∞~∞,因此有 x~N(n)=∑∞r=-∞x(n+rN)(3.3.6) 式(3.3.6)表明: 由X~(k)得到的周期序列x~N(n)是原始序列x(n)以N为周期的周期延拓序列。由时域采样定理可知,时域的采样造成频域的周期延拓,延拓周期为采样频率。这里可以对应地看到,DTFT频域采样会造成时域序列的周期延拓,延拓周期为采样点数。这正是傅里叶变换在时域和频域对称关系的反映。 取x~N(n)的主值序列可得 xN(n)=∑∞r=-∞x(n+rN)RN(n)(3.3.7) 式(3.3.7)表明了有限长序列xN(n)和原始序列x(n)的关系,即xN(n)是x(n)以N为周期进行周期延拓后的主值序列。xN(n)和x(n)的关系可分为以下两种情况: (1) 若x(n)是无限长序列,无论N取何值,周期延拓都会引起时域混叠,不可能从x~N(n)中提取主值区间来不失真地恢复出x(n)。即xN(n)和x(n)始终存在误差,只是随着N的增大,频域采样越密,误差越小。 (2) 若x(n)是有限长序列,长度为M,显然,如果采样点数N≥M,周期延拓将不会产生混叠现象,可以无失真地恢复x(n),即有x(n)=xN(n)。如果N<M,将有时域混叠,无法不失真地恢复x(n)。此时,只有通过增大N,满足N≥M条件。 通过以上分析可以得出以下结论: 对于长度为M的序列x(n),只有当频域采样点数N≥M时,才可由X(ejω)的频域采样值X(k)无失真地恢复x(n),即 x(n)=xN(n)=1N∑N-1k=0X(k)W-knN(3.3.8) 从而可以无失真恢复出x(n)的连续频域特性X(ejω)或X(z),这就是频域采样定理。 3.3.2频域内插公式 根据频域采样定理,既然由频域采样X(k)可以无失真地恢复出序列x(n),而x(n)进行Z变换可得到X(z),因此,可以方便地由X(k)得到X(z)或者X(ejω),相当于由离散的频域采样值表示整个X(z)函数以及连续的频率响应X(ejω)。 假设序列x(n)的长度为N,0≤n≤N-1,X(k)表示X(z)在单位圆的N点等间隔采样,对式(3.3.8)进行Z变换,可得 X(z)=∑∞n=-∞x(n)z-n=∑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(3.3.9) 由于W-kNN=1,因此 X(z)=1N∑N-1k=0X(k)1-z-N1-W-kNz-1(3.3.10) 上式即为用频域采样值X(k)表示X(z)的频域内插公式,该公式将用在第7章中,通过离散的频域采样值设计FIR数字滤波器。 令内插函数为 Φk(z)=1N1-z-N1-W-kNz-1(3.3.11) 则式(3.3.10)的频域内插公式可以重新表述为 X(z)=∑N-1k=0X(k)Φk(z)(3.3.12) 下面讨论频率响应。将z=ejω代入式(3.3.11),可得 Φk(ω)=1N1-e-jωN1-e-j(ω-2πk/N)=1Nsin(ωN/2)sin[(ω-2πk/N)/2]e-jN-12ω+πkN =1Nsin[(ω-2πk/N)N/2]sin[(ω-2πk/N)/2]e-jN-12ω-2πkN=Φω-2πNk(3.3.13) 式中,Φ(ω)与k无关,称为频率响应的内插函数,表达式为 Φ(ω)=1Nsin(ωN/2)sin(ω/2)e-jωN-12(3.3.14) 联立式(3.3.12)和式(3.3.13),频率响应的内插公式可表示为 X(ejω)=∑N-1k=0X(k)Φω-2πNk 3.4离散傅里叶变换的应用 DFT是数字信号处理的重要工具,在数字通信、语音处理、图像处理、雷达等领域得到广泛应用。DFT有快速算法——快速傅里叶变换,能够大大降低运算量,使得DFT在工程实践中应用更加广泛和高效。本节将介绍DFT的两类典型应用: 线性卷积计算以及对连续时间信号和序列进行谱分析。 3.4.1DFT计算线性卷积 在数字信号通过线性时不变系统或者数字滤波器时,其输出等于输入与系统单位冲激响应的线性卷积,如果能够将线性卷积转化为循环卷积,根据DFT的时域循环卷积定理,循环卷积可以用DFT(FFT)来计算,这样就能够用FFT来计算线性卷积,提高运算速度。 下面首先讨论循环卷积和线性卷积的等价条件,然后介绍具体利用DFT计算线性卷积的方法。 1. 线性卷积和循环卷积的等价条件 设h(n)和x(n)均为有限长序列,长度分别为N和M,循环卷积长度L≥max[N,M],则线性卷积和循环卷积分别表示为 yl(n)=h(n)*x(n)=∑N-1m=0h(m)x(n-m)(3.4.1) yc(n)=h(n)*○x(n)=∑L-1m=0h(m)x((n-m))LRL(n)(3.4.2) 式中,x((n))L表示x(n)的周期延拓,表达式为 x((n))L=∑∞q=-∞x(n+qL)(3.4.3) 将上式代入式(3.4.2),可得 yc(n)=∑N-1m=0h(m)∑∞q=-∞x(n-m+qL)RL(n) =∑∞q=-∞∑N-1m=0h(m)x(n+qL-m)RL(n)(3.4.4) 对比式(3.4.1),可以看出 ∑N-1m=0h(m)x(n+qL-m)=yl(n+qL)(3.4.5) 因此,有 yc(n)=∑∞q=-∞yl(n+qL)RL(n)(3.4.6) 上式表明: 循环卷积yc(n)相当于线性卷积yl(n)周期延拓后的主值序列,周期正好是循环卷积长度L。而线性卷积yl(n)的长度为N+M-1,因此,只有当循环卷积长度L≥N+M-1时,yl(n)进行周期延拓才无混叠,此时主值序列就是yl(n),式(3.4.6)即为yc(n)=yl(n),两者等价。由此得出结论: 线性卷积和循环卷积的等价条件是循环卷积长度大于或等于线性卷积的长度,即 L≥N+M-1(3.4.7) 图3.4.1给出了循环卷积和线性卷积的对比图。h(n)和x(n)分别为4点和5点矩形序列,线性卷积长度N+M-1=8,如图3.4.1(c)所示,循环卷积长度L为6、8、10时如图3.4.1(d)~(f)所示。可以看出,只有当L≥8时,循环卷积和线性卷积的结果才相同。 图3.4.1循环卷积和线性卷积的对比图 2. 线性卷积的DFT计算方法 基于线性卷积和循环卷积的等价条件,利用DFT计算线性卷积的框图如图3.4.2所示。h(n)和x(n)需要补零达到循环卷积长度L,图中输出y(n)为 y(n)=h(n)*x(n)=h(n) x(n),L≥N+M-1(3.4.8) 图3.4.2DFT计算线性卷积的框图 在实际应用中,DFT和IDFT通常采用快速傅里叶变换来实现,能够大大减少运算量。当序列h(n)和x(n)的长度比较长,且相差不大时,与直接计算线性卷积耗用的乘法和加法运算相比,采用FFT计算线性卷积的运算量较低,因而也称为快速卷积法。 但是,快速卷积法在某些场合下并不一定“快速”。假设一长一短两个序列进行线性卷积,长度相差较大,如MN,如果选择L≥N+M-1进行快速卷积,则短序列需要补充很多个零,而且FFT算法的点数比较大,运算效率比较低,与直接计算线性卷积相比,运算量不一定小。此时,直接计算或许是一个更好的选择。 此外,某些场合下长序列的长度不固定,甚至接近无限长,如语音信号、数字通信信号等,往往要求持续接收和处理,强调实时性。在这种情况下,不能直接套用快速卷积方法,比较好的解决思路是利用卷积的线性性质,将长序列分段,每一段序列与短序列分别进行卷积,再合成最终结果。这就是分段处理的思想,具体方法包括重叠相加法和重叠保留法。 3. 重叠相加法 设序列h(n)长度为N,x(n)为无限长序列,n≥0,将x(n)进行均匀分段,每段长度为M,则x(n)可以表示为 x(n)=∑∞k=0xk(n)(3.4.9) 式中: 第k段序列为xk(n)=x(n)·RM(n-kM)。 那么,h(n)与x(n)的线性卷积计算为 y(n)=h(n)*x(n)=h(n)*∑∞k=0xk(n) =∑∞k=0h(n)*xk(n)=∑∞k=0yk(n)(3.4.10) 式中: yk(n)=h(n)*xk(n)表示第k段序列的线性卷积结果,其起始时刻为n=kM,长度为N+M-1。 式(3.4.10)表明: 计算有限长序列h(n)与无限长序列x(n)的线性卷积时,可先将x(n)进行分段,计算每一段xk(n)与h(n)的线性卷积,再将分段卷积结果yk(n)重叠相加即可。分段卷积可以采用快速卷积方法或者直接进行计算。 图3.4.3给出了线性卷积的重叠相加法示意图。从图中可以看出,由于分段卷积结果yk(n)长度为N+M-1,而起始时刻为n=kM,因此,yk+1(n)与yk(n)必然有N-1个点发生重叠,必须把yk+1(n)的重叠部分加到yk(n)上,才能得到完整的线性卷积序列y(n)。也就是说: 当y0(n)计算完毕后,能输出M个值; y1(n)计算完毕后,能输出2M个值; 当yk(n)计算完毕后,能输出kM个值,其后续N-1个值等待yk+1(n)重叠相加后才能确定。因此,该卷积方法也称为“重叠相加法”。 图3.4.3重叠相加法计算线性卷积示意图 4. 重叠保留法 与重叠相加法的分段卷积叠加思想不同,重叠保留法是由分段卷积结果衔接而成。在按照式(3.4.9)对x(n)进行分段的同时,将输出y(n)也进行均匀分段,每段长度为M,则有 y(n)=∑∞k=0yk(n)(3.4.11) 式中,yk(n)=y(n)RM(n-kM)。 将y(n)线性卷积表达式(3.4.1)代入yk(n)可得 yk(n)=∑N-1m=0h(m)x(n-m)RM(n-kM)(3.4.12) 由于yk(n)自变量取值为n=kM~(k+1)M-1,而h(m)自变量取值为m=0~N-1,因此,对于x(n)而言,其真正参与yk(n)线性卷积的只是变量n-m所对应的序列,即x(n)中自变量取值为kM-(N-1)~(k+1)M-1的一段子序列,不妨记为x′k(n)。令n-m=l,式(3.4.12)可重新表示为 yk(n)=∑kM+M-1l=kM-(N-1)x(l)h(n-l)RM(n-kM)(3.4.13) 若从序列x(n)分段的角度来看,子序列x′k(n)按自变量取值范围可分为kM-(N-1)~kM-1以及kM~(k+1)M-1两部分,前一部分来自xk-1(n),后一部分正好为xk(n)。这说明: 对于任一分段序列xk(n),不仅要全部参与yk(n)卷积运算,而且要保留部分样点参与yk+1(n)的卷积运算。因此,这种卷积方法称为“重叠保留法”。 图3.4.4给出了线性卷积的重叠保留法示意图。卷积过程描述如下: (1) 将x(n)均匀分段,形成多个长度为M的分段序列xk(n); (2) 保留前一个分段序列xk-1(n)的尾部N-1点,并结合下一个完整的M点分段序列xk(n),形成长度为M+N-1序列x′k(n),与N点序列h(n)进行分段线性卷积。 图3.4.4重叠保留法计算线性卷积示意图 (3) 去掉分段线性卷积结果的头尾各N-1点,将中间M点作为y(n)分段序列yk(n)。实际上,yk(n)就是h(n)与x′k(n)卷积时完全重叠时产生的计算结果。 (4) 将分段序列yk(n)直接衔接起来,即可得到输出y(n)。 3.4.2DFT进行谱分析 DFT是计算机分析信号与系统频域特性的主要工具,其主要应用之一就是对信号进行谱分析,可用于频偏估计、干扰抑制等场合。谱分析是指信号的傅里叶变换,获得信号的频谱。由于实际信号可能为连续时间信号或者序列,对于连续时间信号,通过时域采样才可用DFT进行谱分析。下面分别针对连续时间信号和序列,介绍如何利用DFT进行谱分析,并探讨影响谱分析效果的因素。 1. 连续时间信号的谱分析 1) 谱分析原理 连续时间信号的频谱与信号周期性密切相关,对于非周期性的连续时间信号,其傅里叶变换存在,即频谱函数存在且是连续的; 而对于周期性的连续时间信号,其傅里叶变换为冲激函数,通常采用傅里叶级数来表示频谱,并且频谱是离散的。利用DFT进行频谱分析,这里重点针对非周期性的连续时间信号,通过DFT离散频谱来表征连续的频谱函数。 假设连续时间非周期信号为xa(t),为了突出频谱函数与频率f的关系,直接采用Xa(jf)而非Xa(jΩ)来表示频谱函数,其中Ω=2πf。那么xa(t)与Xa(jf)构成的傅里叶变换对表示为 Xa(jf)=∫∞-∞xa(t)e-j2πftdt(3.4.14) xa(t)=∫∞-∞Xa(jf)ej2πftdf(3.4.15) 利用DFT对xa(t)进行频谱分析时,需要执行以下三个步骤: (1) 时域采样: 对连续时间信号xa(t)进行等间隔时域采样,得到离散时间信号,即序列x(n)。 (2) 时域截短: 将序列x(n)截短为N点有限长序列xN(n),截短可视为与有限长序列w(n)相乘。 (3) 离散傅里叶变换: 计算xN(n)的DFT得到XN(k),由于XN(k)是频域离散的,可理解为连续频谱的频域采样。 从频域角度来看,连续时间信号的DFT谱分析实际上依次通过序列的DTFT、截短序列的DTFT以及离散傅里叶变换来实现,这是对连续时间信号频谱的一个逼近过程,整个过程如图3.4.5所示。 图3.4.5连续时间信号DFT谱分析的逼近过程 通过上述三步,得到的XN(k)可作为连续时间信号xa(t)的谱分析结果。那么,现在的问题是: 频域离散的XN(k)是否能够准确地代表连续频谱Xa(jf),XN(k)是否是Xa(jf)的准确采样?要回答上述问题,需要探讨XN(k)与Xa(jf)的关系,即有限长序列xN(n)的DFT和连续时间信号xa(t)的傅里叶变换到底有何关系。 下面围绕三个步骤,详细讨论XN(k)和Xa(jf)之间的关系。 (1) 时域采样。 设连续时间信号xa(t)的最高频率为fc,在对xa(t)进行时域采样时要满足时域采样定理,即采样频率必须大于或等于最高频率的2倍,否则会引起频谱混叠现象。令fs表示采样频率,T表示采样间隔,则有 fs=1/T≥2fc,T≤1/(2fc)(3.4.16) 采样后得到的序列x(n)为 x(n)=xa(t)|t=nT=xa(nT)(3.4.17) x(n)的离散时间傅里叶变换为 X(ejω)=∑∞n=-∞x(n)e-jωn(3.4.18) 通过时域采样,由连续时间信号xa(t)得到序列x(n),从频域角度来看,连续信号频谱Xa(jf)和X(ejω)也存在相互关系。 由第2章中序列DTFT和连续时间信号的傅里叶变换关系可知 X(ejω)=1T∑∞k=-∞XajΩ-jk2πT=1T∑∞k=-∞Xaj2πf-kT(3.4.19) 上式表明: Xa(jf)以周期fs=1/T进行周期延拓后再乘以1/T即可得到X(ejω)。若Xa(jf)频率范围有限,fs≥2fc,那么X(ejω)将无频谱混叠现象,并且 Xa(jf)=TX(ejω),|f|≤fs/2(3.4.20) (2) 时域截短。 从序列x(n)中截取一段长度为N的有限长序列,记为xN(n),可以表示为 xN(n)=x(n)w(n)(3.4.21) 式中,w(n)是一个有限长序列,0≤n≤N-1,称为窗函数。 两个序列相乘就体现为时域截短,相当于对x(n)进行加窗操作,窗内序列保留,窗外序列置0。典型窗函数为矩形窗w(n)=RN(n),此时xN(n)=x(n)RN(n)。窗函数详细内容在第7章中介绍。 根据DTFT的性质,时域相乘对应于频域卷积,因此,xN(n)的DTFT为 XN(ejω)=X(ejω)*W(ejω)=∑N-1n=0xN(n)e-jωn (3.4.22) 式中,W(ejω)为窗函数w(n)的DTFT。 若w(n)足够长或变化平缓,W(ejω)呈现类似Sa(·)函数的尖峰特性,那么XN(ejω)≈X(ejω)。因此,就可利用XN(ejω)替代X(ejω)来逼近xa(t)频谱Xa(jf),即 Xa(jf)≈TXN(ejω)(3.4.23) (3) 离散傅里叶变换。 对有限长序列xN(n)进行N点DFT,可得 XN(k)=∑N-1n=0xN(n)e-j2πNkn,0≤k≤N-1(3.4.24) XN(k)相当于XN(ejω)在数字频率0~2π之间进行N点等间隔采样。基于式(3.4.19)的周期延拓特性,数字频率-π~π对应模拟频率-fs/2~fs/2,数字频率π~2π与-π~0的频谱相同,因此对应到Xa(jf),相当于其在模拟频率-fs/2~fs/2之间进行N点等间隔采样。 设模拟频率的采样间隔为F,则参数fs、N、T和F的关系为 F=fsN=1NT(3.4.25) 由于NT就是有限长序列xN(n)对应的采样时间,不妨记为Tp=NT,那么 F=1Tp(3.4.26) 将f=kF和ω=2πk/N代入式(3.4.23),可以得连续信号的频谱采样。 当0≤k≤N2-1(N为偶数)或0≤k≤N-12(N为奇数)时,有 Xa(jkF)≈TXN(k)(3.4.27) 当N2≤k≤N-1(N为偶数)或N+12≤k≤N-1(N为奇数)时,有 Xa(j(k-N)F)≈TXN(k)(3.4.28) 上述分析表明: 连续时间信号xa(t)频谱函数Xa(jf)可由时域采样和截短后的序列的DFT来逼近,离散的XN(k)并不是Xa(jf)的准确频域采样,只能近似地表征Xa(jf)。 图3.4.6给出了连续时间信号和序列的时域以及频域示意图,通过比较可以看出它们在时域上的相互关系以及频域上的逼近过程,即左侧在时域上呈现采样、截短,而右侧在频域上体现周期延拓、卷积以及频域采样,最终得到连续时间信号的DFT谱分析结果。 图3.4.6连续时间信号的频谱逼近示意图 2) 谱分析参数选择 在对连续时间信号进行谱分析时,主要考虑以下两个方面的参数: (1) 谱分析范围: 信号最高频率fc代表了谱分析范围。根据时域采样定理,谱分析范围受采样频率fs的限制。为了避免频域混叠现象,信号最高频率fc≤fs/2,即采样间隔T≤1/(2fc)。 若信号频谱为无限宽,可以选取占信号总能量一定百分比的频带宽度(|f|<fc)来确定信号最高频率fc,如百分比为90%或98%等。当信号最高频率已经确定时,选择采样频率要满足fs≥2fc。 (2) 频率分辨率: 频域采样间隔F代表了频率分辨率,表示谱分析中能够分辨的最小频率间隔。F越小,谱分析越接近Xa(jf),频率分辨率越高。 当给定频率分辨率要求时,根据式(3.4.25)中F=fs/N=1/Tp,在保证谱分析范围不变(fs不变)的情况下,采样时间Tp和采样点数N必须满足 Tp≥1F N=TpT≥fsF(3.4.29) 若要提高频率分辨率(F减小),需要增加采样时间Tp和采样点数N。 例3.4.1利用DFT对语音信号进行谱分析,要求频率分辨率F≤10Hz,信号最高频率fc=4kHz,试确定: 最小记录时间Tpmin、最大采样间隔Tmax、最少的采样点数Nmin。如果fc不变,要求频率分频率增加1倍,最少的采样点数和最小的记录时间是多少? 解: Tp≥1F=110=0.1(s) 即Tpmin=0.1s,而fs≥2fc,fsmin=2fc=8000Hz,因此 T≤12fc=12×4000=0.125×10-3(s),N≥2fcF=2×400010=800 即Tmax=0.125×10-3s,Nmin=800。当频率分辨率提高1倍时,F=5Hz,那么 Tpmin=15=0.2(s),Nmin=2×40005=1600 在实际应用中,为了使用DFT的快速算法FFT,通常选取N为2的整数幂。此时,若采样频率fs不变,即采样间隔T不变,那么采样点数N可分别选取1024和2048,采样时间Tp相应增大,F值减小,具有更高的频率分辨能力。 2. 序列的谱分析 利用DFT可以对序列进行谱分析,若从连续时间信号时域采样得到序列的角度来看,序列谱分析只相当于连续时间信号谱分析步骤中的后两步或者第三步,这与序列是无限长或有限长有关。设序列为x(n),其离散时间傅里叶变换存在,即X(ejω)=DTFT[x(n)],且是连续频谱。 (1) 若序列为有限长序列,长度为N,则直接进行N点DFT计算,得到X(k)=DFT[x(n)],X(k)相当于X(ejω)在数字频率[0,2π)上的N点等间隔采样。根据频域采样定理和内插公式,X(ejω)可以由X(k)无失真恢复。 (2) 若序列为无限长序列,需要进行时域截短,得到有限长序列,才可计算DFT。对于周期序列,可采用离散傅里叶级数来表征频谱,也可截取一个或多个周期构成有限长序列,利用DFT来分析频谱。 在实际应用场合中,周期序列由模/数转换器对连续时间周期信号进行采样得到,序列中除有用信号成分以外,还可能包括噪声成分。在利用DFT进行谱分析时,噪声会对有用信号频谱产生影响。为了提高噪声条件下周期序列的谱分析效果,可以考虑截取多个周期进行DFT,相当于通过增加时间积累来抑制噪声的影响。下面讨论周期序列多个周期的DFT,并与单周期DFT进行比较。 假设x~(n)是周期为N的周期序列,其主值序列x(n)=x~(n)RN(n),X(k)=DFT[x(n)]。截取x~(n)的m个周期,m为正整数,截取后序列长度M=mN,序列表示为 xM(n)= x~(n)RM(n) 对xM(n)进行M点DFT,可得 XM(k)=DFT[xM(n)]=∑M-1n=0x~(n)e-j2πMkn =∑mN-1n=0x~(n)e-j2πmNkn,k=0,1,…,mN-1 令n=n′+rN,则r=0,1,…,m-1,n′=0,1,…,N-1。那么 XM(k)=∑m-1r=0∑N-1n′=0x~(n′+rN)e-j2π(n′+rN)kmN=∑m-1r=0∑N-1n′=0x(n′)e-j2πn′mNke-j2πmrk =∑m-1r=0Xej2πN·kme-j2πmrk=Xej2πN·km∑m-1r=0e-j2πmrk 由于 ∑m-1r=0e-j2πmkr=m,k/m为整数 0,k/m不为整数 所以 XM(k)=mXkm,k/m为整数 0,k/m不为整数(3.4.30) 上式表明: XM(k)也能表示x~(n)的频谱特性,当k=rm时,XM(rm)=mX(r),相当于单周期DFT的X(r)幅度扩大m倍,而k取其他值时,XM(k)=0。从频域采样角度来看,单周期DFT的X(r)与多周期DFT的XM(rm)的频率是对应的,即2πNr=2πmNmr; 同时,由于多周期DFT的幅度扩大m倍,使得谱分析时具有更好的抗噪声能力,更容易获得原始周期序列的频谱特性。 3. 谱分析的误差来源及改进措施 从前两节讨论可以看出,DFT谱分析实际上是以离散频谱对原始频谱的一个逼近过程,这种逼近可能会带来一定的误差,而误差来源则与谱分析步骤密切有关。对于连续时间信号,涉及时域采样、时域截短和离散傅里叶变换,而序列谱分析涉及离散傅里叶变换,甚至时域截短。下面针对这些步骤讨论DFT谱分析的误差问题。 1) 混叠现象 针对连续时间信号的时域采样步骤,若时域采样未满足采样定理,则会引起频谱混叠现象,混叠出现在数字频率ω=π和模拟频率f=fs/2附近。解决频谱混叠的方法是提高采样频率,使之满足fs≥2fc,通常为fs=(3~5)fc。在采样频率已确定的情况下,可以在时域采样前对连续时间信号进行预滤波,滤除高于fs/2的频率成分,避免混叠现象。值得说明的是,利用DFT进行谱分析,能够观察的最高频率成分为fs/2。 2) 频谱泄漏 针对时域截短步骤,时域截短是为了得到有限长序列xN(n)=x(n)·w(n),序列时域相乘对应频域卷积,即XN(ejω)=X(ejω)*W(ejω)。如果w(n)频谱W(ejω)为单位冲激函数δ(ω),那么卷积后XN(ejω)=X(ejω),频谱不变。但由于窗函数w(n)是有限长序列,W(ejω)不可能为单位冲激函数,如矩形窗w(n)=RN(n)频谱具有1个主瓣和多个旁瓣。这样,X(ejω)与W(ejω)卷积会带来频谱的变化。 图3.4.7频谱泄漏的示意图 图3.4.7给出了单位冲激函数形式的频谱与矩形窗频谱卷积后的示意图。可以看出,原来的离散谱线通过卷积后已成为具有主瓣和旁瓣的频谱,出现了拖尾和展宽。因此,频域卷积一定会造成频谱X(ejω)的“扩散”,包括拖尾和展宽,这种现象称为频谱泄漏。 频谱泄漏使得频谱变得模糊,给频谱分析带来误差。最直接的影响体现在两个方面: 一是使频率分辨率降低,受窗函数频谱主瓣宽度影响,当待分析的两个频率逐步靠近时,其中间部分频谱叠加将会超过两个频率幅度,呈现单峰特性,使得两个频率无法分辨; 二是造成谱间干扰,同样由于窗函数主瓣和旁瓣影响,待分析的相邻频率分量之间存在相互干扰,如果频谱展宽使最高频率超过fs/2,会引起更大范围的频谱混叠。 谱分析的一种典型应用是干扰抑制,如消除某个幅度很高的干扰频率成分,要求是尽可能减轻其频谱泄漏对有用信号的影响; 否则,在消除干扰频率的同时,有用信号的频率成分也被消除过多,损伤过大,或者弱小有用信号会被干扰的泄漏成分所淹没,无法有效识别,因此,必须重视频谱泄漏问题。 减小频谱泄漏有两种方法: 一种方法是增加窗函数w(n)长度N,在获得更长的数据的同时,使W(ejω)主瓣更窄,提高频率分辨率; 另一种方法是序列不要突然截短,而是缓慢截短,这就要求改变窗函数w(n)的形状,选择其他窗函数,同时加大N,使频谱W(ejω)主瓣能量集中,旁瓣能量更小,衰减更大,降低谱间干扰。典型窗函数有三角形窗、升余弦窗等,关于窗函数的详细内容将在第7章进行介绍。 3) 栅栏效应 针对离散傅里叶变换步骤,连续时间信号和序列都存在这种现象。由于N点DFT是频谱X(ejω)在频率[0,2π)上的等间隔采样,DFT就像一个“栅栏”,只能在离散的频率点上看到谱线,其他频率点的频谱看不到,这种现象称为“栅栏效应”。减轻栅栏效应的思路是增加频域采样点数,使离散谱线更密,就可以看到原来看不到的频谱分量,这些分量并不一定为零。具体做法可以采取序列尾部补零的方式,进行更大点数的DFT来实现。 3.5离散傅里叶变换的MATLAB仿真 本章主要从讨论周期序列的离散傅里叶级数出发,引出有限长序列的离散傅里叶变换,介绍了DFT的定义、基本性质和定理; 讨论了频域采样定理和频域内插公式; 着重探讨了利用DFT的两类典型应用,即计算线性卷积和谱分析。本节将通过MATLAB编程,对前面各节涉及的DFT定义、性质、频域采样定理以及应用示例等进行仿真和对比验证,以便更好地理解掌握相关概念和理论知识。 3.5.1DFT与IDFT计算仿真 1. DFT与IDFT函数的MATLAB编程 (1) 按照DFT定义,利用for循环编写函数DFT_For( ),对应m文件为DFT_For.m; function X=DFT_For(x,N) %利用for循环计算N点DFT; %x为输入序列,列向量; %X为输出列向量,X(k)=DFT[x(n)]; M=length(x); %序列原始长度,M<=N X=zeros(N,1); for k=0:N-1 for n=0:M-1 X(k+1)=X(k+1)+exp(-j*2*pi*n*k/N)*x(n+1); end end (2) 按照DFT矩阵表示形式编写函数DFT_Mat( ),对应m文件为DFT_Mat.m; function X=DFT_Mat(x,N) %利用矩阵形式计算N点DFT; %x为输入序列,列向量; %X为输出DFT,列向量,X(k)=DFT[x(n)]; M=length(x); %序列原始长度,M<=N x=[x;zeros(N-M,1)]; %补零 n=0:N-1;k=0:N-1; %变量范围 kn=k'*n; %生成k*n矩阵 WN=exp(-j*2*pi/N); %复指数WN WNkn=WN.^kn; %生成(WN).^kn X=WNkn*x; %矩阵相乘计算DFT (3) 类似地,可以按照for循环或矩阵形式编写IDFT函数。下面给出函数IDFT_Mat( )作为参考,对应m文件为IDFT_Mat.m; function x=IDFT_Mat(X,N) %利用矩阵形式计算N点IDFT; %X为输入DFT,列向量; %x为输出序列,列向量,x(n)=IDFT[X(k)]; M=length(X); %X(k)原始长度,M<=N X=[X;zeros(N-M,1)]; %补零 n=0:N-1;k=0:N-1; %变量范围 nk=n'*k; %生成n*k矩阵 WN=exp(-j*2*pi/N); %复指数WN WNnk=WN.^nk; %生成(WN).^nk x=1/N*conj(WNnk)*X; %矩阵形式计算IDFT (4) 除自己编写函数外,也可直接调用MATLAB函数fft( )和ifft( )。DFT和IDFT存在快速计算方法,即快速傅里叶变换和快速傅里叶逆变换,相关内容在第4章介绍。函数调用格式: X=fft(x,N); %x为序列,N为FFT点数 x=ifft(X,N); %一般N大于等于序列长度 2. 例3.2.1 计算DFT的MATLAB仿真 为便于对比DFT和DTFT,这里给出DTFT函数DTFT_Mat( ),对应m文件为DTFT_Mat.m; function [Xw,w]=DTFT_Mat(x,I) %利用矩阵形式计算序列DTFT; %x为输入序列,列向量,长度为N; %Xw为输出DTFT,w为频率,均为列向量,X(e^jw)=DTFT[x(n)]; %I为频率插值倍数,w和Xw长度均为N*I。 N=length(x);n=0:N-1; %序列长度和n范围 w=(0:N*I-1)'*2*pi/(N*I); %w范围0~2pi,等间隔取N*I个值 wn=exp(-j*w*n); %复指数exp(-jwn) Xw=wn*x; %矩阵相乘计算DTFT %主程序Ch3_5_1.m clc; clear all; close all; xn=[1 1 1 1]'; N=length(xn); %序列及长度 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %1)调用DTFT_Mat函数计算DTFT,并与MATLAB函数freqz()进行对比 [Xw,w]=DTFT_Mat(xn,100); %w范围0~2pi XwAmp=abs(Xw); %DTFT幅度 XwAng=angle(Xw)*180/pi; %DTFT相位,度数 figure(1); subplot(2,1,1);plot(w/pi,XwAmp); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|'); title('DTFT幅度');grid on subplot(2,1,2);plot(w/pi,XwAng); xlabel('\omega /\pi');ylabel('arg[X(e ^j^\omega)]'); title('DTFT相位');grid on %直接调用MATLAB函数freqz()作为对比 [Xw2,w2]=freqz(xn,1,N*100,'whole'); %频率包含整个周期0~2pi,取N*100个值 [Xw-Xw2 w-w2]; %观察可见: 两种方法的频率取值和DTFT相同 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2) 计算N=4,8,16,32点DFT N=4;Xk4=DFT_Mat(xn,N); %调用DFT_ Mat函数( ) N=8;Xk8=DFT_Mat(xn,N); N=16;Xk16=DFT_Mat(xn,N); N=32;Xk32=DFT_Mat(xn,N); % N=4;Xk4=DFT_For(xn,N); %或调用DFT_For函数( ) % N=8;Xk8=DFT_For(xn,N); % N=16;Xk16=DFT_For(xn,N); % N=32;Xk32=DFT_For(xn,N); figure(2); subplot(2,1,1);stem(0:4-1,abs(Xk4)); xlabel('k');ylabel('|X(k)|'); title('N=4点DFT'); axis([0 4 0 4]);grid on;set(gca,'XTick',0:4); subplot(2,1,2);stem(0:8-1,abs(Xk8)); xlabel('k');ylabel('|X(k)|'); title('N=8点DFT'); axis([0 8 0 4]);grid on;set(gca,'XTick',0:8); figure(3); subplot(2,1,1);stem(0:16-1,abs(Xk16)); xlabel('k');ylabel('|X(k)|'); title('N=16点DFT'); axis([0 16 0 4]);grid on;set(gca,'XTick',0:2:16); subplot(2,1,2);stem(0:32-1,abs(Xk32)); xlabel('k');ylabel('|X(k)|'); title('N=32点DFT'); axis([0 32 0 4]);grid on;set(gca,'XTick',0:4:32); 主程序Ch3_5_1.m运行结果如图3.5.1~图3.5.3所示,与2.2节的图以及图3.2.1一致。 (1) 图3.5.1给出了频率0~2π间的DTFT幅度和相位特性曲线,幅度包含2个主瓣和2个旁瓣,且关于ω=π左右对称,这是由于R4(n)是实数序列,DTFT具有共轭对称特性。相位具有线性特点,这是由于R4(n)时域上是偶对称的,满足线性相位条件,相关内容将在第7章中介绍。 图3.5.1编程计算矩形序列R4(n)的DTFT幅度和相位特性 (2) 观察图3.5.2和图3.5.3的DFT幅度特性,并对比图3.5.1中DTFT幅度曲线可以看出,DFT相当于DTFT在数字频率0~2π间的等间隔采样,像一个“栅栏”,中间频谱无法看到,并不一定是零,这种现象就是栅栏效应。当DFT点数加大时,中间频谱逐步呈现,DFT幅度包络愈发明显,与DTFT幅度曲线一致。 图3.5.2矩形序列R4(n)的4点和8点DFT幅度特性 图3.5.3矩形序列R4(n)的16点和32点DFT幅度特性 3.5.2DFT性质与定理仿真 1. 循环移位与循环卷积函数的MATLAB编程 1) 循环移位函数: y(n)=x((n+m))NRN(n) function y=cirshift(x,m,N) %输入序列为行向量,长度不超过N x=[x zeros(N-length(x))]; %补零,使长度为N n=0:N-1; y=x(mod(n+m,N)+1); %计算x2((n+m))N 2) 循环卷积函数: y(n)=x1(n) x2(n)=∑N-1m=0x1(m)x2((n-m))NRN(n) function y=circonv(x1,x2,N) %输入序列x1,x2为行向量,长度均不超过N x1=[x1 zeros(1,N-length(x1))];%补零,使长度为N x2=[x2 zeros(1,N-length(x2))]; %补零,使长度为N m=0:N-1; %x2flip=x2(mod(-m,N)+1); %计算x2((-m))N x2cir=zeros(N,N); %存储x2循环翻转 for n=0:N-1 x2cir(n+1,:)=x2(mod(n-m,N)+1); %计算x2((n-m))N end y=x1*conj(x2cir'); %序列乘累加 %y=x1*transpose(x2cir); %或使用转置函数transpose 3) 例3.2.3计算循环卷积的MATLAB仿真 %主程序Ch3_5_2_1.m clc;clear all x1=[1 1 1 1 0 0 0 0]; x2=[0 0 1 1 1 1 0 0]; N=8;n=0:N-1; y=circonv(x1,x2,N); %调用循环卷积函数 figure(1); subplot(2,2,1);stem(n,x1); xlabel('n');ylabel('x_1(n)'); title('x_1(n)');grid on subplot(2,2,2);stem(n,x2); xlabel('n');ylabel('x_2(n)'); title('x_2(n)');grid on subplot(2,1,2);stem(n,y); xlabel('n');ylabel('y(n)'); title('x_1(n)与x_2(n)的8点循环卷积');grid on 主程序Ch3_5_2_1.m运行结果如图3.5.4所示,与图3.2.3对比可见,循环卷积计算与仿真结果一致。 图3.5.4序列循环卷积仿真结果 2. 时域和频域循环卷积定理的MATLAB仿真 %主程序Ch3_5_2_2.m: clc;clear all x1=[1 1 1 1 0 0 0 0]; x2=[0 0 1 1 1 1 0 0]; N=8;n=0:N-1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %时域循环卷积定理验证: %先时域循环卷积,后DFT y=circonv(x1,x2,N); %调用循环卷积函数 Y=fft(y,N) %或DFT_Mat(conj(y'),N) %先DFT,后相乘 X1=fft(x1,N); X2=fft(x2,N); Y2=X1.*X2 %DFT[x1(n)].*DFT[x2(n)] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %频域循环卷积定理验证: %先相乘,后DFT x=x1.*x2 %x1(n).*x2(n) X=fft(x,N) %先频域循环卷积,后/N Xcir=circonv(X1,X2,N)/N %调用循环卷积函数 主程序Ch3_5_2_2.m运行后,Y和Y2输出结果一致,X和Xcir输出结果一致,验证了时域循环卷积定理和频域循环卷积定理。 Y = 16.0000-4.8284 + 4.8284i00.8284 + 0.8284i 00.8284 - 0.8284i0-4.8284 - 4.8284i x = 0 0 1 1 0 0 0 0 X = 2.0000-0.7071 - 1.7071i-1.0000 + 1.0000i0.7071 + 0.2929i 00.7071 - 0.2929i-1.0000 - 1.0000i-0.7071 + 1.7071i 3. 共轭对称性的MATLAB仿真 1) 序列实虚部DFT的共轭对称性 %主程序Ch3_5_2_3.m: clc;clear all N=8; %序列实虚部DFT的共轭对称性 xr=[4 3 2 1 0 3 2 1]; %序列实部 xi=[4 3 2 1 0 3 2 1]; %序列虚部 x=xr+j*xi; %实部+j*虚部 Xr=fft(xr,N) %DFT[实部]-共轭对称 Xij=fft(j*xi,N) %DFT[j*虚部]-共轭反对称 X=fft(x,N) %DFT[x(n)] 主程序Ch3_5_2_3.m运行后,输出结果如下: Xr = 16.00004.00000 - 4.0000i4.0000 04.00000 + 4.0000i4.0000 Xij = 0+16.0000i0 + 4.0000i4.00000 + 4.0000i 00 + 4.0000i-4.00000 + 4.0000i X = 16.0000+16.0000i4.0000+4.0000i4.0000-4.0000i4.0000+4.0000i 04.0000+4.0000i-4.0000+4.0000i4.0000+4.0000i 可以看出: Xr作为序列实部的DFT,满足Xr(k)=Xr*(N-k),具有共轭对称特性; Xij作为j乘序列虚部后的DFT,满足Xij(k)=-Xij*(N-k),具有共轭反对称特性; X=Xr+Xij表明共轭对称和共轭反对称两部分相加正好为构成序列的DFT。 2) 序列DFT实虚部的共轭对称性 %主程序Ch3_5_2_4.m: clc;clear all N=8; %序列DFT实虚部的共轭对称性 xep=[43+3j 2+2j 1+j01-j2-2j3-3j];%共轭对称序列 xop=[4j 3+3j 2+2j 1+j0-1+j-2+2j-3+3j]; %共轭反对称序列 x=xep+xop;%合成序列 Xep=fft(xep,N) %DFT[共轭对称序列]-DFT[x]实部 Xop=fft(xop,N) %DFT[共轭反对称序列]-j*DFT[x]虚部 X=fft(x,N) %DFT 主程序Ch3_5_2_4.m运行后,输出结果如下: Xep = 16.000016.48534.00002.8284 0-0.4853-4.0000-2.8284 Xop = 0 +16.0000i0 - 2.8284i0 - 4.0000i0 - 0.4853i 00 + 2.8284i0 + 4.0000i0 + 16.4853i X = 16.0000+16.0000i16.4853-2.8284i4.0000-4.0000i2.8284-0.4853i 0-0.4853+2.8284i-4.0000+4.0000i-2.8284+16.4853i 可以看出: Xep作为序列共轭对称部分的DFT,为实数; Xop作为共轭反对称部分的DFT,为纯虚数; Xep+Xop=X,两者分别对应序列DFT的实部、j乘虚部。 3) 共轭对称实序列、共轭反对称虚数序列DFT的对称性 %主程序Ch3_5_2_5.m: clc;clear all N=8; xr=[4 3 2 1 0 1 2 3];%序列实部 - 且满足共轭对称 xij=[4 3 2 1 0 1 2 3]*j; %序列虚部*j -且满足共轭反对称 x=xr+xij; %实部+j*虚部 Xr=fft(xr,N) %DFT[实部]-共轭对称,且为实数 Xij=fft(xij,N) %DFT[j*虚部]-共轭反对称,且为虚数 X=fft(x,N) %DFT[x] 主程序Ch3_5_2_5.m运行后,输出结果如下: Xr = 16.00006.828401.1716 0 1.1716 0 6.8284 Xij = 0+16.0000i0 + 6.8284i0 0 + 1.1716i 0 0 + 1.1716i0 0 + 6.8284i X = 16.0000+16.0000i6.8284+6.8284i01.1716+1.1716i 01.1716+1.1716i06.8284+6.8284i 可以看出: 由于xr作为序列实部,且满足共轭对称,因而Xr具有共轭对称特性,也为实数; xij为纯虚数,且满足共轭反对称,故Xij具有共轭反对称特性,也为纯虚数; 该仿真同时体现了序列实虚部DFT 和DFT实虚部的共轭对称性。 需要说明的是: 上述程序基于N=8仿真验证了偶数点序列DFT的共轭对称性; 类似地,构造奇数点序列进行仿真,如N=7,剔除N=8点序列中第4点,也可以得到关于共轭对称性的相同结论。 3.5.3频域采样与内插仿真 本节以矩形序列为例,仿真不同采样点数下的频域采样及内插恢复过程,探讨频域采样无失真条件,从而验证与3.3节理论分析的一致性。 %主程序Ch3_5_3.m: clc;clear all xn=[1 1 1 1 1 1]'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %1)调用DTFT_Mat函数计算DTFT M=length(xn);n=0:M-1;%序列长度和n范围 [Xw,w]=DTFT_Mat(xn,100); %w范围0~2pi,插值M*100倍 figure(1); subplot(2,1,1);stem(n,xn);grid on; xlabel('n');ylabel('x(n)'); title('原始序列x(n),长度为6'); axis([0 M 0 1]);grid on;set(gca,'XTick',0:M); subplot(2,1,2);plot(w/pi,abs(Xw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|');title('DTFT幅度'); axis([0 2 0 6]);grid on;set(gca,'XTick',0:2/M:2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %2)频域采样得到X(k) -> IDFT产生XN(n) -> 恢复DTFT并比较 nw=length(w); %频率值个数 N=4; %a)频域采样点数 < 原始序列长度 wstep=round(nw/N); %采样间隔 Xk=Xw(1:wstep:end); xNn4=IDFT_Mat(Xk,N); %时域恢复序列 [XNw wN]=DTFT_Mat(xNn4,100); %恢复频谱 figure(2); subplot(2,1,1);stem(0:N-1,abs(Xk));grid on xlabel('k');ylabel('|X(k)|'); title('N=4点频域采样X(k)幅度'); axis([0 N 0 6]);grid on;set(gca,'XTick',0:N); subplot(2,1,2);plot(wN/pi,abs(XNw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|');title('DTFT幅度'); axis([0 2 0 6]);grid on;set(gca,'XTick',0:2/M:2); N=6; %b)频域采样点数 = 原始序列长度 wstep=round(nw/N); %采样间隔 Xk=Xw(1:wstep:end); xNn6=IDFT_Mat(Xk,N); %时域恢复序列 [XNw wN]=DTFT_Mat(xNn6,100); %恢复频谱 figure(3); subplot(2,1,1);stem(0:N-1,abs(Xk));grid on xlabel('k');ylabel('|X(k)|'); title('N=6点频域采样X(k)幅度'); axis([0 N 0 6]);grid on;set(gca,'XTick',0:N); subplot(2,1,2);plot(wN/pi,abs(XNw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|');title('DTFT幅度'); axis([0 2 0 6]);grid on;set(gca,'XTick',0:2/M:2); N=8; %c)频域采样点数 > 原始序列长度 wstep=round(nw/N); %采样间隔 Xk=Xw(1:wstep:end); xNn8=IDFT_Mat(Xk,N); %时域恢复序列 [XNw wN]=DTFT_Mat(xNn8,100); %恢复频谱 figure(4); subplot(2,1,1);stem(0:N-1,abs(Xk));grid on xlabel('k');ylabel('|X(k)|'); title('N=8点频域采样X(k)幅度'); axis([0 N 0 6]);grid on;set(gca,'XTick',0:N); subplot(2,1,2);plot(wN/pi,abs(XNw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|');title('DTFT幅度'); axis([0 2 0 6]);grid on;set(gca,'XTick',0:2/M:2); 主程序Ch3_5_3.m运行结果显示如下所示: xNn4 =xNn6 =xNn8 = 2.00001.0000 - 0.0000i1.0000 - 0.0000i 2.0000 + 0.0000i1.0000 - 0.0000i1.0000 - 0.0000i 1.0000 - 0.0000i1.0000 - 0.0000i1.0000 - 0.0000i 1.0000 - 0.0000i1.0000 + 0.0000i1.0000 - 0.0000i 1.0000 + 0.0000i1.0000 - 0.0000i 1.0000 - 0.0000i1.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i DTFT频谱如图3.5.5~图3.5.8所示。 (1) 对比时域恢复序列,xNn6和原始序列xn一致,xNn8相当于xn补2个零,±0.0000是计算精度造成的,小于10-15,可以忽略。这与频域采样定理的序列关系式(3.3.7)一致。 图3.5.5矩形序列R6(n)及其DTFT幅度特性 图3.5.6矩形序列R6(n)DTFT经过4点频域采样及内插恢复情况 图3.5.7矩形序列R6(n)DTFT经过6点频域采样及内插恢复情况 图3.5.8矩形序列R6(n)DTFT经过6点频域采样及内插恢复情况 (2) 观察频谱恢复情况可见: 当频域采样点数大于或等于原始序列长度时,即N≥M,通过频域采样值X(k)能够无失真恢复出原始序列x(n)及其DTFT频谱。这与理论分析是一致的。 3.5.4DFT计算线性卷积仿真 1. 线性卷积与循环卷积的等价条件 下面以两个序列为例,通过仿真探讨线性卷积与循环卷积的等价条件。 %主程序Ch3_5_4_1.m: 针对3.4节中图3.4.1 clc;clear all hn=[1 1 1 1];%N=4 xn=[1 1 1 1 1]; %M=5 N=length(hn);M=length(xn);%原始序列长度 yn=conv(hn,xn) %调用线性卷积函数 L=N+M-1;n=0:L-1; %卷积长度: N+M-1 yn6=circonv(hn,xn,6);%调用循环卷积函数 yn8=circonv(hn,xn,8);%调用循环卷积函数 yn10=circonv(hn,xn,10);%调用循环卷积函数 figure(1); subplot(2,2,1);stem(n,yn); xlabel('n');ylabel('y(n)'); title('线性卷积,长度N+M-1');grid on subplot(2,2,2);stem(0:6-1,yn6); xlabel('n');ylabel('y(n)'); title('循环卷积,L=6');grid on subplot(2,2,3);stem(0:8-1,yn8); xlabel('n');ylabel('y(n)'); title('循环卷积,L=8');grid on subplot(2,2,4);stem(0:10-1,yn10); xlabel('n');ylabel('y(n)'); title('循环卷积,L=10');grid on 主程序Ch3_5_4_1.m运行结果如图3.5.9所示,与图3.4.1对比可以看出,序列线性卷积、循环卷积的计算与仿真结果一致; 并且当循环卷积长度大于或等于线性卷积长度时,循环卷积与线性卷积等价。 图3.5.9循环卷积与线性卷积等价条件仿真 2. 线性卷积的DFT计算方法 下面以两个序列线性卷积为例,对直接计算和利用DFT方法计算进行仿真对比。 %主程序Ch3_5_4_2.m: 针对图3.4.2 DFT计算线性卷积 clc;clear all hn=[1 1 1 1];%N=4 xn=[1 1 1 1 1]; %M=5 N=length(hn);M=length(xn);%原始序列长度 yn=conv(hn,xn) %调用线性卷积函数 L=N+M-1;n=0:L-1; %卷积长度: N+M-1 hn1=[hn zeros(1,L-N)]; %补零,使长度为L xn1=[xn zeros(1,L-M)]; %补零,使长度为L %先DFT,后相乘 Hk=fft(hn1,L); Xk=fft(xn1,L); Yk=Hk.*Xk; %DFT[h(n)].*DFT[x(n)] %IDFT yn2=ifft(Yk,L) 主程序Ch3_5_4_2.m运行结果如下: yn = 1 2 3 4 4 3 2 1 yn2 = 1 2 3 4 4 3 2 1 可以看出,采用DFT方法的计算结果yn2与线性卷积结果yn完全一致。这表明,在满足线性卷积与循环卷积等价条件时,采用DFT方法计算线性卷积是可行的。 3.5.5DFT谱分析仿真 1. 连续时间信号谱分析的频谱泄漏仿真研究 下面以不同频率间隔的双音信号为例,仿真呈现谱分析过程中的频谱泄漏现象,并探讨提高频率分辨率、减小谱间干扰的措施。 %主程序Ch3_5_5_1.m: clc;clear all; close all; %1)产生双音序列 f1=800;f2=1200;fs=5000;T=1/fs; %设置两个频率,满足时域采样定理 N=16;n=0:N-1; xn1=cos(2*pi*n*f1*T)'+cos(2*pi*n*f2*T)'; [Xw,w]=DTFT_Mat(xn1,N); %调用DFT_Mat函数( ) figure(1); subplot(2,1,1);stem(0:N-1,xn1); xlabel('n');ylabel('x(n)');axis([0 N -2 2]); title('800/1200Hz双音序列x(n),长度为16');grid on subplot(2,1,2);plot(w/pi,abs(Xw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|'); title('DTFT幅度');grid on %2)产生双音序列,两个频率更近 f1=1000;f2=1200;fs=5000;T=1/fs; %设置两个频率,满足时域采样定理 N=16;n=0:N-1; xn2=cos(2*pi*n*f1*T)'+cos(2*pi*n*f2*T)'; [Xw,w]=DTFT_Mat(xn2,N); %调用DFT_Mat函数( ) figure(2); subplot(2,1,1);stem(0:N-1,xn2); xlabel('n');ylabel('x(n)');axis([0 N -2 2]); title('1000/1200Hz双音序列x(n),长度为16');grid on subplot(2,1,2);plot(w/pi,abs(Xw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|'); title('DTFT幅度');grid on %3)增加序列长度,提高频率分辨率 N=32;n=0:N-1; xn3=cos(2*pi*n*f1*T)'+cos(2*pi*n*f2*T)'; [Xw,w]=DTFT_Mat(xn3,N); %调用DFT_Mat函数( ) figure(3); subplot(2,1,1);plot(w/pi,abs(Xw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|'); title('32点1000/1200Hz双音序列的DTFT幅度');grid on %4)加窗+增加序列长度,进一步减小谱间干扰 N=64;n=0:N-1; xn4=cos(2*pi*n*f1*T)'+cos(2*pi*n*f2*T)'; Win=hanning(N); xNn=xn4.*Win; [Xw,w]=DTFT_Mat(xNn,N); %调用DFT_Mat函数() subplot(2,1,2);stem(0:N-1,Win); xlabel('n');ylabel('w(n)'); axis([0 N 0 1]); title('窗函数w(n),长度为64');grid on figure(4); subplot(2,1,1);stem(0:N-1,xNn); xlabel('n');ylabel('x_N(n)'); title('加窗后序列x(n)*w(n),长度为64');grid on subplot(2,1,2);plot(w/pi,abs(Xw)); xlabel('\omega /\pi');ylabel('|X(e ^j^\omega)|'); title('DTFT幅度');grid on 主程序Ch3_5_5_1.m运行结果如图3.5.10~图3.5.13所示。 对比图3.5.10和图3.5.11可以看出: 尽管单音连续时间信号的傅里叶变换理论上为单位冲激样式的频谱,但经过时域采样和截短后,出现了展宽的主瓣和拖尾的旁瓣,形成频谱泄漏现象。主瓣宽度直接影响了频率分辨率,当双音频率为800Hz、1200Hz时,频谱上呈现双尖峰特性,还可以分辨出两个频率; 但当双音频率靠近为1000Hz、1200Hz时,主瓣部分叠加,呈现单峰特性,无法有效分辨两个频率。 图3.5.10双音序列时域与频域分析结果 图3.5.11双音频率靠近的序列时域与频域分析结果 为了提高频率分辨率,考虑增加时域截短长度,即增加矩形窗长度,从而降低窗函数主瓣宽度。图3.5.12给出了长度增加1倍后的DTFT幅度分析结果。可以看出,随着时域截短序列长度增加,双音频率的主瓣宽度变窄,峰值特性更加尖锐,能够有效分辨两个频率。同时也注意到旁瓣幅度相对较高,谱间干扰比较明显。 图3.5.12增加序列长度改善频率分辨率、窗函数形状仿真 为了降低谱间干扰,同时保证频率分辨率,考虑在增加窗函数长度的基础上,选择更加平缓的窗函数形状。图3.5.12给出了64阶窗函数w(n)样点值,图3.5.13则给出了加窗后序列样点值及谱分析结果。观察DTFT幅度可见,双音频率可以分辨,并且显著降低了谱间干扰。 图3.5.13窗函数改变形状、增加长度降低谱间干扰仿真 2. 周期序列的谱分析仿真 下面针对周期序列,取不同周期数构成长度不同的序列,仿真研究了序列谱分析结果的相互关系。 %主程序Ch3_5_5_2.m: clc;clear all;close all; xn=[0 1 2 3 4 5 6 7]'; %单个周期 N=length(xn); Xk8=DFT_Mat(xn,N); %调用DFT_Mat函数( ) Xk16=DFT_Mat([xn;xn],2*N);%重复,2个周期 Xk24=DFT_Mat([xn;xn;xn],3*N); %重复,3个周期 figure(1); subplot(2,1,1);stem(0:N-1,xn); xlabel('n');ylabel('x(n)'); title('x(n) 单个周期'); axis([0 N 0 8]);grid on;set(gca,'XTick',0:N); subplot(2,1,2);stem(0:N-1,abs(Xk8)); xlabel('k');ylabel('|X(k)|'); title('N=8点DFT'); axis([0 N 0 30]);grid on;set(gca,'XTick',0:N); figure(2); subplot(2,1,1);stem(0:2*N-1,[xn;xn]); xlabel('n');ylabel('x(n)'); title('x(n) 两个周期'); axis([0 2*N 0 8]);grid on;set(gca,'XTick',0:2:2*N); subplot(2,1,2);stem(0:2*N-1,abs(Xk16)); xlabel('k');ylabel('|X(k)|'); title('N=16点DFT'); axis([0 2*N 0 60]);grid on;set(gca,'XTick',0:2:2*N); figure(3); subplot(2,1,1);stem(0:3*N-1,[xn;xn;xn]); xlabel('n');ylabel('x(n)'); title('x(n) 三个周期'); axis([0 3*N 0 8]);grid on;set(gca,'XTick',0:3:3*N); subplot(2,1,2);stem(0:3*N-1,abs(Xk24)); xlabel('k');ylabel('|X(k)|'); title('N=24点DFT'); axis([0 3*N 0 90]);grid on;set(gca,'XTick',0:3:3*N); set(gca,'YTick',0:30:90); 主程序Ch3_5_5_2.m运行结果如图3.5.14~图3.5.16所示。对比DFT幅度可以看出,对周期序列进行谱分析时,取m个周期序列DFT相当于单个周期序列DFT间隔插入m-1个0,且幅度扩大到m倍,m=2,3,…。这与3.4.2节中理论分析是一致的。 图3.5.14取单个周期的序列DFT谱分析 图3.5.15取两个周期的序列DFT谱分析 图3.5.16取三个周期的序列DFT谱分析 习题 3.1设x(n)=R4(n),x~(n)=x((n))6,试求X~(k),并作图表示x~(n)、X~(k)。 3.2已知周期序列x~(n),X~(k)=DFS[x~(n)],试求下列周期序列的DFS: (1) x~(n+m)(2) x~*(n) (3) x~*(-n)(4) Re[x~(n)] 3.3求下列有限长序列的N点离散傅里叶变换: (1) δ(n)(2) δ(n-n0),0<n0<N (3) x(n)=anRN(n)(4) x(n)=Rm(n),0<m<N (5) x(n)=ejω0nRN(n)(6) x(n)=ej2πNmnRN(n),0<m<N (7) x(n)=cos(ω0n)RN(n)(8) x(n)=nRN(n) 3.4已知下列X(k): (1) X(k)=N2ejθ,k=m N2e-jθ,k=N-m 0,其他(2) X(k)=-N2jejθ,k=m N2je-jθ,k=N-m 0,其他 求x(n)=IDFT[X(k)],其中m为正整数,0<m<N/2。 3.5已知序列x(n)=4δ(n)+3δ(n-1)+2δ(n-2)+δ(n-3),画出序列x1(n)和x2(n)的图形: (1) x1(n)=x((n-2))4R4(n)(2) x1(n)=x((2-n))4R4(n) 3.6已知长度N=10的两个有限长序列: x1(n)=1,0≤n≤4 0,5≤n≤9x2(n)=1,0≤n≤4 -1,5≤n≤9 试作图表示x1(n)、x2(n)以及两者循环卷积x(n)=x1(n) x2(n)。 3.7已知两个有限长序列为 x(n)=n+1,0≤n≤3 0,4≤n≤6y(n)=-1,0≤n≤4 1,5≤n≤6 试作图表示x(n)、y(n)以及N=8点循环卷积w(n)=x(n) y(n)。 3.8已知4点序列x(n)=δ(n)+2δ(n-1)+2δ(n-2)+δ(n-3) ; (1) 画出线性卷积x(n)*x(n)的图形 (2) 画出N=4点循环卷积x(n) x(n)的图形 (3) 画出N=7点循环卷积x(n) x(n)的图形,并与线性卷积结果比较,说明线性卷积与循环卷积的关系。 3.9证明DFT的频域循环卷积定理。 3.10如果X(k)=DFT[x(n)],证明DFT的对称定理: DFT[X(n)]=Nx(N-k) 3.11如果X(k)=DFT[x(n)],证明DFT的初值定理: x(0)=1N∑N-1k=0X(k) 3.12证明离散相关定理。若X(k)=X*1(k)X2(k),则有 x(n)=IDFT[X(k)]=∑N-1l=0x1*(l)x2((n+l))NRN(n) 3.13已知N点有限长序列x(n),X(k)=DFT[x(n)],将x(n)补零后形成长度为rN的序列: y(n)=x(n),0≤n≤N-1 0,N≤n≤rN-1 试求y(n)的rN点离散傅里叶变换(Y(k)=DFT[y(n)])与X(k)的关系。 3.14已知N点有限长序列x(n),X(k)=DFT[x(n)],将x(n)每点后插入r-1个零,得到长度为rN的有限长序列: y(n)=x(n/r),n=ir,0≤i≤N-1 0,其他 试求y(n)的rN点DFT与X(k)的关系。 3.15已知序列x(n)=anu(n),0<a<1,对x(n)的Z变换X(z)在单位圆上等间隔采样N点,采样值为 X(k)=X(z)|z=W-kN,k=0,1,…,N-1 求有限长序列IDFT[X(k)]。 3.16用计算机对实序列进行谱分析,要求谱分辨率F≤50Hz,信号最高频率为3kHz,试确定下列参数: (1) 最小记录时间Tpmin; (2) 最大采样间隔Tmax; (3) 最少采样点数Nmin; (4) 若信号带宽不变,频率分辨率提高1倍时的N值。 3.17假设频谱分析时DFT点数必须为2的整数幂,要求分辨率小于或等于10Hz,如果采样时间间隔为0.1ms,试确定: (1) 最小记录长度; (2) 所允许处理的信号最高频率; (3) 在一个记录中的最少点数,此时实际分辨率为多少。 3.18在某数字通信系统,对模拟信号以9.6kHz速率进行采样,频谱分析时计算256点离散傅里叶变换,试确定频谱分辨率和采样数据记录时间。若计算512点DFT,频谱分辨率和记录时间又是多少。 3.19设连续时间信号x(t)=sin(2πfct)+sin(2πfJt),其中有用信号频率fc=1800Hz,干扰频率fJ=2000Hz: (1) 确定合适的采样频率及采样点数N,利用MATLAB编写DFT谱分析程序,并画出DFT幅度和相位特性。 (2) 探讨提高频率分辨率的措施,并编程验证。