1、信号与系统第5版MATLAB仿真举例MATLAB仿真举例一、MATLAB应用概要1. MATLAB的知识与特点 MATLAB是由美国Mathworks公司于1984年推出的一种功能强大的科学计算软件，主要适用于矩阵运算、信息处理和控制等领域的分析与设计。MATLAB是“矩阵实验室”（MATrix LABoratoy）的缩写。它以矩阵运算为基础，以简洁便用的“科学便笺式”编程语言，可交互的集成环境，强大的图形可视化功能，广泛应用于科学、工程和绘图的各种领域，受到科学界和教育界的高度好评。其特点是： (1) 数值计算功能强，每个数值元素都视为复数，用双精度（64位）一种格式，编程简化； (2) 具
2、有符号计算功能；(3) 简单而强大的作图功能； (4) 人机界面适合科技人员；(5) 功能多，可扩展性强，有多种工具箱作专门处理。2 MATLAB的窗口与应用常识(1) 命令窗口以MATLAB6.x为例，其集成环境有桌面平台及8个组成部分：指令窗口、历史指令窗口、工作台及工具箱窗口、工作目录窗口、工作空间窗口、矩阵编辑器、程序编辑器和帮助浏览器。点击桌面上的MATLAB图标进入MATLAB后，即可看到命令窗口（Command Window）,它是主窗口，当显示“”时，表示系统处于准备接受命令的状态。(2) 图形窗口一般，只要执行任一种绘图命令，就会自动产生图形窗口,以后的绘图都在这一个图形窗口
3、中进行。(3) 文本编辑窗口用MATLAB计算有两种方式，一种是简单计算，可在命令窗口一行一行地输入各命令，另一种是对复杂计算时，把多行命令组成一个M文件，让MATLAB自动执行。(4) M文件采用MATLAB所规定的一套语言及语法编写的源文件记作.m文件。其文件名不能以数字开头，也不能包含文字，但可以用下划线“_”。(5) Help命令Help命令是查询函数相关信息的最基本方式，信息会直接显示在命令窗口中。(6) 应用常识 MATLAB使用双精度数据，所有系统命令都是小写形式。 矩阵是MATLAB进行数据运算的基本元素，矩阵中的下标从1开始，而不是0，标量是作为11的矩阵来处理的。 语句或命
4、令的结尾的分号“；”会屏蔽当前结果的显示。 注释位于%之后，不被执行。 使用上下箭头实现命令的滚动显示，可用于再编辑和再执行。 变量名必须以字母开头，可以由字母、数字和下划线组成。 要画出一条平滑的连续曲线，最少需要200个数据点。 即使预先对t做了定义，也应该使用数组操作，如f(t)=(2.(-t).*sin(2*t); 若写作2(-t)*sin(2*t)时会提示出错。 要想在同一个图形坐标上绘制多条曲线，应使用plot(t,f)，hold on，plot(t1,f1)，hold off或plot(t,f,t1,f1)均可。 为了在图形中标出，等特殊符号，应采用MATLAB提供的专门字母表示
5、。如表1所示。表1 特殊符号的表示法字符串 符号字符串 符号字符串 符号alpha beta gamma delta epsilon theta lambda tao int omega phi Phi pi infty Omega 3 MATLAB的常用函数(1) 基本数学函数常用的基本数学函数如表2所示。表2 基本数学函数正弦 sin余弦 cos正切 tan指数 exp自然对数 log矩阵指数函数 expm绝对值 abs求实部 real求虚部 imag 平方根 sqrt 求数组的最大值 max 求数组的最小值 min 求数组的平均值 mean 求和 sum 求矩阵的特征值 eig 求复数的
6、相角 angle 求共轭 conj 辛格函数 sinc (2) 基本作图函数作图函数如表3所示表3 常用的作图函数 plot 绘制连续波形 stem 绘制离散波形 loglog 绘双对数坐标图 subplot 分割图形窗口 hold 保留当前曲线 figure 定义图形窗口 line 画直线 title 为图形加标题 grid 画网格线 xlabel 为X轴加上轴标 ylabel 为Y轴加上轴标 text 在图上加文字说明 gtext 用鼠标在图上加说明 axis 定义x,y坐标轴标度 zzplot 画符号函数的图形4 矩阵的创建在MATLAB中，矩阵是运算的基本单元，所以创建矩阵非常重要。创
7、建矩阵有三种方法：直接输入法、内部函数创建法和利用矩阵编辑器和修改矩阵。(1) 直接输入法实数矩阵：对于维数较小的实数矩阵，可以从键盘直接输入。具体方法是：将矩阵的元素用方括号“”括起来，按行输入各元素，各行用分号“；”或Enter隔开。例如矩阵a= 可写成指令：a=1 2 3; 4 5 6; 7 8 9复数矩阵：矩阵的各元素为复数时，与上述实数矩阵类似。不过写法有规定，如6+4j不可写成6+j4。应该写为6+j*4或6+i*4。(2) 内部函数创建法利用MATLAB提供的许多函数可以创建一些特殊矩阵。例如： ones(m,n) 生成一个元素全部为1的m行n列矩阵 eye(m,n) 生成一个主
8、对角线元素为1的其余为0的矩阵 pascal(n) 生成一个n维帕斯卡矩阵(3) 向量的创建创建向量最常用的方法是利用冒号“：”运算符生成。例如：a=m:p:n是生成起始值为m，终值为n，步长为p的均匀等分行向量a。例如：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 05 信号波形的产生和运算函数这里列出常用的波形产生和运算的函数如下： sawtooth 产生周期锯齿波或三角波 cumsum 信号累加 square 产生周期方波 sum 信号求和 pulstran 脉冲串 diff 信号微分 sinc 辛格函数 quad 信
9、号积分 fliplr 信号翻转 conv 信号卷积 6 连续系统时域分析的函数在时域分析中，常用的符号形式如表4所示表4 时域分析符号名称符号形式说明阶跃函数冲激函数定义符号变量化简微分二阶导数定积分求一阶微分方程全解求二阶微分方程全解求冲激响应求阶跃响应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)impulse(b,a)step(b,a)使之包含的字符最少f被积函数，a,b为积分的上下限a,b,k均为数值向量同上7 连续信号频域分析函数进行周期信
10、号的频谱分析，其符号形式为：int(f,t,a,b)其中f为符号表达式，t为积分变量，a和b分别为积分的下限和上限例如，计算傅里叶系数的被积函数，则有：f=A*exp(-j*n*2*pi/T*t);Fn=int(f,t,-tao/2,tao/2)/T;求非周期信号的傅里叶变换时，调用格式为：F=fourier(f,t,w)其中f为信号表达式，t为积分变量，w为角频率，F表示f(t)的傅氏变换。傅里叶反变换指令格式为：f=ifourier(F,w,t)f=ifourier(F)例如，对，则有f=3*exp(-2*t)*sym(Heaviside(t);F=fourier(f);分析系统的频率响应
11、，调用格式为：H=fregs(b,a,w)式中，b为H(w)中的分子多项式的系数向量，a为分母多项式系数向量，w为计算H(w)的取样点数。8系统的复频域分析函数 信号的拉氏变换指令为：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)其中，为F(s)的分子多项式系数向量；为分母多项式的系数向量，为留数列向量，即展开的各分式的系数；p为各极
12、点列向量，k为直接项系数行向量。 由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)绘制阶跃响应曲线9 离散系统分析的基本函数 求差分方程零状态响应时，其调用格式为：y=filter(b,a,x)其中，x为输入信号行向量，；为差分方程输出变量的系数向量，为输入变量的系数向量。求单位响应h(n)是调用格式为：h=impz(b,a)
13、求离散卷积和时，调用格式为：y=conv(h,x)其中，h为单位响应序列值，x为输入序列值，分别存储h(n)和x(n)的行向量。 z变换，调用格式为：Fz=ztrans(fn,n,z)其中，fn为f(n)的符号表达式；n为序号；z为复频率。用长除法进行z反变换的指令为：impz(b,a,N)， N为样点的长度也可以用部分分式展开法，即r,p,k=residue(b,a) 计算频率响应，调用格式为：freqz(b,a) freqz(b,a,n)其中，b和a分别为系统函数H(z)的分子、分母之系数向量，n为频率的计算点数，常取z的整数次幂；绘制频率特性的横坐标（即T）的范围为0到。二、MATLAB
14、方法实现信号波形例1 应用MATLAB方法实现单位阶跃信号和矩形脉冲。解：对于阶跃函数，MATLAB中有专门的stairs绘图命令。例如，实现和矩形脉冲的程序如下：t=-1:2; % 定义时间范围向量tx=(t=0); subplot(1,2,1),stairs(t,x);axis(-1,2,-0.1,1.2); grid on % 绘制单位阶跃信号波形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); grid on % 绘制矩形脉冲波形运行结果如图1所示
15、。图1 例1图例2 应用MATLAB方法生成信号和的波形。解：为生成函数可直接调用MATLAB中的专门命令，程序如下：t=-5:0.01:5; % 定义时间范围向量tf=sinc(t); % 计算Sa(t)函数plot(t,f); grid on % 绘制Sa(t)的波形运行结果如图2所示。图2 例2程序运行结果一和的关系如下：生成信号波形的MATLAB程序如下：t=-3*pi:0.01*pi:3*pi; % 定义时间范围向量tf=sinc(t/pi); % 计算Sa(t)函数plot(t,f); grid on % 绘制Sa(t)的波形运行结果如图3所示。图3 例2程序运行结果二例3 应用M
16、ATLAB方法生成相加信号和相乘信号的波形。解：对相加信号，程序如下：syms t; % 定义符号变量tf=cos(18*pi*t)+cos(20*pi*t); % 计算符号函数f(t)=cos(18*pi*t)+cos(20*pi*t)ezplot(f,0 pi); grid on % 绘制f(t)的波形运行结果如图4所示。图4 例3程序运行结果一对相乘信号，程序如下：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
17、)*cos(20*pi*t); grid on % 加注波形标题运行结果如图5所示。图5 例3程序运行结果二例4 应用MATLAB方法生成调制信号的波形。解：对调制信号，程序如下：syms t; % 定义符号变量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,0 pi); grid on % 绘制f(t)的波形运行结果如图6所示。图6 例4图三、MATLAB方法用于时域分析例5设方程，试求零状态响应。解：程序如下：yzs=dsolve(D2y+5*Dy+6*y=2*exp(
18、-t),y(0)=0,Dy(0)=0) % 利用dslove命令求解零状态响应ezplot(yzs,0 8); grid on % 绘制零状态响应曲线运行结果：yzs =exp(-t)+exp(-3*t)-2*exp(-2*t)即：图形如图7所示。图7 例5图例6 已知二阶系统方程对下列情况分别求，并画出其波形。a. b. c. d. 解：程序如下：R=input(电阻R=); % 以交互方式输入电阻R的值L=input(电感L=); % 以交互方式输入电阻L的值C=input(电容C=); % 以交互方式输入电阻C的值b=1/(L*C);a=1 R/L 1/(L*C);impulse(b,a
19、); % 绘制脉冲响应h(t)的波形运行结果如图8至图11所示。a. 电阻R=4 电感L=1 电容C=1/3图8 例6程序运行结果一b. 电阻R=2 电感L=1 电容C=1图9 例6程序运行结果二c. 电阻R=1 电感L=1 电容C=1图10 例6程序运行结果三d. 电阻R=0 电感L=1 电容C=1图11 例6程序运行结果四例7 实现卷积，其中：。解：主程序如下：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
20、)的值y,k=sconv(f,h,nf,nh,p); % 计算y(t)=f(t)*h(t)subplot(3,1,1),stairs(nf,f); grid on % 绘制f(t)的波形 title(f(t);axis(0 3 0 2.1);subplot(3,1,2),stairs(nh,h); grid on % 绘制h(t)的波形title(h(t);axis(0 3 0 2.1);subplot(3,1,3),plot(k,y); grid on % 绘制y(t)=f(t)*h(t)的波形title(y(t)=f(t)*h(t);axis(0 3 0 2.1);子程序sconv如下：%
21、 此函数用于计算连续信号的卷积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(n
22、h)-2; % 计算序列y(n)非零样值的终点位置k=p*(left:right); % 确定卷积和y(n)非零样值的时间向量运行结果如图12所示。图12 例7图例8 实现卷积，其中：。解：主程序如下：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); grid on % 绘制f(
23、t)的波形 title(f(t);axis(0 6 0 2.1);subplot(3,1,2),plot(nh,h); grid on % 绘制h(t)的波形title(h(t);axis(0 6 0 2.1);subplot(3,1,3),plot(k,y); grid on % 绘制y(t)=f(t)*h(t)的波形title(y(t)=f(t)*h(t);axis(0 6 0 2.1);子程序sconv同例7。运行结果如图13所示。图13 例8图四、MATLAB方法用于频域分析例9 求图14(a)所示周期矩形脉冲信号的傅里叶级数表示式，并用MATLAB方法求出N=7和N=21时的合成图。
24、解：该信号的系数:前N项的合成表达式为： 为奇数利用MATLAB工具分析的程序如下：t=-3:0.001:3; % 定义时间范围向量tN=input(N=); % 以交互方式输入N的值F0=0.5;fN=F0*ones(1,length(t);for n=1:2:N fN=fN+cos(pi*n*t)*sinc(n/2);endplot(t,fN); % 绘制fN的波形title(N= num2str(N);axis(-3 3 -0.2 1.2); grid on运行结果如图14(b)和图14(c)所示。 (a)(a) (b) (c)图14 周期矩形脉冲的合成例10 如图15所示周期矩形脉冲，
25、试绘出其频谱图。 图15 例10图解：程序如下：clear allsyms t n T tao A T=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; % 给定频谱的整数自变量,eps代表0fn=subs(fn,n,n); % 计算傅里叶系数对应各个n的值subplot(2,1,1),stem(n,fn,filled); % 绘制频谱line(-20 20,0 0); % 在图形中添加坐标线title(周期矩形脉冲的频谱);
26、 grid onsubplot(2,1,2),stem(n,abs(fn),filled); % 绘制频谱title(周期矩形脉冲的幅度谱); grid onaxis(-20 20 0 0.3);运行结果如图15所示。图15 例10频谱图例11 如图16(a)所示三角波信号，即：，试求其频谱。 (a) (b)图16 例11图解：程序如下：syms t w f ft; % 定义符号变量f=(1-(abs(t)/2); % 三角波信号ft=f*exp(-j*w*t); % 计算被积函数 F=int(ft,t,-2,2); % 计算傅里叶变换F(w)F=simple(F);F % 化简axis(-3
27、 3 0 1.1);title(三角波信号);ezplot(abs(F),-8:0.01:8); grid on % 绘制三角波信号的频谱title(三角波信号的频谱);运行结果：F =-(cos(2*w)-1)/w2即：频谱如图16(b)所示。例12 二阶低通滤波器特性为即：和令和1时，分别求幅频特性和相频特性。解：程序如下：Q=input(输入Q=); % 以交互方式输入Qnormalizedw=linspace(0.1,10,100);H=1./(1-normalizedw.2+j*normalizedw/Q); % 二阶低通滤波器的频率特性表达式subplot(1,2,1),plot(
28、normalizedw,abs(H); % 绘制幅频特性曲线title(幅频特性曲线);gridsubplot(1,2,2),plot(normalizedw,angle(H); % 绘制相频特性曲线title(相频特性曲线);grid运行结果如图17和图18所示。输入Q=1/sqrt(2)图17 例12程序运行结果一输入Q=1图18 例12程序运行结果二例13三阶低通滤波器特性为试求：a. 该系统的幅频特性和相频特性，b. 该系统的冲激响应。解：求幅频特性和相频特性的程序如下：w=0:0.01:5; H=1./(j*w).3+3*(j*w).2+2*j*w+1); % 三阶低通滤波器的频率特
29、性表达式subplot(1,2,1),plot(w,abs(H); % 绘制幅频特性曲线title(幅频特性曲线);grid;axis tight;subplot(1,2,2),plot(w,angle(H); % 绘制相频特性曲线title(相频特性曲线);grid;axis tight;运行结果如图19所示。图19 例13程序运行结果一求该系统的冲激响应的程序如下：b=1; % 分子多项式系数a=1 3 2 1; % 分母多项式系数impulse(b,a); % 冲激响应h(t)运行结果如图20所示。图20 例13程序运行结果二例14 应用MATLAB方法生成信号，其中的波形如图21所示。 图21 例14图解：程序如下：t=-
