第3章 CHAPTER 3 离散傅里叶变换 计算机只能计算有限长离散序列,因此有限长序列在数字信号处理中就显得很重要,虽然可以用Z变换和序列傅里叶变换研究它,但是这两种变换无法直接利用计算机进行数值计算。针对序列“有限长”这一特点,可以导出一种更有效的变换——离散傅里叶变换(Discrete Fourier Transform,DFT)。作为有限长序列的一种表示方法,离散傅里叶变换除了在理论上相当重要之外,由于存在有效的快速算法——快速傅里叶变换(Fast Fourier Transform,FFT),因而在各种数字信号处理的算法中起着核心作用。本章主要讨论离散傅里叶级数和离散傅里叶变换定义性质,圆周卷积,利用DFT计算线性卷积和以及频域采样理论等内容。 视频讲解 3.1引言 傅里叶变换是以时间t,n为自变量的“信号”与以频率(Ω、f或ω )为自变量的“频谱”函数之间的一种变换关系。当自变量“时间”与“频率”为连续形式和离散形式的不同组合时,就形成了各种不同形式的傅里叶变换对,即“信号”与“频谱”的对应关系。 1. 连续时间、非周期信号的傅里叶变换 连续时间、非周期信号通过连续傅里叶变换(FT)得到非周期连续频谱密度函数。 正变换 Xa(jΩ)= ∫∞-∞xa(t)e-jΩtdt(31) 逆变换 xa(t)=12π∫∞-∞Xa(jΩ)ejΩtdΩ(32) 从图31的矩形脉冲及其频谱可以看出: (1) 时域连续函数造成频域是非周期的谱; (2) 时域的非周期造成频域是连续的谱。 图31矩形脉冲及其频谱 2. 连续时间周期信号的傅里叶级数 周期为T的周期性连续时间函数xa(t)可展开成傅里叶级数Fn,是离散非周期性频谱。 正变换 Fn=1T∫T2-T2xa(t)e-jnΩtdt(33) 逆变换 xa(t)=∑∞n=-∞FnejnΩt(34) 由图32的周期矩形脉冲及其频谱可以看出: (1) 时域的连续函数造成频域是非周期的频谱函数; (2) 频域的离散频谱与时域的周期时间函数对应(频域采样,时域周期延拓)。 图32周期矩形脉冲及其频谱 3. 非周期离散信号的傅里叶变换 非周期离散的时间信号得到周期性连续的频率函数。 正变换 X(ejw)=∑∞n=-∞x(n)e-jωn(35) 逆变换 x(n)=12π∫π-πX(ejω)ejωndω(36) 由图33的非周期离散时间信号及其频谱可以看出 (1) 时域的离散造成频域的周期延拓; (2) 时域的非周期对应于频域的连续。 图33非周期离散时间信号及其频谱 4. 傅里叶变换的4种形式 通过上面分析,可以得出一般性规律: 一个域的离散对应另一个域的周期延拓,一个域的连续必定对应另一个域的非周期。上面的3种傅里叶变换对,都不适用在计算机上运算,因为至少在一个域(时域或频域)中函数是连续的。因为从数字计算角度,感兴趣的是时域及频域都是离散的情况,就是这里要谈到的离散傅里叶变换。 傅里叶变换的4种形式如下: (1) 连续时间、连续频率——傅里叶变换(FT),是非周期连续信号傅里叶变换,如图34(a)所示; (2) 连续时间、离散频率——傅里叶级数(FS); 是周期连续信号傅里叶级数,如图34(b)所示; (3) 离散时间、连续频率——序列傅里叶变换(DTFT); 是序列傅里叶变换,如图34(c)所示; (4) 离散时间、离散频率——离散傅里叶变换(DFT),是离散周期序列频谱,如图34(d)所示。 图34各种形式的傅里叶变换 视频讲解 3.2周期序列的离散傅里叶级数 3.2.1离散傅里叶级数定义 若离散时间序列x(n)为周期序列,则一定满足 x(n)=x(n+rN) 其中,N(正整数)为信号的周期,r为任意整数。为了与非周期序列区分,周期序列记作: x~(n)。 因为周期序列不是绝对可和,因此周期序列不能用傅里叶变换表示,但是周期序列可以用傅里叶级数(Discrete Fourier Series,DFS)表示,离散傅里叶级数定义为 x~(n)=1N∑N-1k=0X~(k)ej2πNnk=IDFS[X~(k)](37) 其中,X~(k)为周期序列傅里叶级数的系数,表示为 X~(k)=∑N-1n=0x~(n)e-j2πNnk=DFS[x~(n)](38) 为了书写方便,常令符号WN=e-j2πN,这样周期序列傅里叶变换对可以再次写为 正变换 X~(k)=DFS[x~(n)]=∑N-1n=0x~(n)WnkN 反变换 x~(n)=IDFS[X~(k)]=1N∑N-1k=0X~(k)W-nkN MATLAB程序实现DFS过程如下: function [Xk]=dfs(xn,N)%计算(DFS)系数 %[Xk]=dfs(xn,N) %Xk=在0<=n<=N-1之间的一个单周期信号 %N=xn的基本周期 n=[0:1:N-1]; %n的行向量 k=[0:1:N-1]; %k的行向量 WN=exp(-j*2*pi/N); %Wn因子 nk=n'*k; %产生一个含nk值的N×N维矩阵 Xk=xn* WN.^nk; %DFS 系数行向量 MATLAB程序实现IDFS过程如下: function [xn]=idfs(Xk,N) %计算逆(DFS)系数 %[xn]=dfs(Xk,N) % xn =周期信号在0<=n<=N-1之间的一个单周期信号 % Xk =在0<=k<=N-1间的DFS系数数组 %N=Xk的基本周期 n=[0:1:N-1]; %n的行向量 k=[0:1:N-1]; %k的行向量 WN=exp(-j*2*pi/N); %Wn因子 nk=n'*k; %产生一个含nk值的N×N维矩阵 xn=Xk* (WN.^(-nk))/N; %DFS 系数行向量 【例31】设x~(n)为周期脉冲串,求x~(n)=∑∞r=-∞δ(n+rN)。 【解】因为对于0≤n≤N-1,x~(n)=δ(n),x~(n)的DFS系数为 X~(k)=∑N-1n=0x~(n)WnkN=∑N-1n=0δ(n)WnkN=1 在这种情况下,对于所有k值,X~(k)均相同,于是有 x~(n)=∑∞r=-∞δ(n+rN)=1N∑N-1k=0W-nkN=1N∑N-1k=0ej2πNnk 当n为N整数倍时结果为1,这正好是周期性脉冲串。 说明: 1N∑N-1n=0ej2πNkn=1N1-ej2πNkN1-ej2πNk=1,k=mN,m为整数 0,其他(39) 式(39)一般称为复正弦序列正交特性。 【例32】已知周期序列x~(n)如图35所示,其周期N=10,求解其傅里叶级数系数X~(k)。 图35周期脉冲序列(N=10) 【解】 X~(k)=∑10-1n=0x~(n)Wnk10=∑4n=0e-j2π10nk 这一有限求和有闭合形式 X~(k)=1-W5k101-Wk10=e-j4πk10sin(5πk/10)sin(5πk/10) 周期脉冲序列DFS如图36所示, MATLAB程序实现DFS过程如下: >>xn=[1,1,1,1,1,0,0,0,0,0];N=10; >>Xk=dfs(xn,N) Xk = Columns 1 through 5 5.0000 + 0.0000i1.0000 - 3.0777i-0.0000 + 0.0000i1.0000 - 0.7265i -0.0000 + 0.0000i Columns 6 through 10 1.0000 - 0.0000i -0.0000 - 0.0000i1.0000 + 0.7265i-0.0000 - 0.0000i 1.0000 + 3.0777i 图36周期脉冲序列DFS 下面给出L=5,N=20周期方波的离散傅里叶级数的MATLAB程序: L=5;N=20;k=[-N/2:N/2]; %方波参数 xn=[ones(1,L),zeros(1,N-L)]; %方波x(n) Xk=dfs(xn,N); %DFS magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]); %DFS幅度 subplot(2,2,1);stem(k,magXk); axis([-N/2,N/2,-0.5,5.5]); xlabel('k');ylabel('Xtilde(k)'); title('方波的DFS:L=5,N=20') 上述程序产生的图及其他情形如图37所示,注意到X~(k)是周期信号,图中只画出了从-N/2到N/2的部分。从图37可以看出,方波的DFS系数包络看起来像sinc函数,k=0时的幅度为L,同时函数的零点位于N/L(占空比的倒数)的整数倍处。 图37不同L和N周期方波的离散傅里叶级数幅度 正变换定义公式(38)中的周期序列X~(k)可看成是对x~(n)的第一个周期x(n)作Z变换,然后将Z变换在Z平面单位圆上按等间隔角2π/N采样而 图38Z平面单位圆上等间隔采样 得到,如图38所示。令 x(n)= x~(n)·RN(n)=x~(n),0≤n≤N-1 0,其他(310) 通常称x(n)为x~(n)的主值区序列,则x(n)的Z变换为 X(z)=∑∞n=-∞x(n)z-n=∑N-1n=0x~(n)z-n(311) X~(k)=X(z)z=W-kN=ej 2πNk(312) 由于单位圆上的Z变换为序列傅里叶变换,周期序列X~(k)也可以解释为x~(n)的一个周期x(n)的傅里叶变换的等间隔采样。因为 X(ejω)=∑N-1n=0x(n)e-jωn=∑N-1n=0x~(n)e-jωn(313) 比较式(38)和式(313),可以看出 X~(k)= X(ejω)ω=2πk/N=∑N-1n=0x~(n)e-j2πNkn,k=0,1,2,…,N-1(314) 这相当于以2π/N的频率间隔对傅里叶变换进行采样。也就是说,非周期离散时间信号经序列傅里叶变换(DTFT),得周期连续谱函数,再经采样得周期离散频谱函数(DFS),过程如图39所示。 图39序列傅里叶变换与离散傅里叶级数关系 由于频域N点取样使得频域离散从而形成时域序列的周期化。采样频点间隔为 ω0=2πN(315) 数字频率为 ω=kω0=2πNk,k=0,1,2,…,N-1(316) 【例33】傅里叶级数系数X~(k)和周期信号x~(n)的一个周期的傅里叶变换之间的关系。 【解】图310为周期序列x~(n),它的一个周期为 x(n)=1,0≤n≤4 0,其他 x~(n)的一个周期的傅里叶变换为 X(ejω)=∑4n=0e-jωn=1-e-j5ω1-e-jω=e-j2ωsin(5ω/2)sin(ω/2) 图310周期序列 根据|x(ejω)|绘制一个周期的DTFT幅度谱如图311所示。 根据上面分析,傅里叶级数为 X~(k)=X(ejω)ω=2πk/10=e-j4πk10sin(5πk/10)sin(πk/10) 图311x~(n)一个周期的DTFT幅度谱 根据X~(k)绘制周期序列傅里叶级数DFS如图312所示,可以看出X~(k)为|X(ejω)|在[0,2π]的N点等间隔采样。 图312x~(n)一个周期的DTFT幅度谱和DFS 3.2.2离散傅里叶级数的性质 由于可以用采样Z变换解释DFS,因此它的许多性质与Z变换性质非常相似。但是,由于x~(n)和X~(k)两者都具有周期性,这就使它与Z变换性质还有一些重要的差别。此外,DFS在时域和频域之间具有严格的对偶关系,这是序列的Z变换表示所不具有的。设x~1(n)和x~2(n)皆是周期为N的周期序列,各自的DFS分别为X~1(k)、X~2(k)。 1. 线性 DFS[ax~1(n)+bx~2(n)]=aX~1(k)+bX~2(k)(317) 2. 序列的移位 DFS[x~(n+m)]=W-mkNX~(k)=ej2πNmkX~(k)(318) DFS[WnlNx~(n)]= X~(k+l)(319) 或 IDFS[X~(k+l)]=WnlNx~(n)=e-j2πNnlx~(n)(320) 3. 周期卷积 如果 Y~(k)= X~1(k)X~2(k)(321) 则 y~(n)=IDFS[Y~(k)]=∑N-1m=0x~1(m)x~2(n-m)(322) 或 y~(n)=∑N-1m=0x~2(m)x~1(n-m)(323) 两个周期都为N的周期序列x~1(n)和x~2(n),其卷积的结果也是周期为N的周期序列,求和只在一个周期上进行,即m为0到N-1,所以称为周期卷积。注意: n的取值不在N的范围内,得到结果后进行周期延拓,如图313所示。 图313两个周期序列(N=6)的周期卷积过程 视频讲解 3.3有限长序列离散傅里叶变换 3.3.1离散傅里叶变换定义 由于长度为N的有限长序列可以看作周期为N的周期序列的一个周期,利用DFS计算周期序列的一个周期,就可以得到有限长序列的离散傅里叶变换。设x(n)是长度为N的有限长序列,可以把它看作是周期为N的周期序列x(n)的一个主周期,而将x~(n)看作x(n)以N为周期进行周期延拓得到,即 x(n)=x~(n),0≤n≤N-1 0,其他 = x~(n)RN(n)(324) 同理 X~(k)=X((k))N(325) X(k)= X~(k)RN(k)(326) 通常,把x~(n)的第一个周期n=0到n=N-1定义为主值区间,故x(n)是x~(n)的主值序列。而称x~(n)为x(n)的周期延拓(图314)。 图314周期序列、主值区间、主值序列之间关系 DFS与IDFS变换前后都是周期的、无限长的,但这里一个周期信息与其他周期信息相同,因而可以得到有限长序列的离散傅里叶变换的定义: X(k)=DFT[x(n)]=∑N-1n=0x(n)WnkN,0≤k≤N-1 x(n)=IDFT[X(k)]=1N∑N-1k=0X(k)W-nkN,0≤n≤N-1(327) (328) 注意,所处理的有限长序列都是作为周期序列的一个周期表示的。换句话说,离散傅里叶变换隐含周期性。 【例34】已知序列x(n)=δ(n),求它的N点DFT。 【解】单位脉冲序列的DFT很容易由DFT的定义式(327)得到 X(k)=∑N-1n=0δ(n)WnkN=W0N=1,k=0,1,2,…,N-1 δ(n)的X(k)如图315所示。这是一个很特殊的例子,它表明对序列δ(n)来说,不论对它进行多少点的DFT,所得结果都是一个离散矩形序列。 图315δ(n)及其DFT 【例35】已知x(n)=cos(nπ/6)是一个长度N=12的有限长序列,求它的N点DFT。 【解】由DFT的定义式(327)可得 X(k)=∑11n=0cosnπ6Wnk12=∑11n=012ejnπ6+e-jnπ6e-j2π12nk =12∑11n=0e-j2π12n(k-1)+∑11n=0e-j2π12n(k+1) 利用复正弦序列正交特性式(39),则 ∑N-1n=0ej2πN(k-m)n=N,k=m 0,k≠m 再考虑k的取值区间,可得 X(k)=6,k=1,11 0,其他 有限长余弦序列及其DFT如图316所示。 图316有限长余弦序列及其DFT 【例36】已知X(k)=3,k=0 1,1≤k≤9,求k=0,1,2,…,9点IDFT。 【解】X(k)可以表示为X(k)=1+2δ(k),0≤k≤9。 由于一个单位脉冲序列的DFT为常数 x1(n)=δ(n)X1(k)=DFT[x1(n)]=1 同样,一个常数的DFT是一个单位脉冲序列 x2(n)=1X2(k)=DFT[x2(n)]=Nδ(k) X(k)=∑10-1n=0x(n)e-j2π10nk=1-e-j2πk1-e-j2π10k=10,k=0 0,k=1,2,…,9 =10δ(k) 由于 10δk→1 2δk→15 所以 x(n)=15+δ(n) 3.3.2DFT与序列傅里叶变换、Z变换的关系 若x(n)是一个有限长序列,长度为N,Z变换X(z)=∑N-1n=0x(n)z-n,则 X(z)z=W-kN=∑N-1n=0x(n)WnkN=DFT[x(n)]=X(k)(329) z=W-kN=ej2πNk表明W-kN是Z平面单位圆上幅角为ω=2πNk的点,即将Z平面单位圆N等分后的第k点,如图317所示。所以X(k)是对X(z)在Z平面单位圆上N点等间隔采样值。此外,由于序列的傅里叶变换X(ejω)即是单位圆上的Z变换,根据式(329),DFT与序列傅里叶变换的关系为 X(k)=X(ejω)ω=2πNk=X(ejkωN)(330) 其中,ωN=2πN。 DFT的物理意义: 说明X(k)可看作序列x(n)的傅里叶变换X(ejω)在区间[0, 2π]上的N点等间隔采样,其采样间隔为ωN=2π/N。显而易见,DFT的变换区间长度N不同,表示对X(ejω)在区间[0, 2π]上的采样间隔和采样点数不同,所以DFT的变换结果也不同。 图317DFT物理意义说明图示 【例37】有限长序列x(n)为 x(n)=1,0≤n≤4 0,其他 求其N=5点离散傅里叶变换X(k)。 【解】将有限长序列x(n),如图318(a)所示,以N=5为周期将x(n)延拓成周期序列x~(n),如图318(b),x~(n)的DFSX~(k)与x(n)的DFT X(k)相对应如图318(c)和图318(d)所示。因为图318(b)中的序列在区间0≤n≤N-1上为常数值,所以可以得出 图318例37图 X~(k)=∑N-1n=0e-j(2πk/N)n=1-e-j2πk1-e-j(2πk/N)=N,k=0, ±N, ±2N, … 0,其他 在k=0和k=N的整数倍处才有非零的DFS系数X~(k)值。就是X(ejω)在频率ωk=2πk/N处的样本序列。x(n)的DFT对应于取X~(k)的一个周期而得到的有限长序列X(k)。 X(k)=∑5-1n=0x(n)e-j2π5nk=1-e-j2πk1-e-j2π5k=5,k=0 0,k=1, 2, 3, 4 视频讲解 3.4离散傅里叶变换的性质 本节讨论DFT的一些性质,它们本质上与周期序列的DFS概念有关,而且是由有限长序列及其DFT表示式隐含的周期性得出的。以下讨论的序列都是N点有限长序列,用DFT[·]表示N点DFT,且设DFT[x1(n)]=X1(k)、DFT[x2(n)]=X2(k)。 1. 线性 若两个有限长序列x1(n)和x2(n)的线性组合x3(n)=ax1(n)+bx2(n),则有 DFT[ax1(n)+bx2(n)]=aX1(k)+bX2(k)(331) 其中,a,b为任意常数。 说明: (1) 若x1(n)和x2(n)的长度均为N,则x3(n)的长度为N; (2) 若x1(n)和x2(n)的长度不等,x1(n)的长度为N1,x2(n)的长度为N2,则x3(n)的长度为N=max[N1,N2],离散傅里叶变换的长度必须按N计算。 2. 圆周移位 一个长度为N的有限长序列x(n)的圆周移位定义为y(n)=x((n+m))NRN(n),将图319(a)所示x(n)以N为周期进行周期延拓,得到周期序列x~(n)=x((n))N,如图319(b)所示; 再将x~(n)加以移位x((n+m))N=x~(n+m)如图319(c)所示,然后,再对移位的周期序列x~(n+m)取主值区间(n为0到N-1)上的序列值,即x((n+m))NRN(n)。所以,一个有限长序列x(n)的圆周移位序列y(n)仍然是一个长度为N的有限长序列。 图319圆周移位过程 3. 时域圆周移位定理 设x(n)是长度为N的有限长序列,y(n)为x(n)圆周移位,即 y(n)=x((n+m))NRN(n)(332) 则圆周移位后的DFT为 Y(k)=DFT[y(n)]=DFT[x((n+m))NRN(n)]=W-mkNX(k)(333) 4. 频域圆周移位定理 若X(k)=DFT[x(n)],则 IDFT[X((k+l))NRN(k)]=WnlNx(n)=e-j2πNnlx(n)(334) 上式称为频率移位定理,也称为调制定理,此定理说明时域序列的调制等效于频域的圆周移位。 5. 圆周卷积 设x1(n)和x2(n)都是点数为N的有限长序列(0≤n≤N-1),且有 Y(k)=X1(k)X2(k) 则 y(n)=IDFT[Y(k)]=∑N-1m=0x1(m)x2((n-m))NRN(n) =∑N-1m=0x2(m)x1((n-m))NRN(n)(335) 卷积过程: 圆周卷积流程如图320所示,过程如图321所示。先将x1(n)和x2(n)补零,使得长度均为N点,并将变量n变成m,x2(m)周期化x2((m))N,再反转x2((-m))N,取主值序列x2((-m))NRN(m),对x2(m)圆周右移n,形成x2((n-m))NRN(m),当n=0,1,2,…,N-1时,分别将x1(m)与x2((n-m))NRN(m)相乘,并在m为0到N-1区间内求和,便得到圆周卷积y(n)。 图320圆周卷积流程 图321圆周卷积过程示意图 特别要注意: 两个长度小于或等于N的序列的N点圆周卷积长度仍为N,这与一般的线性卷积不同。为了区别线性卷积,用*表示线性卷积,用 表示N点圆周卷积。 y(n)=x1(n)x2(n)=∑N-1m=0x1(m)x2((n-m))NRN(n)(336) 利用时域与频域的对称性,可以证明频域圆周卷积定理。 若y(n)=x1(n)x2(n), x1(n),x2(n)皆为N点有限长序列,则 Y(k)=DFT[y(n)]=1N∑N-1l=0X1(l)X2((k-l))NRN(k) =1N∑N-1l=0X2(l)X1((k-l))NRN(k)=1NX1(k)X2(k)(337) 即时域序列相乘,乘积的DFT等于各个DFT的圆周卷积再乘以1N。 6. 线性卷积与圆周卷积关系 时域圆周卷积在频域上相当于两序列的DFT的乘积,而计算DFT可采用它的快速算法——快速傅里叶变换(FFT),因此圆周卷积与线性卷积相比,计算速度可以大大加快。若序列xi(n)的长度为Ni(i=1,2),当N1与N2大小相当时,它们的线性卷积可以用L点DFT快速实现,其中 L≥N1+N2-1,如图322所示。 图322用圆周卷积和计算出线性卷积和的过程 实际问题大多总是要求解线性卷积。例如,信号通过线性时不变系统,其输出就是输入信号与系统的单位脉冲响应的线性卷积,如果信号以及系统的单位脉冲响应都是有限长序列,那么是否能用圆周卷积运算来代替线性卷积运算而不失真呢?下面就来讨论这个问题。 设x1(n)是N1点的有限长序列(0≤n≤N1-1),x2(n)是N2点的有限长序列(0≤n≤N2-1)。 1) 线性卷积 y1(n)=x1(n)x2(n)=∑∞m=-∞x1(m)x2(n-m)=∑N1-1m=0x1(m)x2(n-m)(338) y1(n)是N1+N2-1点有限长序列,即线性卷积的长度等于参与卷积的两序列的长度之和减1。 2) 圆周卷积 先假设进行L点的圆周卷积,再讨论L取何值时,圆周卷积才能代表线性卷积。设y(n)=x1(n)x2(n)是两序列的L点圆周卷积,L≥max[N1, N2],这就要将x1(n)与x2(n)都看成是L点的序列。在这L点的序列值中,x1(n)只有前N1个是非零值,后L-N1个均为补充的零值。同样,x2(n)只有前N2个是非零值,后L-N2个均为补充的零值。则 y(n)=x1(n)x2(n)=∑L-1m=0x1(m)x2((n-m))LRL(n)(339) 可以证明 y~(n)=∑∞r=-∞y1(n+rL)(340) y~(n)为x1(n)与x2(n)线性卷积的周期延拓,周期也为L,定义为周期卷积。 前面已经分析了y1(n)具有N1+N2-1个非零值。因此可以看到,如果周期卷积的周期L<N1+N2-1,那么y1(n)的周期延拓必然有一部分非零序列值会交叠起来,从而出现混叠现象。只有在L≥N1+N2-1时,才没有交叠现象,而圆周卷积正是周期卷积取主值序列。 y(n)=x1(n)x2(n)= y~(n)RL(n)(341) 因此 y(n)=∑∞r=-∞y1(n+rL)RL(n)(342) 所以要使圆周卷积等于线性卷积而不产生混叠的必要条件为L≥N1+N2-1。线性卷积与圆周卷积关系如表31所示。 表31线性卷积与圆周卷积的关系 圆 周 卷 积 线 性 卷 积 1. 针对DFT引出的一种表示方法 1. 信号通过线性系统时,信号输出等于输入与系统单位脉冲响应的卷积 2. 两序列长度必须相等,若不等则按要求补零 2. 两序列长度可以相等,也可以不相等 3. 卷积结果长度与两信号长度相等皆为N 3. 卷积结果长度N=N1+N2-1 下面以具体示例说明,两个有限长矩形序列x1(n)与x2(n),它们长度分别为N1=4、N2=5,分别计算线性卷积和圆周卷积结果如图323所示。其中,图322(c)为线性卷积结果,图322(d)~(f)分别为长度为6、8、10的圆周卷积结果,根据L≥N1+N2-1关系,确实可以看到当L=8时线性卷积与圆周卷积结果相等。 图323线性卷积与圆周卷积示例 【例38】一个有限长序列为xn=δ(n)+2δ(n-5), (1) 计算序列x(n)的10点离散傅里叶变换。 (2) 若序列y(n)的DFT为Yk=ej2k2π10Xk,式中,X(k)是x(n)的10点离散傅里叶变换,求序列y(n)。 (3) 若10点序列y(n)的10点离散傅里叶变换是Yk=XkWk,式中, X(k)是序列x(n)的10点DFT,W(k)是序列w(n)的10点DFT。 wn=1,0≤n≤6 0,其他 求序列y(n)。 【解】(1) 由式(327)可求得x(n)的10点DFT: X(k)=∑N-1n=0x(n)WnkN=∑10-1n=0[δ(n)+2δ(n-5)]Wnk10 =1+2W5k10=1+2e-j2π105k=1+2(-1)k (2) 根据移位性质DFT[x((n+m))NRN(n)]=W-mkNX(k),X(k)乘以一个W-mkN形式的复指数相当于是x(n)圆周移位m点。本题中m=2, x(n)向左圆周移位了2点,就有 yn=xn+210R10n=2δ(n-3)+δ(n-8) (3) X(k)乘以W(k)相当于x(n)与w(n)的圆周卷积。为了进行圆周卷积,可以先计算线性卷积再将结果周期延拓并取主值序列。x(n)与w(n)的线性卷积为 zn=xnwn=1,1,1,1,1,3,3,2,2,2,2,2 根据式(342),圆周卷积为 y(n)=∑∞r=-∞z(n+10r)R10(n) 在0≤n≤9求和中,仅有序列z(n)和z(n+10)有非零值,用表列出z(n)和z(n+10)的值,对n=0, 1, 2, …, 9求和,得到 n0123456789 1011 Z(n) 111113322222 z(n+10)2200000000 00 y(n) 3311133222 —— 所以10点圆周卷积为: yn=3,3,1,1,1,3,3,2,2,2,由于6+7-1=12> 10所以线性卷积不等于圆周卷积。 【例39】已知x1(n)={1,1,1,1,0,2},x2(n)={0,1,2,1,0,3},n=0,1,…,5。求圆周卷积y1(n)=x1(n)⑥x2(n),y2(n)=x1(n)⑨x2(n),y3(n)=x1(n)x2(n)。 【解】首先计算线性卷积和yl(n)={0,1,3,4,4,6,6,7,5,0,6},n=0,1, …,10。yl(n)长度为11。 然后,将yl(n)以6为周期进行周期延拓得到y~(n)。由于周期延拓的周期小于yl(n)的长度,所以,在周期延拓时每个周期会有(N1+N2-1)-L个混叠点。当L=6,N1+N2-1=11,每个周期有5个混叠点,得到 y~(n)={…,6,8,8,4,10,6,…} 最后取y~(n)的主值序列即为圆周卷积和 y1(n)=x1(n)⑥x2(n)={6,8,8,4,10,6},n=0,1,…,5 同理,可以得到 y2(n)=x1(n)⑨x2(n)={0,7,3,4,4,6,6,7,5},n=0,1,…,8 当L≥N1+N2-1时,y(n)=yl(n),所以 y3(n)=x1(n)x2(n)=yl(n)={0,1,3,4,4,6,6,7,5,0,6},n=0,1,…,10 7. DFT形式下的帕塞瓦尔定理 ∑N-1n=0x(n)y(n)=1N∑N-1k=0X(k)Y(k)(343) 证明: ∑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.5频域采样理论 考虑一个任意的绝对可和的非周期序列x(n),它的Z变换为 X(z)=∑∞n=-∞x(n)z-n 由于绝对可和,所以其傅里叶变换存在且连续,故Z变换收敛域包括单位圆。如果对X(z)在单位圆上进行N点等距采样: X(k)=X(z)z=W-kN=∑∞n=-∞x(n)WnkN,k=0,1,2,…,N-1(344) 问题在于,这样采样以后是否仍能不失真地恢复原序列x(n)。也就是说,频率采样后从X(k)的反变换中所获得的有限长序列,即xNn=IDFT[Xk],能不能代表原序列x(n)?为此,先来分析X(k)的周期延拓序列X~(k)的离散傅里叶级数的反变换,令其为x~N(n)。 x~N(n)=IDFS[X~(k)]=1N∑N-1k=0X~(k)W-nkN=1N∑N-1k=0X(k)W-nkN(345) 将式(327)代入式(345),可得 x~N(n)=1N∑N-1k=0∑∞m=-∞x(m)WmkNW-nkN=∑∞m=-∞x(m)1N∑N-1k=0W(m-n)kN(346) 由于 1N∑N-1k=0W(m-n)kN=1,m=n+rN,r为任意整数 0,其他(347) 所以 x~N(n)=∑∞r=-∞x(n+rN)(348) 这说明由X~(k)得到的周期序列x~N(n)是原非周期序列x(n)的周期延拓,其时域周期为频域采样点数N。在第1.2节中已经知道,时域采样造成频域的周期延拓,这里又看到一个对称的特性,即频域采样同样会造成时域的周期延拓。 (1) 如果x(n)是有限长序列,点数为M,则当频域采样不够密,即当N<M时,x(n)以N为周期进行延拓,就会造成混叠。这时,从x~N(n)就不能不失真地恢复出原信号x(n)。因此,频域采样不失真的条件是频域采样点数N大于或等于时域采样点数M(时域序列长度),即满足 N≥M 此时可得到 xN(n)= x~N(n)RN(n)=∑∞r=-∞x(n+rN)RN(n)=x(n),N≥M(349) 也就是说,点数为N(或小于N)的有限长序列,可以利用它的Z变换在单位圆上的N个等间隔点上的采样值精确地表示。 (2) 如果x(n)不是有限长序列(即无限长序列),时域周期延拓后,必然造成混叠现象,因而一定会产生误差; 当n越大时信号衰减得越快,或频域采样越密(即采样点数N越大),则误差越小,即xN(n)越接近x(n)。 【例310】已知一个序列x(n)为5点矩形序列,其序列傅里叶变换如图324所示,分析频域上的5点抽样。 图324矩形序列及其DTFT 【解】序列x(n)时域是有限长非周期的,所以频域是连续信号。现在频域上进行抽样处理,使其频域离散化。按N=5点进行频域抽样,由于频域抽样会造成时域延拓相加,时域延拓的周期个数等于频域的抽样点数N=5,由于N=M,所以时域延拓恰好无混叠现象,如图325所示。 图325频域抽样时域周期延拓(N=M) 按N=4时进行抽样,由于N=4,序列长度为M=5,N<M,时域延拓后产生混叠现象如图326所示。 图326频域抽样时域周期延拓(N<M) 既然N个频域采样X(k)能不失真地代表N点有限长序列x(n),那么这N个采样值X(k)也一定能够完全地表达整个X(z)及频率响应X(ejω)。讨论如下: X(z)=∑N-1n=0x(n)z-n 由于 x(n)=1N∑N-1k=0X(k)W-nkN(350) X(z)=∑N-1n=01N∑N-1k=0X(k)W-nkNz-n=1N∑N-1k=0X(k)∑N-1n=0W-nkNz-n =1N∑N-1k=0X(k)1-W-NkNz-N1-W-kNz-1(351) 由于W-NkN=1,因此 X(z)=1-z-NN∑N-1k=0X(k)1-W-kNz-1(352) 这就是用N个频率采样X(k)表示X(z)的内插公式。它可以表示为 X(z)=∑N-1k=0X(k)Φk(z)(353) Φk(z)=1N1-z-N1-W-kNz-1(354) 式(354)称为内插函数。令其分子为零,得 z=ej2πNr,r=0,1,2,…,k,…,N-1(355) 图327内插函数的零极点分布 即内插函数在单位圆的N等分点上(也即采样点上)有N个零点。而分母为零,则有z=W-kN=ej2πNk的一个极点,它将与第k个零点相抵消。因而,插值函数Φk(z)只在本身采样点r=k处不为零,在其他(N-1)个采样点r上(r=0,1,2,…, N-1,r≠k)都是零点(有N-1个零点)。而它在z=0处还有(N-1)阶极点,如图327所示。 现在讨论频率响应,即求单位圆上z=ejω的Z变换。 X(ejω)=∑N-1k=0X(k)Φk(ejω)(356) 将z=ejω代入式(354),Φk(ejω)可以表示成更方便的形式: Φk(ejω)=1N1-e-jωN1-e-jω-k2πN=1NsinωN2sinω-2πNk2e-jN-12ω+kπN =1NsinNω2-πNksinω2-πNkejkπN(N-1)e-jN-12ω 这样 Φk(ejω)=Φω-k2πN(357) 其中, Φ(ω)=1Nsin(ωN/2)sin(ω/2)e-jN-12ω(358) 式(358)称为频域响应的内插函数。而式(356)又可改写为 X(ejω)=∑N-1k=0X(k)Φω-2πNk(359) Φω-k2πN满足以下关系 Φω-k2πN=1,ω=k2πN=ωk 0,ω=i2πN=ωi,i≠k(360) 也就是说,函数Φω-k2πN在本采样点ωk=k2πN上,Φωk-k2πN=1,而在其他采样点ωi=i2πN,i≠k上,函数Φωi-k2πN=0。整个X(ejω)就是由N个Φω-k2πN函数分别乘上X(k)后求和的。所以很明显,在每个采样点上X(ejω)就精确地等于X(k)(因为其他点的插值函数在这一点上的值为零,没有影响)即 X(ejω)ω=2πNk=X(k),k=0,1,2,…,N-1(361) 就是说,各采样点之间的X(ejω)值由各采样点的加权插值函数X(k)Φω-2πNk在所求ω点上的值的叠加得到的,如图328所示。 图328内插函数幅度特性与相位特性(N=5) 【例311】频域采样定理的验证,给定信号如下: xn=n+1,0≤n≤13 27-n,14≤n≤26 0,其他 编写程序分别对频谱函数X(ejω)=FT[x(n)]在区间0,2π上等间隔采样32和16点,得到 X32(k)=X(ejω)ω=2π32k,k=0,1,2,…,31 X16(k)=X(ejω)ω=2π16k,k=0,1,2,…,15 再分别对X32k和X16k进行32点和16点IFFT,得到 x32(n)=IFFT[X32(k)]32,n=0,1,2,…,31 x16(n)=IFFT[X16(k)]16,n=0,1,2,…,15 分别画出X(ejω)、X32k和X16k的幅度谱,并绘图显示x(n)、x32n和x16n的波形,进行对比和分析,验证总结频域采样理论。 提示: 频域采样用以下方法编程实现。 (1) 直接调用MATLAB函数FFT计算X32(k)=FFT[x(n)]32得到X(ejω)在0,2π的32点频率域采样。 (2) 抽取X32(k)的偶数点即可得到X(ejω)在0,2π的16点频率域采样X16(k),即X16(k)=X32(2k),k=0,1,2,…,15。 (3) 当然也可以按照频域采样理论,先将信号x(n)以16为周期进行周期延拓,取其主值区(16点),再对其进行16点DFT(FFT),得到的就是X(ejω)在0,2π的16点频率域采样X16(k)。 MATLAB代码如下: M=27;N=32;N=0:M; %产生M长三角波序列x(n) xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=[xa,xb]; Xk=fft(xn,1024); %1024点FFT[x(n)], 用于近似序列x(n)的TF X32k=fft(xn,32); %32点FFT[x(n)] x32n=ifft(X32k); %32点IFFT[X32(k)]得到x32(n) X16k=X32k(1:2:N); %隔点抽取X32k得到X16(K) x16n=ifft(X16k,N/2); %16点IFFT[X16(k)]得到x16(n) subplot(3,2,2);stem(n,xn,'.');box on title('(b)三角波序列x(n)');xlabel('n');ylabel('x(n)'); axis([0,32,0,20]);k=0:1023;wk=2*k/1024; subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]'); xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,1,0,200]) k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box on title('(c) 16点频域采样');xlabel('k');ylabel('|X_1_6(k)|'); axis([0,8,0,200]);n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box on title('(d) 16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)'); axis([0,32,0,20]);k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box on title('(e) 32点频域采样');xlabel('k');ylabel('|X_3_2(k)|'); axis([0,16,0,200]);n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box on title('(f) 32点IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)'); axis([0,32,0,20]) function tstem(xn,yn) %时域序列绘图函数 % xn:信号数据序列,yn:绘图信号的纵坐标名称(字符串) n=0:length(xn)-1; stem(n,xn,'.');box on xlabel('n');ylabel(yn); axis([0,n(end),min(xn),1.2*max(xn)]) 结果如图329所示。图329(a)和(b)分别为X(ejω)和xn的波形; 图329(c)和(d)分别为X(ejω)的16点采样X16k和x16n=IDFTX16k波形; 图329(e)和(f)分别为X(ejω)的32点采样X32k和x32n=IDFTX32k波形图; 由于实序列DFT满足共轭对称性,因此频域图仅画出0,π上的幅频特性波形。本例中x(n)的长度M=26。从图中可以看出,当采样点数N=16<M时,x16n确实等于原三角序列xn以16为周期的周期延拓序列的主值序列。由于存在时域混叠失真,因而x16n≠xn; 当采样点数N=32>M时,无时域混叠失真,x32n=IDFTX32k=xn。 图329例311图 视频讲解 3.6MATLAB应用实例 【例312】x[k]=cos(2πrk/N), N=16, r=4,利用MATLAB计算16点序列x[k]的512点DFT,如图330所示。 图330x[k]的512点DFT MATLAB代码如下: N = 16;k = 0:N-1; L = 0:511;x = cos(2*pi*k*4./16); X = fft(x);plot(k/16,abs(X),'o'); hold on;XE = fft(x,512); plot(L/512,abs(XE)) xlabel('归一化频率');ylabel('幅度'); 【例313】离散傅里叶变换X(k)=∑N-1n=0x(n)WnkN(矩阵相乘的方法)。 MATLAB代码如下: function [Xk]=dft(xn,N) n=[0:1:N-1]; k=[0:1:N-1]; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^(nk); Xk=xn*WNnk; 【例314】逆离散傅里叶变换x(n)=1N∑N-1k=0X(k)W-nkN。 MATLAB代码如下: function [xn]=idft(Xk,N) n=[0:1:N-1]; k=[0:1:N-1]; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^(-nk); 【例315】信号的傅里叶分解与合成。 MATLAB代码如下: clear all N = 256; dt = 0.05; % data numbers and sampling intervel,sampling frequence is 20Hz n=0:N-1;t=n*dt; %序号序列和时间序列 x1=sin(2*pi*t);x2=0.5*sin(2*pi*5*t); x=sin(2*pi*t)+0.5*sin(2*pi*5*t); %signals add m=floor(N/2)+1; %down for integer a=zeros(1,m);b=zeros(1,m); for k=0:m-1 for ii=0:N-1 a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%matlab's array index must be increase from 1 b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N); end c(k+1)=sqrt(a(k+1).^2+b(k+1).^2); end if(mod(N,2)~=1)a(m)=a(m)/2;end for ii=0:N-1 xx(ii+1)=a(1)/2; for k=1:m-1; xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k* ii/N); end end subplot(2,2,1),plot(t, x1); title('正弦信号1'), xlabel('时间/s'); subplot(2,2,2),plot(t, x2); title('正弦信号2'), xlabel('时间/s'); %subplot(4,1,1),plot((0:N-1)*dt,xx);title('Composed sitgnal'); subplot(2,2,3),plot(t, x); title('合成信号'), xlabel('时间/s'); subplot(2,2,4),plot((0:m-1)/(N*dt),c),title('傅里叶变换'),xlabel('频率/Hz'), ylabel('幅度'); 信号傅里叶合成与分解如图331所示。 图331例315图 此处的1Hz和5Hz的振幅与原来信号振幅不完全一致,是由于数据采样点较少导致的,即N较小。 【例316】补零序列的离散傅里叶变换。 MATLAB代码如下: n=0:4; x=[ones(1,5)]; %产生矩形序列 k=0:999;w=(pi/500)*k; X=x*(exp(-j*pi/500)).^(n'*k); %计算离散时间傅里叶变换 Xe=abs(X); %取模 subplot(3,2,1);stem(n,x);ylabel('x(n)'); %画出矩形序列 subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %画出离散时间傅里叶变换 N=10;x=[ones(1,5),zeros(1,N-5)]; %将原序列补零为10长序列 n=0:1:N-1; X=dft(x,N); %进行DFT magX=abs(X); k=(0:length(magX)'-1)*N/length(magX); subplot(3,2,3);stem(n,x);ylabel('x(n)'); %画出补零序列 subplot(3,2,4);stem(k,magX); %画出DFT结果 axis([0,10,0,5]);ylabel('|X(k)|'); N=20;x=[ones(1,5),zeros(1,N-5)]; %将原序列补零为20长序列 n=0:1:N-1; X=dft(x,N); %进行DFT magX=abs(X); k=(0:length(magX)'-1)*N/length(magX); subplot(3,2,5);stem(n,x);ylabel('x(n)'); %画出补零序列 subplot(3,2,6);stem(k,magX); %画出DFT结果 axis([0,20,0,5]);ylabel('|X(k)|'); 补零序列的离散傅里叶变换如图332所示。 图332补零序列的离散傅里叶变换 序列末端补零后,尽管信号的频谱不会变化,但对序列做补零后再做L点DFT,计算出的频谱实际上是原信号频谱在[0,2π)区间上L个等间隔采样,从而增加了对真实频谱采样的点数,并改变了采样点的位置,将会显示出原信号频谱更多的细节。故而数据后面补零可以克服栅栏效应。在第4章快速傅里叶变换中将详细介绍。 【本章习题】 3.1填空题 (1) 已知一个长度为N的序列xn,它的离散时间傅里叶变换为Xejω,它的N点离散傅里叶变换Xk是关于Xejω的点等间隔。 (2) DFT与DFS有密切关系,因为有限长序列可以看成周期序列的,而周期序列可以看成有限长序列的。 (3) 对长度为N的序列xn圆周移位m位得到的序列用xmn表示,其数学表达式为xmn= 。 (4) 设序列x(n)的N点DFT为X(k),则x((n+m))NRN(n)的N点DFT为。 (5) 某序列的DFT表达式为X(k)=∑N-1n=0x(n)WknN,由此可以看出,该序列时域的长度为,变换后数字频域上相邻两个频率样点的间隔是。 (6) 圆周卷积可被看作是周期卷积的; 圆周卷积的计算是在区间中进行的,而线性卷积不受此限制。 (7) 有限长序列Xz与Xk的关系,Xk与X(ejω)的关系。 3.2选择题 (1) 序列x1n的长度为4,序列x2n的长度为3,则它们线性卷积的长度是(),5点圆周卷积的长度是()。 A. 5,5B. 6,5C. 6,6D. 7,5 (2) 下面描述中最适合离散傅里叶变换DFT的是()。 A. 时域为离散序列,频域也为离散序列 B. 时域为离散有限长序列,频域也为离散有限长序列 C. 时域为离散无限长序列,频域为连续周期信号 D. 时域为离散周期序列,频域也为离散周期序列 (3) 若序列的长度为M,要能够由频域抽样信号X(k)恢复原序列,而不发生时域混叠现象,则频域抽样点数N需满足的条件是()。 A. N≥MB. N≤MC. N≤2MD. N≥2M 3.3x(n)和h(n)是如下给定的有限序列x(n)={5,2,4,-1,2},h(n)={-3,2,-1}。 (1) 计算x(n)和h(n)的线性卷积y(n)= x(n)*h(n); (2) 计算x(n)和h(n)的6点循环卷积y1(n)= x(n)⑥h(n); (3) 计算x(n)和h(n)的8点循环卷积y2(n)= x(n)⑧h(n); 比较以上结果,有何结论? 3.4证明W(N-n)kN=W-nkN=(WnkN)*。 3.5对有限长序列x(n)=1,0,1,1,0,1的Z变换X(z)在单位圆上进行5等分取样,得到取样值X(k),即X(k)=X(z)z=W-k5,k=0,1,2,3,4求X(k)的逆傅里叶变换x1(n)。 3.6试用定义计算周期为5,且一个周期内x(n)={2-,1,3,0,4}的序列x~(n)的DFS。 3.7设 x(n)=1,n=0,1 0,其他 将x(n)以4为周期进行周期延拓,形成周期序列x~(n),画出x(n)和x~(n)的波形,求出x~(n)的离散傅里叶级数X~(k)和傅里叶变换。 3.8已知序列x(n)=δ(n)+2δ(n-2)+δ(n-3),若y(n)是x(n)与其本身的4点循环卷 积,求y(n)及其4点DFT Y(k)。 3.9序列x(n)为 x(n)=2δ(n)+δ(n-1)+δ(n-3) 计算x(n)的5点DFT,然后对得到的序列求平方: Y(k)=X2(k) 求Y(k)的5点DFT反变换y(n)。 3.10设序列x(n)={1,2,3,2,1,0},v(n)={3,2,1,0,1,2}。 (1) 求x(n)的傅里叶变换X(ejω); (2) 求v(k)=DFT[v(n)]6; (3) 请解释v(k)与X(ejω)之间的关系。