1、目 录预备篇 MATLAB 基础及其在信号与系统中应用简介1 MATLAB 简介 .12 MATLAB 应用开发环境 13 MATLAB 基本运算 .44 MATLAB 的可视化绘图 .85 常用信号的 MATLAB 表示 13实验篇 信号与系统实验指导实验一、MATLAB 编程基础及典型实例 .20实验二、连续时间系统的时域分析 .23实验三、连续时间系统的频域分析 .27实验四、连续时间系统的复频域分析 .36实验五、离散时间系统的时域、Z 域分析 .45实验六、连续信号的采样与恢复 .571预备篇 MATLAB 基础及其在信号与系统中应用简介本篇简介：本篇简单介绍了 MATLAB 软件应
2、用基础 及其在信号与系统中经常用到的一些知识。由于同学们已经学习过 MATLAB 软件课程，所以这里只是提纲携领地提供 MATLAB 的相关内容以供复习；若有不懂之 处， 还请课下自己查看有关书籍予以解决。注意：本篇内容是后面实验部分的基础，非常重要，要 务必牢固地加以掌握。1 MATLAB 简介在科学研究与工程应用过程中，人们往往会遇到大量繁重的数学运算和数值分析，传统的高级语言如 FORTRAN、C 等虽然能够在一定的程度上减轻运算量，但它们均要求应用人员具有较强的编程能力和对算法有深入的研究。此外，对大多数科学工作者而言，若要运用这些高级语言对计算结果进行可视化分析以及对计算的图形进行处
3、理，也不是一件轻松的事情。MATLAB 正是在这一应用要求背景下产生的数学类科技应用软件。它具有强大的数值计算和图形可视化功能、简洁易学的工作环境和编程语言，从根本上满足了高校师生和科技人员对工程计算的要求，将他们从繁重的数学运算中解放出来，现已风靡世界，受到越来越多科技人员的喜爱和欢迎。2 MATLAB 应用开发环境2.1 命令窗口点击桌面上的 MATLAB 图标，或点击“开始程序 MATLAB 组中的MATLAB 程序项” ，即可运行 MATLAB，出现的界面就是命令窗口（Command Window） ，其提示符为“ ”。 在命令窗口中，可以直接输入所编写的命令，然后按回车键运行。例如输
4、入：x1 = sqrt(9), x2 = 5/x1, 按回车键后，命令窗口显示：x1 = 3 2x2 = 1.66672.2 M 文件对于一些比较简单的问题，可从命令窗口中输入直接输入指令加以执行。但随着指令数的增加、控制流复杂度的增加，以及重复计算的要求，直接在命令窗口中进行计算就显得很烦琐。MATLAB 提供的 M 文件很容易地解决了这个问题。M 文件实际上是由一系列 MATLAB 指令所构成的集合或程序文件，以“文件名.m”形式存放。文件名以英文字母开头，一般只能由英文字母、数字以及下划线符号“_”组成。在本讲义中所有的示范程序都是以 M 文件的形式给出；同学们在做实验时所编写的程序也必
5、须以 M 文件形式存放。M 文件有两种形式：脚本文件（Script file）和函数文件（ Function file） 。下面通过例子说明。示例 1：脚本文件 Ex_1x1 = 5;x2 = 6;x3 = 7;average_x = (x1+x2+x3)/3;运行 Ex_1 后，average_x 等于 6。示例 2：函数文件 Ex_2function y = aver(x1,x2,x3)y = (x1+x2+x3)/3;然后在主程序所在的脚本文件中或直接在命令窗口中执行：average_x = Ex_2(5,6,7), 或者先给 x1,x2,x3 赋值，再 average_x = Ex_2
6、(x1,x2,x3)，得到同样的结果。M 文件的创建与保存： 在命令窗口中点击工具栏“New M-File”按钮或者菜单“FileNew M-File”出现一个 M 文件编辑窗口（ Editor）,即创建了一个新的空白 M 文件； 输入相应的内容后，点击编辑窗口（Editor）工具栏中的 “Save”按钮或者菜单“FileSave” 即可保存。3M 文件的运行： 在命令窗口中点击“FileOpen”,选择 M 文件并打开； 在该 M 文件的编辑窗口中点击菜单“DebugRun” ；或直接在命令窗口中键入该 M 文件的文件名（如 Ex_1，不要带.m 后缀） ，按回车，即可运行该 M 文件。2.
7、3 MATLAB 的文件管理在编辑和运行 M 文件之前，还有一个很重要的工作要做，就是设置MATLAB 的当前工作目录和搜索路径。MATLAB 中的目录概念实际上就是文件夹，所有 M 文件都存放在某一目录中。当前工作目录的设置MATLAB 的使用者首先要新建一个目录（既文件夹）用来专门存放自己所编制的 M 文件，包括脚本文件和函数文件；然后，再将该目录设置成当前工作目录，方法如下：在命令窗口中右侧点击工具栏“Browse for folder”按钮，出现“浏览文件夹”选项，选中待确定为当前工作目录的文件夹，按“确定”即可。将自己的文件夹设置成当前工作目录后，就可以在 MATLAB 中运行自己所
8、编制的 M 文件了；也可以调用在搜索路径目录中所列出的所有文件夹中的M 函数文件，这些 M 函数文件可以是 MATLAB 中其它工具箱中自带的，也可以是你的同学编制的。这样，就达到了资源共享的目的。注意：你自己编制的 M 函数文件若不放在当前工作目 录中，而单独另外放在一个文件夹中，只要将此文件 夹设置为搜索路径， 则 你同样可调用这些 M 函数文件。搜索路径的设置在命令窗口中点击菜单“FileSet Path”，在出现的界面左上角点击“Add Folder” 按钮，出现“浏览文件夹”选项，选择待加入搜索路径的的目录，按“确定” ，再“Save ”、 “Close”即可。2.4 其它通过 MA
9、TLAB 的工作空间浏览器可以直观地查看 MATLAB Workspace 中包含的所有元素，这对于检查程序的运行状态、调试程序等是非常方便的。在4命令窗口中点击菜单“DesktopWorkspace ”，即出现 Workspace 界面，可以在其中查看内存变量的取值、维数等，也可以对变量做保存、删除等操作。MATLAB 还提供了强大的帮助系统。例如，点击工具栏 “？”按钮，出现“help”界面，点击其中的“Search”按钮，在“ Search for”框中输入所要查寻的关键字后回车，就能得到详细的帮助信息。在命令窗口中点击菜单“HelpDemos” ，再根据出现的界面做相应的选择，就可进入
10、 MATLAB 自带的各种演示程序。这些演示程序对初学者是一个很好的学习工具，可以方便地在不同条件下完成算法的仿真，并显示形象的可视化结果。3 MATLAB 基本运算3.1 基本概念变量：和其他高级语言一样，MATLAB 也是使用变量来保存信息。变量名由英文字母开头，一般只能由英文字母、数字以及下划线“_”组成。矩阵：矩阵是（Matrix）是 MATLAB 进行数据处理的基本单元，矩阵也经常与数组（Array ）不加区分。MATLAB 中的大部分运算或命令都是在矩阵运算的意义下执行的，矩阵运算是 MATLAB 最重要的运算。通常意义上的数量（也称为标量）在 MATLAB 系统中是作为 11 的
11、矩阵来处理的。 矩阵的输入在命令窗口或 M 文件编辑窗口中直接输入 a = 1,2,3; 4,5,6; 7,8,9 或 a = 1 2 3; 4 5 6; 7 8 9 ，就可生成 33 矩阵 a。 利用 MATLAB 函数创建矩阵MATLAB 为用户提供了创建基本矩阵的函数，如ones(m, n )：产生 mn 全 1 矩阵；zeros(m, n)：产生 mn 全 0 矩阵。向量：向量实际也是一种矩阵，是仅有一行或者一列的矩阵；它在基于 MATLAB5的信号与系统分析中发挥着重要作用。除了利用前面介绍的创建矩阵的方法来生成向量外，下面再介绍两种常用的方法。 利用冒号“：”运算生成向量如： x
12、= -2 : 4，则 x = -2, -1, 0, 1, 2, 3, 4;又如：y = 0 : 0.2: 1, 则 y = 0, 0.2, 0.4, 0.6, 1。要特别注意的是向量中元素的序号是从 1 开始的，例如上面的 x = -2, -1, 0, 1, 2, 3, 4 中，x(1) = -2, x(2) = -1, x(3) = 0 等；在 y = 0, 0.2, 0.4, 0.6, 1 中，y(1) = 0, y(2) = 0.2 等。 利用 linspace( )函数生成向量linspace( )函数用于生成线性等分向量。调用格式 x = linspace(m, n, s)表示生成从
13、起始值 m 开始到终止值 n 之间的 s 个线性等分点的行向量。例如 x = linspace(0, 10, 5)，则 x = 0, 2.5, 5, 7.5 10。3.2 矩阵的算术运算首先要说明的是在这里我们将矩阵视为数组，涉及运算的两个矩阵维数相同。在这个条件下，两个矩阵的加、减、乘、除均指的是两个矩阵相对应位置上的元素进行加、减、乘、除运算。例如设 A = , B = , 则在15926643.MATLAB 中四种运算的表示及结果为：A+B ，AB ，A.*B ，A./B .215069. 95783. 9051348. 5263注意：上面乘法运算符号“*”前要加“.” ,即这种矩阵乘法
14、是一种点乘运算,它不同于线性代数中两个矩阵之间的所定义的乘法；在线性代数中要求 A 的列数与 B 的行数相同才能进行乘法运算，而 这里仅 要求 A、B 的维数相同。点乘运算是 MATLAB 所特有的，极大地 简化了编程，要予以高度重视。同理，上面例子中的除法用的也是点除。6MATLAB 还提供了点幂运算“.” ，如X = ，则 X3 为6543212165478另外，一个矩阵还可与一个数进行加、减、乘、除运算，其结果是该矩阵中的每一个元素与这个数进行相应的运算。如：A = ，则 A10 ，A*2 15926512451243.3 关系运算MATLAB 的基本关系运算符为：（大于） ， =（大于
15、等于） ，= Y ，X = Y 01001Y n2)|(n1n2)error(输入不正确,输入参数要应满足 n1=”,就可得到单位阶跃函数 、单位冲激序列 ，如图 6 所示。(t)(n)-5 0 500.51似似似似似似nu(n)-5 0 500.51似似似似似似tu(t)图 619我们编制了函数文件 stepseq.m 来生成单位阶跃序列 。)n-u(0function x, n = stepseq(n1,n2,n0)% 产生序列 u(n-n0),其中 n1n2)|(n1n2)error(输入不正确 ,输入参数要应满足 n1=0);示例 8：绘图表示(1) 门函数 ；1t0,(t)g2(2)
16、 序列 。5)u(n3xnt = -3 :0.05: 3;z1 = (t+1) = 0);z2 = (t-1) = 0);g = z1 - z2; % 门函数figure;subplot(221)plot(t,g,r);axis(-3 3 0 1.1)x1 = stepseq(-5,10,-3); % 调用函数 stepseqx2 = stepseq(-5,10,5); % 调用函数 stepseqx = x1 - x2;n = -5:10;subplot(222)stem(n,x);axis(-5,10,0,1.1)运行结果如图 7 所示。20-2 0 200.20.40.60.81-5 0
17、 5 1000.20.40.60.81图 75.3 其他典型的信号 实指数信号 nax()其 MATLAB 实现为：n = n1: n2; x = a.n; 复指数信号 )j(ep()其 MATLAB 实现为：n = n1: n2; x = exp(sigma+jw)*n; 正(余)弦信号 )cos(x)其 MATLAB 实现为：n = n1: n2; x = cos(w*n+sita)5.4 工具箱中的信号产生函数利用 MATLAB 信号处理工具箱提供的一些函数,可以很方便地产生三角波、方波等函数波形。 周期性三角波或锯齿波函数 sawtooth调用格式为：x = sawtooth(t, w
18、idth)功能：产生一个周期为 2、幅度在-1 到+1 之间的 周期性三角波信号。其中 width 表示最大幅度出现的位置：即在一个周期内，信号从 t=0到 width2 时函数值从-1 到+1 线性增加，而从 width2 到2 又是从+1 到-1 线性下降。 width 取值在 0 1 之间。示例 9：产生周期为 0.2 的三角波，width 取值分别为 0、1、0.5。td = 1/100000; % td 为时间间隔t = 0 : td : 1;21x1 = sawtooth(2*pi*5*t,0);x2 = sawtooth(2*pi*5*t,1);x3 = sawtooth(2*p
19、i*5*t,0.5);subplot(311); plot(t,x1);subplot(312); plot(t,x2);subplot(313); plot(t,x3);运行结果如图 8 所示。0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1010 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1010 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-101图 8 周期性方波信号 square调用格式为：x = square (t, duty)功能：产生一个周期为 2、幅度为1 的 周期性方波信号。其
20、中 duty 表示占空比，即在信号的一个周期中正值所占的百分比。例如产生频率为 40Hz、占空比为 75%的 周期性方波所调用的语句为x = square (2*pi*40*t, 75); (非周期)三角波脉冲信号 tripuls调用格式为：x = tripuls (t, width, skew)功能：产生一个最大幅度为 1、宽度为 width、斜率为 skew 的三角脉冲信号。该函数横坐标范围由向量 t 决定，其三角波形是以 t=0 为中心向22左右各展开 width/2 的范围；斜率 skew 在-1 到+1 之间取值，它决定了最大幅度 1 所对应的横坐标位置：width/ 2skew。示
21、例 10：仔细观察由下面代码产生的图 9 中 3 个三角波信号之间的区别，自己对 tripuls 函数的使用做一个总结。t = -3:0.001:3;x1 = tripuls(t,4,0);subplot(131); plot(t,x1);gridt = -6:0.001:6;x2 = tripuls(t,4,0.5);subplot(132);plot(t,x2);axis(-4 4 0 1); gridx3 = tripuls(t+2,4,0.5);subplot(133);plot(t,x3);axis(-4 4 0 1);grid-5 0 500.20.40.60.81-4 -2 0
22、2 400.20.40.60.81-4 -2 0 2 400.20.40.60.81图 9 (非周期)矩形脉冲信号 rectpuls调用格式为：x = rectpuls (t, width)功能：产生一个幅度为 1、宽度为 width、以 t=0 为中心左右对称的矩形波信号。该函数横坐标范围由向量 t 决定，其矩形波形是以 t=0 为中心向左右各展开 width/2 的范围。width 的默认值为 1。23示例 11：生成幅度为 2，宽度 T = 4、中心在 t = 0 的矩形波 x(t)以及 x(t-T/2).t = -4 : 0.0001 : 4;T = 4;x1 = 2*rectpuls
23、(t, T);subplot(121);plot(t, x1);axis(-4 6 0 2.2)grid;x2 = 2*rectpuls(t-T/2,T);subplot(122);plot(t, x2);axis(-4 6 0 2.2)grid;运行结果如图 10 所示。-4 -2 0 2 4 600.511.52-4 -2 0 2 4 600.511.52图 10 取样函数 sinc详见示例 5。24实验篇 信号与系统实验指导实验 一、 MATLAB 编程基础及典型实例一、实验目的(1) 熟悉 MATLAB 软件平台的使用；(2) 熟悉 MATLAB 编程方法及常用语句；(3) 掌握 MA
24、TLAB 的可视化绘图技术；(4) 结合信号与系统的特点，编程实现常用信号及其运算。二、实验原理连续信号是指自变量的取值范围是连续的，且对于一切自变量的取值，除了有若干个不连续点以外，信号都有确定的值与之对应。严格来说，MATLAB并不能处理连续信号，而是用等时间间隔点的样值来近似表示连续信号。当取样时间间隔足够小时，这些离散的样值就能较好地近似连续信号。矩阵是 MATLAB 进行数据处理的基本单元，矩阵运算是 MATLAB 最重要的运算。通常意义上的数量（也称为标量）在 MATLAB 系统中是作为 11 的矩阵来处理的，而向量实际上是仅有一行或者一列的矩阵。通常用向量表示信号的时间取值范围，
25、如 n = -5:5，但信号 x(n)、向量 n 本身的下标都是从 1 开始的，因此必须用一个与向量 x 等长的定位时间变量 n，以及向量 x，才能完整地表示序列 x(n)。这一点详情可参考预备篇示例 7 的程序说明。三、程序示例在预备篇中，我们对 MATLAB 应用开发环境、基本运算、 MATLAB 的可视化绘图方法、常用信号的 MATLAB 表示等做了清楚的介绍，并给出了很多具体实例，包括程序代码和产生的各种图形。请同学们认真参考这些内容，完成下面的相关实验。注意：在以后的实验中，对于以 t 为自变量的连续信号，在绘图时统一用 plot函数；而对 n 为自变量的离散序列，在绘图时统一用 s
26、tem 函数。25示例 1：在两个信号进行加、减、相乘等运算时，参于运算的两个向量要有相同的维数，并且它们的时间变量范围要相同，即要对齐。编制一个函数型 m 文件，实现这个功能。function f1_new, f2_new, n = duiqi(f1,n1,f2,n2)% 功能：将两个序列对齐，以实现两个序列之间的运算% 输入：% （1） f1,f2: 原来的两个序列；% （2） n1,n2: f1,f2 所对应的时间变量范围；% 输出：% f1_new, f2_new ：对齐后的两个序列% n : 对齐后的两个序列的时间变量范围%-a = min(min(n1), min(n2);b =
27、max(max(n1), max(n2);n = a : b;f1_new = zeros(1, length(n);f2_new = zeros(1, length(n);tem1 = find(n=min(n1)f2_new(tem2) = f2;26实验二、 连续时间系统的时域分析一、实验目的(1) 深刻理解卷积运算，掌握离散线性卷积、连续线性卷积的计算方法；(2) 加深对线性时不变系统中零状态响应概念的理解，掌握其求解方法；(3) 掌握给定连续系统的冲激响应和阶跃响应。二、实验原理(1)线性时不变 (LTI) 连续时间系统用常系数线性微分方程进行描述，系统的零状态响应就是在系统初始状态
28、为零条件下微分方程的解。MATLAB 控制系统工具箱提供了一个 lsim 函数来求解连续时间系统的零状态响应。设系统方程为：,f(t)btf(t)fb(t)fy(t)at(t)ya(t) )()( 01230123 该方程左边、右边的系数向量分别为 ， ，所对0123a, 123,应的系统模型 sys 可借助 MATLAB 中的 tf 函数得到：sys = tf(b, a) .这样，系统的零状态响应为：y = lsim(sys, f, t) ，其中 f 是输入信号向量，t 是与 f 对应的时间变量。(2)连续系统的冲激响应、阶跃响应分别是输入信号为 和 所对应的零)t(tu状态响应。MATLA
29、B 控制系统工具箱专门提供了两个函数求解连续系统的冲激响应和阶跃响应。冲激响应：y = impulse(sys, t) ;阶跃响应：y = step(sys, t) .其中 sys, t 的含义同上。(3)卷积是信号与系统中一个最基本、也是最重要的概念之一。在时域中，对于 LTI 连续时间系统，其零状态响应等于输入信号与系统冲激响应的卷积；而利用卷积定理，这种关系又对应频域中的乘积。如实验一所述，我们用离散卷积来代替连续卷积，只要取样时间间隔足够小时，就可得到满意的效果。27MATLAB 信号处理工具箱提供了一个计算两个离散序列卷积和的函数conv。设向量 a、b 代表待卷积的两个序列，则 c
30、 = conv(a, b)就是 a 与 b 卷积后得到的新序列。我们知道两个序列卷积以后，一般而言所得新序列的时间范围、序列长度都会发生变化。例如设 f1(n)长度为 5，3n1；f 2(n)长度为 7，2n8；则卷积后得到的新序列长度为 11，1n9。但是用 conv 函数求出卷积后没有给出新序列所对应的时间变量。为此，我们在下面的程序示例中给出了一个函数文件 dconv，它在完成 conv 函数功能的同时，还产生了一个对应新序列的时间变量。(4)对于连续卷积， )kt(f)kfd)t(f)t(ft)t(f 2102121 lim令 （ 为整数） ，则n(*) kk )kn(f)(f)n(f
31、)f)(f 2121由(*)式，连续卷积积分可由离散卷积和近似代替，只要取样时间间隔 足够小，就可以得到高精度卷积积分的数值计算。在示例 3 中给出了一个函数文件 cconv 来完成该功能。三、程序示例示例 1：已知系统的微分方程为， 。求零状态响应 。f(t)fy(t)(t)y)()( 3412 )t(ey(t)a = 1 4 4;b = 1 3;sys = tf(b, a);td = 0.01;t = 0 : td : 10;f = exp(-t);y = lsim(sys, f, t);plot(t, y);xlabel(t(sec);ylabel(y(t);grid on程序运行结果见
32、下图。280 1 2 3 4 5 6 7 8 9 1000.050.10.150.20.250.30.35t(sec)y(t)示例 2：利用 conv 函数，编制一个函数文件 dconv，其输出为两个序列卷积后的新序列以及与该新序列对应的时间变量。function f, k = dconv(f1, f2, k1, k2)% 计算 f1 与 f2 的卷积，并返回与得到的新序列相对应的时间变量f = conv(f1, f2);k_start = k1(1) + k2(1);k_end = length(f1) + length(f2) - 2;k = k_start : (k_start + k_
33、end);示例 3：在 dconv 函数和(*)式的基础上，编制一个函数文件 cconv，利用离散卷积和来近似计算连续卷积积分。function f, k = cconv(f1,f2,k1,k2,td)% 计算 f1 与 f2 的连续卷积，并返回与得到的新序列相对应的时间变量f = td*conv(f1,f2); % 实验讲义中的 (*)式计算k_start = k1(1) + k2(1);k_end = length(f1) + length(f2) - 2;k = k_start :td: (k_start + k_end*td);29四、实验内容与步骤(1) 已知系统的微分方程为， 。计
34、算系统的零状态响应 、(t)fyt()(t)y1122)t(fy(t)冲激响应 和阶跃响应 ，并画出相应的图形。tg(2) 用 MATLAB 计算如下两序列的卷积和，绘出它们的时域波形。其 它,k,)(f0121 其 它,k)k(f0212(3) 编程实现如下图所示的两个波形；并利用 cconv 函数计算这两个信号的卷积、画出卷积后的波形。t-1 0 12f1(t) -2 0 2t1f2(t)五、实验报告要求整理并给出“实验内容与步骤” （1） 、 （2） 、 （3）中的程序代码与产生的图形；并回答下面的问题。(1) 在“实验内容与步骤” （1） ，零状态响应 和阶跃响应 是否相同？为y(t)t(g什么？(2) 两序列进行卷积后得到新的序列，说明新序列在时域长度、时域区间上与原来两序列的关系。
