余成波matlab本书稿有关程序文档.doc

1.6.1 离散时间信号的MATLAB实现1正弦序列离散正弦序列的MATLAB表示与连续信号类似，只不过是用stem函数而不是用plot函数来画出序列的波形。下面就是正弦序列的MATLAB源程序。程序运行结果如图1.19所示。%正弦序列实现程序k=0:39;fk=sin(pi/6*k);axes(handles.axes1)stem(k,fk)图1.19 正弦序列波形2指数序列离散指数序列的一般形式为，可用MATLAB中的数组幂运算（即点幂运算）c*来实现。下面为用MATLAB编写绘制离散时间实指数序列波形的函数。function dszsu(c,a,k1,k2)%c：指数序列的幅度%a：指数序列的底数%k1：绘制序列的起始序号%k2：绘制序列的终止序号k=k1:k2;x=c*(a.k);stem(k,x,filled)hold onplot(k1,k2,0,0)hold off利用上述函数，实现实指数波形MATLAB程序如下（其中值分别为）。%离散时间实指数序列实现程序subplot 221;dszsu(1,5/4,0,20);xlabel(k);title(f1k);subplot 222dszsu(1,3/4,0,20);xlabel(k);title(f2k);subplot 223;dszsu(1,-5/4,0,20);xlabel(k);title(f3k);subplot 224;dszsu(1,-3/4,0,20);xlabel(k);title(f4k);程序运行结果如图1.20所示。如图可知，对于离散时间实指数序列，当的绝对值大于1时，序列为随时间发散的序列，当的绝对值小于1时，序列为随时间收敛的序列。同时可见，当的值小于零时，其波形在增长或衰减的同时，还交替地改变序列值的符号。图1.20 不同底数的实指数序列对于离散时间虚指数序列，可用通过调用下列绘制虚指数序列时域波形的MATLAB函数。function=dxzsu(n1,n2,w)%n1:绘制波形的虚指数序列的起始时间序号%n2：绘制波形的虚指数序列的终止时间序号%w：虚指数序列的角频率k=n1:n2;f=exp(i*w*k);Xr=real(f)Xi=imag(f)Xa=abs(f)Xn=angle(f)subplot(2,2,1), stem(k,Xr,filled),title(实部);subplot(2,2,3), stem(k,Xi,filled),title(虚部);subplot(2,2,2), stem(k,Xa,filled),title(模);subplot(2,2,4), stem(k,Xn,filled),title(相角);利用上述函数，实现虚指数波形MATLAB程序如下（其中虚指数分别为）%离散时间虚指数实现程序figure(1);dxzsu(0,20,pi/4);figure(2);dxzsu(0,20,2);程序运行结果如图1.21（a）、（b）所示。由图可见，只有当虚指数序列的角频率满足为有理数时，信号的实部和虚部和相角都为周期序列，否则为非周期序列。 （a）波形 (b) 波形图1.21 虚指数序列波形对于复指数序列，其一般形式为可以通过调用下面绘制复指数序列时域波形的MATLAB函数。function dfzsu(n1,n2,r,w)%n1:绘制波形的虚指数序列的起始时间序号%n2：绘制波形的虚指数序列的终止时间序号%w：虚指数序列的角频率%r: 指数序列的底数k=n1:n2;f=(r*exp(i*w).k;Xr=real(f);Xi=imag(f);Xa=abs(f);Xn=angle(f);subplot(2,2,1), stem(k,Xr,filled),title(实部);subplot(2,2,3), stem(k,Xi,filled),title(虚部);subplot(2,2,2), stem(k,Xa,filled),title(模);subplot(2,2,4), stem(k,Xn,filled),title(相角);利用上述函数，实现复指数序列波形MATLAB程序如下。%复指数序列实现程序(r1)figure(1);dfzsu(0,20,1.2,pi/4);%复指数序列实现程序(0r1时，复指数序列的实部和虚部分别为幅度按指数增长的正弦序列;当0r1 (b) 0r=0; %序列的起点为ns，终点为nf，在n=n0点处生成单位阶跃。（3）两个信号相加的生成函数sigadd.mfunction y,n=sigadd(x1,n1,x2,n2)n=min(min(n1),min(n2):max(max(n1),max(n2)y1=zeros(1,length(n);y2=y1;y1=(find(n=min(n1)&(n=min(n2)&(n=min(n1)&(n=min(n2)&(n=max(n2)=1)=x2;y=y1.*y2;（5）序列移位的生成函数sigshift.mfunction y,n=sigshift(x,m,n0)n=m+n0;y=x（6）序列翻折的生成函数sigfold.mfunction y,n=sigfold(x,n)y=fliplr(x);n=-fliplr(n)（7）奇偶综合函数evenodd.m此函数可以将任一给定的序列x(n)分解为xe(n) 和xo(n)两部分。function xe,xo,m=evenodd(x,n)if (imag(x)=0) error(x is not a real sequence)； endm=-fliplr(n);m1=min(m,n);m2=max(m,n);m=m1:m2;nm=n(1)-m(1);n1=1:length(n);x1=zeros(1,length(m);x1(n1+nm)=x;x=x1;xe=0.5*(x+fliplr(x); xo=0.5*(x-fliplr(x);(8）求卷积和求卷积和直接采用MATLAB中的函数conv，即y=conv(x,h)；它默认序列从n=0开始。但是如果序列是从一负值开始，即其中nx10或nh1=min(k1)&(k=min(k2)&(k=min(k1)&(k=min(k2)&(k=0;subplot 221stem(n,x)xlabel(n);ylabel(x(n);title(Step Sequence);grid on;%Decomposition of the Sequencexeven,xodd,m=sigevenodd(x,n);subplot 223stem(m,xeven);xlabel(m);ylabel(x even(n);title(Even Part);grid on;subplot 224stem(m,xodd);xlabel(m);ylabel(x odd(n);title(Odd Part);grid on;程序运行结果如图1.30所示。例1.11 已知x(n)=u(n)-u(n-10),要求将序列分解为奇偶序列。解：程序清单如下。n=0:10;x=stepseq(0,0,10)- stepseq(10,0,10);xe,xo,m=evenodd(x,n);subplot(2,2,1);stem(n,x);ylabel(x(n);title(矩形序列);axis(-10,10,-1.2,1.2)subplot(2,2,2);stem(m,xe);ylabel(xe(n);title(奇序列);axis(-10,10,-1.2,1.2)subplot(2,2,4);stem(m,xo);ylabel(xo(n);title(偶序列);axis(-10,10,-1.2,1.2)程序运行结果如图1.31所示。图1.31 序列分解为奇偶序列1.6.4 常系数线性差分方程解的MATLAB实现例1.12 设线性时不变系统的抽样响应h(n)=(0.9)nu(n),输入序列x(n)=u(n)-u(n-10),求系统的输出y(n)。解：系统的输出y(n)为输入x(n)与系统的单位抽样响应h(n的卷积，即y(n)=x(n)* h(n),可直接采用con_m函数求输出序列。程序清单如下n=-5:50;x=stepseq(0,-5,50)-stepseq(10,-5,50);h=(0.9).n).* stepseq(0,-5,50);subplot(3,1,1);stem(n,x);axis(-5,50,0,2);ylabel(x(n);subplot(3,1,2);stem(n,h);axis(-5,50,0,2);ylabel(h(n);y,ny=conv_m(x,n,h,n);subplot(3,1,3);stem(ny,y);axis(-5,50,0,8); xlabel(n);ylabel(y(n);程序运行结果如图1.32所示。图1.32 线性时不变系统的输出例1.13 给定因果稳定线性时不变系统的差分方程对下列输入序列x(n),求出输出序列y(n)。（1）x(n)=(n)；（2）x(n)=u(n)； （3）x(n)=R30(n)； （4）解：程序清单如下。N=100;n=0:N-1;b=1;a=1,-1,0.9;x1=n=0;y1=filter(b,a,x1);subplot(3,2,1);stem(n,y1,.);axis(0,N,-1,2);ylabel(y1(n);x2=n=0;y2=filter(b,a,x2);subplot(3,2,2);stem(n,y2,.);axis(0,N,-1,2);ylabel(y2(n);x3=(n=0)&(n2 fh，因此采样信号的频谱发生混叠造成取样信号不能无失真的恢复出原始信号。图1.37 从x1(n)和x2(n)分别重构模拟信号2.6 用MATLAB进行离散系统的Z域分析2.6.1 MATLAB几个信号处理工具箱函数1residuez:求极点留数分解前面已经介绍，离散系统的输出Y(z)=X(z)H(z)=B(z)/A(z)他必然是z的有理分式从而得出其时域信号为其中k是当MN时的直接项，即有限序列，而其余的都是无限序列。其中的r,p,k向量可以由下列语句求得r,p,k=residuez(B,A)2impz：求H(z)的Z反变换h(n)h,T=impz(B,A,N)h为存放h(n)的列向量，时间变量N存放在列向量T中，当N为标量时，表示T=0：N-1，计算h(n)，n=0，1，2，N-1；当N为向量时，T=N，仅计算N指定的整数点上的h(n)。3freqz：求数字滤波器H(z)的频率响应函数H=freqzB,A,w计算由向量w指定的数字频率点上数字滤波器H(z)的频率响应，结果存在H向量中。H,w=freqzB,A,M计算出M个频率点上的频率响应，存放在H向量中，M个频率存放在w向量中，freqz函数自动将这M个频率点均匀设置在频率范围0，之间。若缺省w和M时，函数自动选取512个频率点计算。不带输出向量的freqz函数将自动绘制幅频和相频曲线。4zplane:绘制H(z)的零极点图zplane(z,p)绘制出列向量z中的零点（以符号“”表示）和列向量p中的极点（以符号“”表示）以及参考单位圆，并在多阶零点和极点的右上角标出其阶数。如果z和p为矩阵，则会以不同颜色绘出z和p各列中的零点和极点。zplane(B,A)绘制出系统函数H(z)的零极点图。 2.6.2 利用MATLAB绘制离散系统的零极图对于离散系统其系统函数可由式（2.46）表示，则系统函数的零点和极点可以用MATLAB的多项式求根函数roots()来实现，调用该函数的命令格式为：p=roots(A)其中A为待求根的多项式的系数构成的行向量，返回向量p则包含该多项式所有的根位置列向量。值得注意的是：求系统函数零极点时，离散系统的系统函数可能有两种形式，一种是分子和分母多项式按Z的降幂次序排列，另一种是分子和分母多项式按的升幂次序排列。若是以Z的降幂次序排列，则系数向量一定要由多项式的最高幂次开始，一直到常数项，缺项要用0补齐; 若是以按的升幂次序排列，则分子和分母多项式系数向量的维数一定要相同，不足的要用0补齐，不则的零点或极点就可能被漏掉。下面是求系统函数零极点，并绘制其零极点图的MATLAB实用函数ljdt(),该函数在绘制系统零极点图的同时，还绘出了Z平面的单位圆。function ljdt(A,B)% The function to draw the pole-zero diagram for discrete systemp=roots(A); %求系统极点q=roots(B); %求系统零点p=p; %将极点列向量转置为行向量q=q; %将零点列向量转置为行向量x=max(abs(p q 1); %确定纵坐标范围x=x+0.1;y=x; %确定横坐标范围clfhold onaxis(-x x -y y) %确定坐标轴显示范围w=0:pi/300:2*pi;t=exp(i*w);plot(t) %画单位园axis(square)plot(-x x,0 0) %画横坐标轴plot(0 0,-y y) %画纵坐标轴text(0.1,x,jImz)text(y,1/10,Rez)plot(real(p),imag(p),x) %画极点plot(real(q),imag(q),o) %画零点title(pole-zero diagram for discrete system) %标注标题hold off例2.24 已知某离散系统的系统函数为试用MATLAB求出该系统的零极点，并画出零极点分布图，判断系统是否稳定。图2.10 例2.24系统零极点图解 利用ljdt()函数，MATLAB程序为% 绘制零极点分布图的实现程序a=3 -1 0 0 0 1;b=1 1;ljdt(a,b)p=roots(a)q=roots(b)pa=abs(p)程序运行如果为：p = 0.7255 + 0.4633i 0.7255 - 0.4633i -0.1861 + 0.7541i -0.1861 - 0.7541i -0.7455 q = -1pa = 0.8608 0.8608 0.7768 0.7768 0.7455绘制的系统零极点如图2.10所示。由图可知，该系统的所有极点均位于Z平面的单位圆内，故该系统为稳定系统。例2.25 已知某离散系统的系统函数为试用MATLAB求出该系统的零极点，并画出零极点分布图，求系统的单位冲激响应和幅频响应，并判断系统的是否稳定。解 MATLAB实现程序为% 由系统函数求解系统脉冲响应,频率响应实现程序b=0 1 2 1;a=1 -0.5 -0.005 0.3;subplot 311zplane(b,a);xlabel(实部);ylabel(虚部);num=0 1 2 1;den=1 -0.5 -0.005 0.3;h=impz(num,den);subplot 312stem(h);xlabel(k);title(单位脉冲响应);H,w=freqz(num,den);subplot 313plot(w/pi,abs(H);xlabel(频率omega);title(频率响应)程序运行结果如图2.11所示。该系统是稳定的。图2.11 零极点分布图，求系统的单位冲激响应和幅频响应2.6.3 利用MATLAB分析离散系统的零极图分布与系统单位响应时域特性的关系例2.26 已知离散系统的零极分布分别如图2.12所示，其中虚线表示单位圆，试用MATLAB分析系统单位响应的时域特性。解 因系统的零极点分布图已知，则系统的系统函数就可知，故可以利用MATLAB函数impz()函数可求出系统单位响应。对于图2.12所示的系统，其系统函数分别为、则其单位响应的时域波形MATLAB程序为% 零极点分布与单位响应的关系实现程序a=1,-1;b=1;subplot 321impz(b,a);a1=1,-0.8;b1=1;subplot 322impz(b1,a1,10);a2=1,-2;b2=1;subplot 323impz(b2,a2,10);a3=1,-2*0.8*cos(pi/4),0.82;b3=1;subplot 324impz(b3,a3,20);a4=1,-2*0.8*cos(pi/8),1;b4=1;subplot 325impz(b4,a4,20);a5=1,-2*1.2*cos(pi/4),1.22;b5=1;subplot 326impz(b5,a5,20);程序运行结果如图2.13所示。图2.12 离散系统的零极点分布图由此可知：离散系统的单位响应的时域特性完全由系统函数的极点位置决定，位于Z平面单位圆内的极点决定了随时间衰减的信号分量，位于Z平面单位圆上的极点决定了单位响应的稳定信号分量，位于Z平面单位圆外的极点决定了单位响应随时间增长的信号分量。图2.13 系统单位响应的时域波形图例2.27 求一因果线性移不变系统y(n)=0.81y(n-2)+x(n)-x(n-2)的单位抽样响应h(n)，单位阶跃响应g(n)，及绘出H(ej)的幅频和相频特性。解:程序清单如下b=1,0,-1;a=1,0,-0.81;figure(1)subplot(2,1,1);dimpulse(b,a,50);ylabel(h(n);subplot(2,1,2);dstep(b,a,50);ylabel(g(n);figure(2)w=0:1:500*pi/500;freqz(b,a,w)程序运行结果如图2.14所示。图2.14 系统的抽样响应、阶跃响应及幅频、相频曲线例2.28 梳状滤波器零极点和幅频特性。梳状滤波器系统函数有如下两种类型：FIR型： IIR型：分别令N=8，a=0.8,0.9,0.98计算并图示和的零、极点图及幅频特性，说明极点位置的影响。解:程序清单如下b=1,0,0,0,0,0,0,0,-1; %H1(z)和H2(z)的分子多项式系数向量a0=1; %H1(z)分母多项式系数向量a1=1,0,0,0,0,0,0,0,-(0.8)8; %H2(z)分母多项式系数向量（a=0.8）a2=1,0,0,0,0,0,0,0,-(0.9)8; %H2(z)分母多项式系数向量（a=0.9）a3=1,0,0,0,0,0,0,0,-(0.98)8; %H2(z)分母多项式系数向量（a=0.98）H,w=freqz(b,a0); %H1(z)的频响函数H1,w1=freqz(b,a1);H2,w2=freqz(b,a2);H3,w3=freqz(b,a3);subplot(4,2,1);zplane(b,a0); xlabel(实部);ylabel(虚部);title(FIR梳状滤波器零点图)subplot(4,2,2);zplane(b,a1);xlabel(实部);ylabel(虚部);title(IIR梳状滤波器零极点图 a=0.8)subplot(4,2,3);plot(w/pi,abs(H);title(FIR梳状滤波器幅频响应曲线)subplot(4,2,4); plot(w/pi,abs(H1);title(IIR梳状滤波器幅频响应曲线 a=0.8).（作图部分的程序省略，可参考前面的程序）程序运行结果如图2.15所示。由图看出，阶数相同的时候，IIR滤波器具有更平坦的通带特性和更窄的过渡带，极点距离单位圆越近，这个特性就越明显。图2.15 八阶梳状滤波器零、极点图及幅频响应曲线2.6.4 利用系统函数求解离散系统差分方程的MATLAB例2.29 求解差分方程，其中，初始状态，。解 将方程两边进行变换得用MATLAB编程如下% 求差分方程的解的实现程序num=0.45 0.4 -1;den=1 -0.4 -0.45;x0=1 2;y0=0 1;N=50;n=0:N-1;x=0.8.n;Zi=filtic(num,den,y0,x0);y,Zf=filter(num,den,x,Zi);plot(n,x,r-,n,y,b-);title(响应);xlabel(n);ylabel(x(n)-y(n);legend(输入 x,输出 y,1);grid;程序运行结果如图2.16所示。图2.16 离散系统差分方程解2.6.5 利用MATLAB实现域的部分分式展开式MATLAB的信号处理工具箱提供了一个对进行部分分式展开的函数residuez()，其调用形式为r,p,k=r