信号与系统(第5版)MATLAB仿真举例

最好的沉淀整理MATLAB仿真举例一、MATLAB应用概要MATLAB的知识与特点MATLABMathworks1984要适用于矩阵运算、信息处理和控制等领域的分析与设计。是“编程语言，到科学界和教育界的高度好评。其特点是：数值计算功能强，每个数值元素都视为复数，用双精度位）化；具有符号计算功能；简单而强大的作图功能；人机界面适合科技人员；功能多，可扩展性强，有多种工具箱作专门处理。MATLAB的窗口与应用常识命令窗口以MATLAB6.x为例，其集成环境有桌面平台及8个组成部分：指令窗口、历史指令窗口、工作台及工具箱窗口、工作目录窗口、工作空间窗口、矩阵编辑器、程序编辑器和帮助浏览器。点击桌面上的MTLAB图标进入MATLABCommandWindo,它是主窗口，当显示“〉”时，表示系统处于准备接受命令的状态。图形窗口一般，只要执行任一种绘图命令，就会自动产生图形窗口,以后的绘图都在这一个图形窗口中进行。文本编辑窗口另一种是对复杂计算时，把多行命令组成一个M自动执行。M文件MATLAB.m“_”。Help命令Help命令是查询函数相关信息的最基本方式，信息会直接显示在命令窗口中。应用常识①MATLAB使用双精度数据，所有系统命令都是小写形式。进行数据运算的基本元素，矩阵中的下标从1开始，而不是01×1的矩阵来处理的。③语句或命令的结尾的分号“；”会屏蔽当前结果的显示。④注释位于%之后，不被执行。⑤使用上下箭头实现命令的滚动显示，可用于再编辑和再执行。⑥变量名必须以字母开头，可以由字母、数字和下划线组成。⑦要画出一条平滑的连续曲线，最少需要200个数据点。⑧即使预先对t做了定义，也应该使用数组操作，如f(t)=(2.^(-t)).*sin(2*t);若写作2^(-t)*sin(2*t)时会提示出错。⑨要想在同一个图形坐标上绘制多条曲线，应使用plot(t,f)，holdon，plot(t1,f1)，holdoff或plot(t,f,t1,f1)均可。⑩为了在图形中标出，，等特殊符号，应采用MATLAB提供的专门字母表示。如表1所示。表1 特殊符号的表示法字符串\alpha\beta符号字符串\theta\lambda符号字符串\phi\Phi符号\gamma\tao\pi\delta\int\infty\epsilon\omega\OmegaMATLAB的常用函数基本数学函数常用的基本数学函数如表2所示。表2 基本数学函数正弦sin平方根sqrt余弦cos求数组的最大值max正切tan求数组的最小值min指数exp求数组的平均值mean自然对数log求和sum矩阵指数函数expm求矩阵的特征值eig绝对值abs求复数的相角angle求实部real求共轭conj求虚部imag辛格函数sinc基本作图函数3表3常用的作图函数plot绘制连续波形title为图形加标题stem绘制离散波形grid画网格线loglog绘双对数坐标图xlabel为X轴加上轴标subplot分割图形窗口ylabel为Y轴加上轴标hold保留当前曲线text在图上加文字说明figure定义图形窗口gtext用鼠标在图上加说明line画直线axiszzplotx,y画符号函数的图形4．矩阵的创建在MATLAB中，矩阵是运算的基本单元，所以创建矩阵非常重要。创建矩阵有三种方法：直接输入法、内部函数创建法和利用矩阵编辑器和修改矩阵。直接输入法实数矩阵：对于维数较小的实数矩阵，可以从键盘直接输入。具体方法是：将矩阵的元素用方括号“[]”括起来，按行输入各元素，各行用分号“；”或Enter隔开。1 2 3

