MATLAB讲义第二章

第二章 序列的付氏变换（DTFT ）和Z 变换一、序列的DTFT ：X （e jw ）序列x(n)的付氏变换如下：∑+∞∞--=n j j e n x e X ωω)()( ∞≤≤∞-ω ωπωππωd e e X n x n j j ⎰-=)(21)(因为MATLAB 无法计算连续变量w,只能在πωπ<≤-（或考虑对称性，在0<≤-ωπ）范围内，把w 赋值为很密的、长度很长的向量来近似连续变量。

通常最简单的就是赋以K 个等间隔的值： Kk d k πωω2=⋅= 1:0-=K k 它表示把基本数字频率的范围π2均分成K 份后，每一份的大小。

则频率向量表示为w=[w1,w2,…,wk]=k ∙dw,序列的位置向量为n=[n1,n2,…,nN]则DTFT 的计算式可用一个向量与矩阵相乘的运算来实现：[X(w 1),X(w 2),…,X(w k )]=[x(n 1),x(n 2),…,x(n N )]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡∙N k N N k k n j n j n j n j n j n j n j n j n j e e e e e e e e e ωωωωωωωωω 212222111211 =[x(n 1),x(n 2),…,x(n N )]∙w *jn T e又因为 (jn T *w)=jn T *kdw=j*dw*n ’*k 所以Xw=xn*exp(j*dw*n ’*k)例1：求有限长序列xn=[1 3 5 3 1]，n=-1:3的DTFT,画出它在w=-8~8rad/s 范围内的频率特性，讨论其对称性。

再把xn 左右移动，讨论时移对DTFT 的影响。

解:根据DTFT 定义得ωωωωωω32353)()(j j f f n n j j e e e e e n x e X ---∞-∞=-++++==∑ 88≤≤-ω将w 在-8~8rad/s 之间分为1000份。

>> xn=[1 3 5 3 1];nx=-1:3;>> w=linspace(-8,8,1000);%设定频率向量¿>> Xw=xn*exp(-j*nx'*w);%用定义计算DTFT>>subplot(3,5,1);stem(nx,xn);axis([-2,6,0,6]);title('Ô-ʼÐòÁÐ');ylabel('xn');%画序列图>> subplot(3,5,2);plot(w,abs(Xw));title('·ù¶È');% 画幅频曲线 >> subplot(3,5,3);plot(w,angle(Xw));title('Ïà½Ç')% 画相频曲线 >> subplot(3,5,4);plot(w,real(Xw));title('ʵ²¿')% 画实频曲线 >> subplot(3,5,5);plot(w,imag(Xw));title('Ð鲿')% 画虚频曲线 >> nx_2=nx+2;%使xn 右移2位>> Xw_2=xn*exp(-j*nx_2'*w);>>subplot(3,5,6);stem(nx_2,xn);axis([-2,6,0,6]);ylabel('xn_2');title('ÓÒÒÆ2λ')>> subplot(3,5,7);plot(w,abs(Xw_2));>> subplot(3,5,8);plot(w,angle(Xw_2));>> subplot(3,5,9);plot(w,real(Xw_2));>> subplot(3,5,10);plot(w,imag(Xw_2));>> nx_1=nx-1;%使xn 左移1位>> Xw_1=xn*exp(-j*nx_1'*w);>>subplot(3,5,11);stem(nx_1,xn);axis([-2,6,0,6]);ylabel('xn_1');title('×óÒÆ1λ')>> subplot(3,5,12);plot(w,abs(Xw_1));>> subplot(3,5,13);plot(w,angle(Xw_1));>> subplot(3,5,14);plot(w,real(Xw_1));>> subplot(3,5,15);plot(w,imag(Xw_1));由上述结果得出以下结论：（1） 序列的DTFT 是连续函数。

（2） 序列的DTFT 是周期函数，周期为2pi 。

（3） 本题给出的是实序列，实序列的DTFT 具有对称性，幅频和实频偶对称;相频和虚频奇对称。

（4） 信号的时移不影响幅频特性，只影响相频。

MATLAB 工具箱中有计算DTFT 的专用函数freqz 。

它的调用方式为：[H,w]=freqz(b,a,N)其中用于计算DTFT 时，b ——输入序列xn,其向量限定为n=0:(length(b)-1)(要从0开始)。

a 取1;N ——是把πω<≤0分割的份数，计算的是天频率部分的特性。

若省略N ，默认值为512。

要计算全部正负频率特性的响应，则用[H,w]=freqz(b,a,N ，‘whole ‘)二、Z 变换的计算（1）两个序列的Z 变换例1：x1(n)=[2,3,4],ns1=0;x2(n)=[3 4 5 6],ns2=0;求两序列的Z 变换，及两序列的卷积y(n)及Y(z).由Z 变换的定义∑∞-∞=-=n n z n x z X )()(知， X1(z)=2z 0+3z -1+4z-2 X2(z)=3z0+4z -1+5z -2+6z -3∵y(n)=x1(n)*x2(n)∴Y(z)=X1(z)X2(z)∴Y(z)= (2z 0+3z -1+4z -2)（3z0+4z -1+5z -2+6z -3）＝6＋17z -1+34z -2+43z -3+38z -4+24z -5则，y(n)=[6 17 34 43 38 24]y(n)= x1(n)*x2(n)=conv(x1,x2)求卷积程序：>> x1=[2 3 4];nx1=0:2;x2=[3 4 5 6];nx2=0:3;>> y=conv(x1,x2)y =6 17 34 43 38 24即y(n)=[6 17 34 43 38 24]从而Y(z)=＝6＋17z -1+34z -2+43z -3+38z -4+24z -5 由以上看出，两序列的卷积和它们的Z 变换的相乘结果是一致的。

（时域卷积定理）三、Z 反变换(1) 部分分式法MATLAB 的内部函数residuez 可计算离散系统的极点留数。

NN MM z a z a z b z b b z B z A z H ----++++++== 111101)()()( 调用格式：[r,p,C]=residuez(b,a).b ,a ――H （z ）中，分子分母多项式以z －1的升幂排列的系数向量。

p ――极点向量,即分母的根r ――部分分式分解后，各极点所对应的留数A kC ――无穷项多项式系数向量，仅在M ≥N 时存在。

则+++-++-+-==----1111)2()1()(1)()2(1)2()1(1)1()()()(z C C zN p n r z p r z p r z B z A z H则Z 反变换为：y(n)=r(1)p(1)n u (n)+…+r(N)p(N)nu(n)+C(1)δ(n)+ C(2)δ(n-1)+… 若H(z)是系统函数，则其反变换h(n)就是它的单位脉冲响应，可用函数impz 来计算，其调用方式为：h=impz(b,a,L)其中，b,a ——H(z)的按z-1的升幂排列的分子分母系数向量；L ——待求脉冲序列h(n)的长度可用该命令来核对计算。

实例1：计算的)7.01()9.01(1)(121--+-=z z z X 反变换 先用函数poly 求出分母多项式的系数，再用residuez 求X （z ）的极点和留数程序：>> b=1;a=poly([0.9,0.9,-0.7]);>> [r,p,C]=residuez(b,a)r =0.24610.56250.1914p =0.90000.9000-0.7000C =[]四、用Z 变换求系统输出(1)方法1：采用residuez 函数实例1：系统函数2112523)(---+--=z z z z H ，输入信号x=[2 3 4 5],用Z 变换求出输出y(n)。

解：∵Y （z ）=H(z)·X(z)= 211252)(3---+-⋅-z z z X z =)()(z A z B ∴ B(z)是X （z ）与－3z -1的卷积，用函数conv 可求得求出Y(z)后，用函数residuez 可求其反变换，即得y(n)。

>> x=[2 3 4 5];nsx=0;>> nfx=length(x)-1;>> b=-3;a=[2 -5 2];>> B=conv(-3,x);A=a;>> [r,p,k]=residuez(B,A)r =-10.250032.0000p =2.00000.5000k =-24.7500 -7.5000>> nf=input('nf');nf8>> n=0:nf;>> yd=k(1)*impseq(0,0,nf)+k(2)*impseq(1,0,nf) ;>> y=yi+yd >> subplot(1,2,1),stem([nsx:nfx],x);title('x'>> subplot(1,2,2),stem([0:nf],y);title('y'),line([0,10],[0,0]);(2)方法2：采用filter 函数若不需求出表达式，可用函数filter 来解。