4 5 6 7 8 9 可写成指令：a=[123;456;789]复数矩阵：矩阵的各元素为复数时，与上述实数矩阵类似。不过写法有规定，如6+4j不可写成6+j4。应该写为6+j*4或6+i*4。内部函数创建法利用MATLAB提供的许多函数可以创建一些特殊矩阵。例如ones(m,n) 生成一个元素全部为1的m行n列矩阵eye(m,n) 生成一个主对角线元素为1的其余为0的矩pascal(n) 生成一个n维帕斯卡矩阵向量的创建创建向量最常用的方法是利用冒号“：”运算符生成。例如：a=m:p:nmnp例如：>>a=1:-0.1:0则输出结果为：a=1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0信号波形的产生和运算函数这里列出常用的波形产生和运算的函数如下：sawtooth产生周期锯齿波或三角波cumsum信号累加square产生周期方波sum信号求和pulstran脉冲串diff信号微分sinc辛格函数quad信号积分fliplr信号翻转conv信号卷积连续系统时域分析的函数在时域分析中，常用的符号形式如表4所示表4 时域分析符号名称阶跃函数(t)冲激函数(t)定义符号变量化简微分定积分求一阶微分方程全解求二阶微分方程全解

符号形式Heaviside(t)Dirac(t)symssimple(f)diff(f)diff(f,2)int(f,a,b)dsolve(‘Dy+a*y=b’,’y(0)=k’)dsolve(‘D2y+Dy……’)

说明使之包含的字符最少f被积函数，a,b为积分的上下限a,b,k均为数值向量求冲激响应

impulse(b,a)

a[a0

,a, ,a1

],b[b,b1 2

, ,b]M求阶跃响应

step(b,a) 同上连续信号频域分析函数进行周期信号的频谱分析，其符号形式为：int(f,t,a,b)fbAe，则有：f=A*exp(-j*n*2*pi/T*t);Fn=int(f,t,-tao/2,tao/2)/T;求非周期信号的傅里叶变换时，调用格式为：F=fourier(f,t,w)ff(t)令格式为：f=ifourier(F,w,t)f=ifourier(F)例如，对f(t)3e2t(t)，则有f=3*exp(-2*t)*sym(‘Heaviside(t)’);F=fourier(f);分析系统的频率响应，调用格式为：H=fregs(b,a,w)H(w)H(w)的取样点数。系统的复频域分析函数①信号的拉氏变换指令为：

Fs=laplace(ft,t,s)式中，ft为信号f(t)的符号表达式，t为积分变量，s为复频率，Fs为拉氏变换F(s)。②信号的拉氏反变换格式为：ft=ilaplace(fs,s,t)或者 ft=ilaplace(fs)若利用部分分式展开求反变换，则展开式的指令格式：[r,p,k]=residue(b,a)[b,b,,bF(s)[aa,a

]为分母多项式的0 1 m 0 1 n系数向量，r[r,r,,r]T为留数列向量，即展开的各分式的系数；p为各极点列向量，12 nk为直接项系数行向量。③由H(S)确定零、极点，指令格式为：sys=tf(b,a) 指定系统函数sys=zpk(z,p,k) 由零、极点构成系统函z=tzero(sys),p=pole(sys) 求系统的零、极点pzmap(sys) 绘制零、极点图impulse(sys) H(s)绘制冲激响应曲线step(sys) H(s)绘制阶跃响应曲线离散系统分析的基本函数①求差分方程零状态响应时，其调用格式为：y=filter(b,a,x)其中，xa[a0

,a, ,a1

]为差分方程输出变量的系数向量，b[b,b,,b

]为输入变量的系数向量。0 1 M求单位响应h(n)是调用格式为：h=impz(b,a)②求离散卷积和时，调用格式为：y=conv(h,x)其中，h为单位响应序列值，x为输入序列值，分别存储h(n)和x(n)的行向量。③z变换，调用格式为：Fz=ztrans(fn,n,z)其中，fnf(n)z反变换的指令为：impz(b,a,N)，N为样点的长度也可以用部分分式展开法，即[r,p,k]=residue(b,a)④计算频率响应，调用格式为：freqz(b,a)freqz(b,a,n)其中，baH(z)z的整数次幂；绘制频率特性的横坐标（即T）0到。二、MATLAB方法实现信号波形1方法实现单位阶跃信号和矩形脉冲。stairs(t和矩形脉冲的程序如下：t=-1:2;x=(t>=0);%定义时间范围向量tsubplot(1,2,1),stairs(t,x);axis([-1,2,-0.1,1.2]);gridon%绘制单位阶跃信号波形t=-1:0.001:1;%定义时间范围向量tg=(t>=(-1/2))-(t>=(1/2));subplot(1,2,2),stairs(t,g);axis([-1,1,-0.1,1.2]);运行结果如图1所示。gridon%绘制矩形脉冲波形图1 例1图例2应用MATLAB方法生成信号f(t)sinc(t)和f(t)S(t)的波形。a解：为生成函数sinc(t)sintt可直接调用MATLAB中的专门命令，程序如下：t=-5:0.01:5; %tf=sinc(t); %Sa(t)函数plot(t,f);gridon %Sa(t)的波形运行结果如图2所示。图2 例2程序运行结果一S(t和sinc(t的关系如下：at) f(t)S

(t)sint

t)sinc(t)a t

t t生成信号f(t)Sa

程序如下：t=-3*pi:0.01*pi:3*pi; %定义时间范围向量f=sinc(t/pi); %计算Sa(t)函数plot(t,f);gridon %Sa(t)的波形运行结果如图3所示。图3 例2程序运行结果二例3应用方法生成相加信号f(tcos和相乘信号f(t)sinc(t)cos(20t)的波形。解：对相加信号f(t)cos18tcos20t，程序如下：symst; %tf=cos(18*pi*t)+cos(20*pi*t); %计算符号函数f(t)=cos(18*pi*t)+cos(20*pi*t)ezplot(f,[0pi]);gridon %绘制f(t)的波形运行结果如图4所示。图4 例3程序运行结果一对相乘信号f(t)sinc(t)cos(20t)，程序如下：t=-5:0.01:5; %定义时间范围向量f=sinc(t).*cos(20*pi*t); %f(t)=sinc(t)*cos(20*pi*t)plot(t,f); %f(t)的波形title('sinc(t)*cos(20*pi*t)');gridon %加注波形标题运行结果如图5所示。图5 例3程序运行结果二4f(t(22sincos的波形。解：f(t(22sincos，程序如下：symst; %tf=(2+2*sin(4*pi*t))*cos(50*pi*t); %计算符号函数f(t)=(2+2*sin(4*pi*t))*cos(50*pi*t)ezplot(f,[0pi]);gridon %绘制f(t)的波形运行结果如图6所示。图图6 例4图三、MATLAB方法用于时域分析5设方程y(t5y(t6y(t2et(ty(t。解：程序如下：yzs=dsolve('D2y+5*Dy+6*y=2*exp(-t)','y(0)=0,Dy(0)=0') %dslove命令求解零状态响应ezplot(yzs,[08]);gridon %绘制零状态响应曲线运行结果：即：7

yzs=exp(-t)+exp(-3*t)-2*exp(-2*t)y(t)(et2e2te3t)(t)图7 例5图例6已知二阶系统方程u(t)Ru(t)

1u(t)1(t)c Lc

LCc LC对下列情况分别求h(t)，并画出其波形。a.R4,L,C1/3F b.R2,L,Cc.RL,Cd.R0,L,C解：程序如下：R=input('电阻R=');%以交互方式输入电阻R的值L=input('电感L=');%以交互方式输入电阻L的值C=input('电容C=');%以交互方式输入电阻C的值b=[1/(L*C)];a=[1R/L1/(L*C)];impulse(b,a); %绘制脉冲响应h(t)的波形运行结果如图8至图11所示。R=4电感L=1C=1/3图8 例6程序运行结果一R=2L=1C=1图9 例6程序运行结果二R=1电感L=1C=1图10 例6程序运行结果三R=0L=1C=1图11 例6程序运行结果四7f(t*h(tf(t(t(th(t(t(t2。解：主程序如下：p=0.01; %取样时间间隔nf=0:p:1; f(t)对应的时间向量f=2*((nf>=0)-(nf>=1)); %f(n)的值nh=0:p:2; h(t)对应的时间向量h=(nh>=0)-(nh>=2); %h(n)的值[y,k]=sconv(f,h,nf,nh,p); %计算y(t)=f(t)*h(t)subplot(3,1,1),stairs(nf,f);gridon %绘制f(t)的波形title('f(t)');axis([0302.1]);subplot(3,1,2),stairs(nh,h);gridon %h(t)title('h(t)');axis([0302.1]);subplot(3,1,3),plot(k,y);gridon y(t)=f(t)*h(t)title('y(t)=f(t)*h(t)');axis([0302.1]);子程序sconv如下：%此函数用于计算连续信号的卷积y(t)=f(t)*h(t)function[y,k]=sconv(f,h,nf,nh,p)%y:卷积积分y(t)对应的非零样值向量%k:y(t)对应的时间向量%f:f(t)对应的非零样值向量%nf:f(t)对应的时间向量%h:h(t)对应的非零样值向量%nh:h(t)对应的时间向量%p:取样时间间隔y=conv(f,h); %f(n)h(n)y(n)y=y*p; %y(n)y(t)left=nf(1)+nh(1); %计算序列y(n)非零样值的起点位置right=length(nf)+length(nh)-2; %计算序列y(n)非零样值的终点位置k=p*(left:right); %确定卷积和y(n)非零样值的时间向运行结果如图12所示。图12 例7图例8实现卷积f(t)*h(t)，其中：f(t)2[(t)(t2)],h(t)et(t)。解：主程序如下：p=0.01; %取样时间间隔nf=0:p:2; f(t)对应的时间向量f=2*((nf>=0)-(nf>=2)); %f(n)的值nh=0:p:4; h(t)对应的时间向量h=exp(-nh); %h(n)的值[y,k]=sconv(f,h,nf,nh,p); %计算y(t)=f(t)*h(t)subplot(3,1,1),stairs(nf,f);gridon %绘制f(t)的波形title('f(t)');axis([0602.1]);subplot(3,1,2),plot(nh,h);gridon %h(t)title('h(t)');axis([0602.1]);subplot(3,1,3),plot(k,y);gridon y(t)=f(t)*h(t)title('y(t)=f(t)*h(t)');axis([0602.1]);子程序sconv同例7。运行结果如图13所示。图13 例8图四、MATLAB方法用于频域分析914(a)MATLABN=7N=21时的合成图。解：该信号的系数:F0.5Sa(n)n 2前N项的合成表达式为：ft)0.5NNn1

Sancosn为奇数2利用MATLAB工具分析的程序如下：t=-3:0.001:3; %tN=input('N='); %NF0=0.5;fN=F0*ones(1,length(t));forn=1:2:NfN=fN+cos(pi*n*t)*sinc(n/2);endplot(t,fN); %fN的波形title(['N='num2str(N)]);axis([-33-0.21.2]);gridon运行结果如图14(b)和图14(c)所示。f(t)11…………-2.5-0.500.51.522.5 t(a)(a)(a)(b) (c)图14 周期矩形脉冲的合成例10如图15所示周期矩形脉冲f(t)，试绘出其频谱图。1……1…………t-4.5-4-3.5 -0.500.5 3.544.5图15 例10图解：程序如下：clearallsymstnTtaoT=4;A=1;tao=1;f=A*exp(-j*n*2*pi/T*t);fn=int(f,t,-tao/2,tao/2)/T; %计算傅里叶系数fn=simple(fn); %化简n=[-20:-1,eps,1:20]; %给定频谱的整数自变量,eps0fn=subs(fn,n,'n'); %n的值subplot(2,1,1),stem(n,fn,'filled'); %绘制频谱line([-2020],[00]); %在图形中添加坐标title('周期矩形脉冲的频谱');gridonsubplot(2,1,2),stem(n,abs(fn),'filled'); %绘制频谱title('周期矩形脉冲的幅度谱');gridonaxis([-202000.3]);运行结果如图15所示。图15 例10频谱图1116(a)f(t11f(t)1

,2t2，试求其频谱F()。t2tt-2 0 2(b)图16 例11图解：程序如下：symstwfft; %定义符号变量f=(1-(abs(t)/2)); %三角波信号ft=f*exp(-j*w*t); %计算被积函数F=int(ft,t,-2,2); %F(w)F=simple(F);F %化简axis([-3301.1]);title('三角波信号');ezplot(abs(F),[-8:0.01:8]);gridon %绘制三角波信号的频谱title('三角波信号的频谱');运行结果：F=-(cos(2*w)-1)/w^2即：1cos(2) 2sin2()F()2S2()2频谱如图16(b)所示。

2 a例12二阶低通滤波器特性为H()12

11jQ00 00即：1

1 1 Q 1211212 0 0Q 2令Q

121时，分别求幅频特性和相频特性。2

1 0解：程序如下：Q=input('输入Q='); %以交互方式输入normalizedw=linspace(0.1,10,100);H=1./(1-normalizedw.^2+j*normalizedw/Q); %二阶低通滤波器的频率特性表达式subplot(1,2,1),plot(normalizedw,abs(H)); %绘制幅频特性曲线title('幅频特性曲线');gridsubplot(1,2,2),plot(normalizedw,angle(H)); %绘制相频特性曲线title('相频特性曲线');grid1718Q=1/sqrt(2)Q=1

图17 例12程序运行结果一图图18 例12程序运行结果二13三阶低通滤波器特性为H()1(j)33(j)22(j)1试求：a.该系统的幅频特性H)和相频特性()，b.该系统的冲激响应h(t)。解：求幅频特性H)和相频特性()的程序如下：w=0:0.01:5;H=1./((j*w).^3+3*(j*w).^2+2*j*w+1); %三阶低通滤波器的频率特性表达式subplot(1,2,1),plot(w,abs(H)); %绘制幅频特性曲线title('幅频特性曲线');grid;axistight;subplot(1,2,2),plot(w,angle(H)); %绘制相频特性曲title('相频特性曲线');grid;axistight;图19 图19 例13程序运行结果一求该系统的冲激响应h(t)的程序如下：b=[1];%分子多项式系数a=[1321]; %分母多项式系数impulse(b,a); %h(t)运行结果如图20所示。图20 例13程序运行结果二例14应用MATLAB方法生成信号ft)St)pt)pt)的波形如图21所示。a1……1…………t-1.4 -1-0.8-0.4-0.200.20.40.81 1.4图21 例14图解：程序如下：t=-3*pi:0.01:3*pi; %定义时间范围向量s=sinc(t/pi); %计算Sa(t)函数subplot(3,1,1),plot(t,s);gridon %绘制Sa(t)的波形p=zeros(1,length(t)); %预定义p(t)的初始值为fori=16:-1:-16p=p+rectpuls(t+0.6*i,0.4); %rectpuls0.4,1p(t)endsubplot(3,1,2),stairs(t,p); %用阶梯图形表示矩形脉冲axis([-101001.2]);gridf=s.*p;subplot(3,1,3),plot(t,f);gridon %f(t)=Sa(t)*p(t)的波形运 行 结 果 如 图 22 所 示 。图22 例14程序运行结果1523f(tf(ty(ty(t的频谱1 2图。f(t)Hf(t)H)低通4y(t)1c(t)T(a)(t)T…………11-0.5 0 0.5 f(t)12|t|

-Ts 0 Ts Ts=0.2s(c)图23 例15图解：程序如下：symstwf; %定义符号变量f=(1-2*abs(t))*exp(-j*w*t); %计算被积函数F=int(f,t,-1/2,1/2); %F(w)F=simple(F);F %化简subplot(3,1,1), %F(w)low=-26*pi;high=-low; %wezplot(abs(F),[low:0.01:high]);gridonaxis([lowhigh-0.10.5]);xlabel('\omega');title('三角波的频谱');subplot(3,1,2), %4*pi低通滤波器后的频谱Y1(w)ezplot(abs(F),[-4*pi:0.01:4*pi]);gridonaxis([lowhigh-0.10.5]);xlabel('\omega');title('低通滤波后的频谱');%取样信号的频谱是原信号频谱的周期延拓，延拓周期为(2*pi)/Ts%利用频移特性F[f(t)*exp(-j*w0*t)]=F(w+w0)来实现subplot(3,1,3); %Y(w)Ts=0.2; %取样信号的周期w0=(2*pi)/Ts; %10*pifork=-2:2ft=f*exp(-j*w0*k*t);FT=int(ft,t,-1/2,1/2);ezplot((1/Ts)*abs(FT),[(-4*pi-k*w0):0.01:(4*pi-k*w0)]);holdonendaxis([lowhigh-0.12.5]);xlabel('\omega');gridtitle('取样后的频谱');运行结果如图24所示。图24 例15各频谱图五、MATLAB方法用于复频域分析2s1例16 部分分式展开：F(s) s32s25s解：程序如下：b=[21]; %分子多项式系数a=[1250]; %分母多项式系数[rpk]=residue(b,a) %部分分式展开运行结果：r=-0.1000-0.4500i-0.1000+0.4500i0.2000p=-1.0000+2.0000i-1.0000-2.0000i0k=[

0.10.45j 0.10.45j 0.2故F(s) s(12j) s(12j) s例17求拉氏变换：a.f(t)etcost b.f(t)3e2t(t)解：程序如下：symstw %twfat=exp(-t)*cos(w*t);fbt=3*exp(-2*t);fas=laplace(fat) %拉氏变换fbs=laplace(fbt) %拉氏变换运行结果：fas=(s+1)/((s+1)^2+w^2)fbs=3/(s+2)即：s1 3F(s)，F(s) a (s1)22 b s2例18求拉氏反变换：2s1 s2a.F(s)b.F(s) s27s10 s23s2解：程序如下：symss %s为符号变量fas=(2*s+1)/(s^2+7*s+10);fbs=s^2/(s^2+3*s+2);fat=ilaplace(fas) %拉氏反变换fbt=ilaplace(fbs) %拉氏反变换运行结果：fat=3*exp(-5*t)-exp(-2*t)fbt=Dirac(t)-4*exp(-2*t)+exp(-t)即：f(t)(3e5te2t(t) fa

(t)(t)(4e2tet)(t)例19零极点分析s2H(s)，求零极点并画出零极点图，并求阶跃响应 s(t)和冲击响应s24s5h(t)。解：程序如下：b=[12]; %系统函数分子多项式系数a=[145]; %系统函数分母多项式系数sys=tf(b,a); %传递函数subplot(1,3,1),pzmap(sys); %绘制零极点图subplot(1,3,2),step(b,a); %阶跃响应s(t)subplot(1,3,3),impulse(b,a); %冲激响应h(t)运行结果如图25所示。图25 例19程序运行结果一注：将鼠标移到零极点上即能显示其位置坐标。H(s) s2 H(s的零极点分布。s33s22s1程序如下：b=[12]; %系统函数分子多项式系数a=[1321]; %系统函数分母多项式系数sys=tf(b,a); %H(s)pzmap(sys); %绘制零极点图运行结果如图26所示。图图26 例19程序运行结果二例20设有系统函数

H(s)1s32s22s1试画出零、极点分布图，并求该系统的冲激响应和频率特性H()。解：程序如下：num=[1]; %系统函数分子多项式系数den=[1221]; %系统函数分母多项式系数sys=tf(num,den);poles=roots(den)subplot(2,2,1),pzmap(sys);title('零、极点分布图'); %绘制零、极点分布图t=0:0.02:10;h=impulse(num,den,t);subplot(2,2,2),plot(t,h);title('冲激响应');gridon %绘制冲激响应曲线[Hw]=freqs(num,den);subplot(2,2,3),plot(w,abs(H));title('幅频特性');gridon %绘制幅频特性曲subplot(2,2,4),plot(w,angle(H));title('相频特性');gridon %绘制相频特性曲线运行结果为：poles=-1.0000-0.5000+0.8660i-0.5000-0.8660i各图形如图27所示。图图27 例20图2128所示，F(s)+F(s)+.G(s)_其系统函数：

H(s)

图28 例21图34.5s2119.7s98.1s335.714s2119.741s98.1试求其阶跃响应s(t)。解：程序如下：num=[34.5119.798.1]; %系统函数分子多项式系数den=[135.714119.74198.1]; %系统函数分母多项式系数step(num,den);gridon %s(t)运行结果如图29所示。图29 例21的阶跃响应22某卫星角度跟踪天线控制系统的系统函数为：H(s)

1375020s4174s32268s213400s13750试画出其零极点图，并求其冲激响应h(t)。解：程序如下：num=[13750]; %系统函数分母多项式系数den=[2017422681340013750]; %系统函数分子多项式系数sys=tf(num,den); %传递函数H(s)subplot(1,2,1),pzmap(sys); %绘制零极点图subplot(1,2,2),impulse(num,den);gridon %h(t)运行结果如图30所示。图30 例22图六、MATLAB方法用于离散系统时域分析23已知差分方程y(ny(n0.8y(n2f(n)f(n0.5n(ny(n；f(n(n时，求单位响应h(n。解：程序如下：b=[1];a=[1-10.8];%差分方程的系数n=0:15; %序列的个数fn=0.5.^n; %输入序列y1=filter(b,a,fn); %零状态响y2=impz(b,a,16); %单位响应subplot(1,2,1),stem(n,y1,'filled');title('零状态响应');gridonsubplot(1,2,2),stem(n,y2,'filled');title('单位响应');gridon运行结果如图31所示。图31 例23图例24求卷积和：若f(n)0.8n5(n5)，到n30，h(n)R10

y(nf(n)*h(n。解：程序如下：nf=5:30;Nf=length(nf); %f(n)的序号向量和区间长度f=0.8.^(nf-5); %f(n)序列值nh=0:9;Nh=length(nh); %h(n)的序号向量和区间长度h=ones(1,Nh);; %h(n)序列值left=nf(1)+nh(1); %确定卷积序列的起点right=nf(Nf)+nh(Nh); %确定卷积序列的终点y=conv(f,h); %f(n)h(n)的卷积subplot(3,1,1),stem(nf,f,'filled');gridon %f(n)axis([04001]);title('f(n)');subplot(3,1,2),stem(nh,h,'filled');gridon %h(n)axis([04001.1]);title('h(n)');subplot(3,1,3),stem(left:right,y,'filled');gridon %y(n)axis([04005]);title('y(n)=f(n)*h(n)');运行结果如图32所示。图32 例24图25数学模型为：

{0,1,1,2,3,5,8,13,……}y(n)y(n1)y(n2)0试求n0~30时y(n)的值，并分别画出n0~20和n0~30曲线y(n)。解：程序如下：y(1)=0;y(2)=1;fori=3:30y(i)=y(i-1)+y(i-2);enddisp(y);subplot(1,2,1),stem(y(1:20),'filled');gridonn=0~20subplot(1,2,2),stem(y(1:30),'filled');gridonn=0~30时菲波那契数列的图形运行结果为：0112358132134558914423337761098715972584418167651094617711286574636875025121393196418317811514229图形如图33所示。图图33 例25图七、MATLAB方法用于z域分析26z解：程序如下：

f(n)cosnsymsnzOMEGA %定义符号变量fn=cos(n*OMEGA); %定义f(n)fz=ztrans(fn,n,z) %对f(n)进行z变运行结果：fz=(z-cos(OMEGA))*z/(z^2-2*z*cos(OMEGA)+1)即：z(zcos)F(z) z22zcos127解：程序如下：

z