12-13-1信号与系统matlab实验指导书

实验一 基本信号的波形

一、实验目的：

1.掌握matlab软件的基本操作。 2.熟悉matlab的基本命令的使用。 3.掌握用matlab绘出基本信号。 二、实验原理：

1. 信号的时域表示方法

1.1将信号表示成独立时间变量的函数

例如 x(t)=sin(ωt) 和 x[n]=n(0.5)nu[n]

分别表示一个连续时间信号和一个离散时间信号。在MATLAB中有许多内部函数，可以直接完成信号的这种表达，例如：

sin()：正弦信号 cos()：余弦信号 exp()：指数信号

1.2用信号的波形图来描述信号

用函数曲线表示一个信号，图1.1就是一个连续时间信号和一个离散时间信号的波形图。

图1.1 连续时间信号与离散时间信号的波形图

1.3将信号用一个数据序列来表示

对于离散时间信号，还可以表示成一个数的序列，例如： x[n]={...., 0.1, 1.1, -1.2, 0, 1.3, ….} ↑n=0

在《信号与系统》和《数字信号处理》课程中，上述三种信号的描述方法是经常要使用的。

2 用MATLAB仿真连续时间信号和离散时间信号

在MATLAB中，无论是连续时间信号还是离散时间信号，MATLAB都是用一个数字

序列来表示信号，这个数字序列在MATLAB中叫做向量(vector)。通常的情况下，需要与时间变量相对应。

如前所述，MATLAB有很多内部数学函数可以用来产生这样的数字序列，例如sin()、

cos()、exp()等函数可以直接产生一个按照正弦、余弦或指数规律变化的数字序列。

2.1连续时间信号的仿真

程序Program1_1是用MATLAB对一个正弦信号进行仿真的程序，请仔细阅读该程

序，并在计算机上运行，观察所得图形。

% Program1_1

% This program is used to generate a sinusoidal signal and draw its plot clear, % Clear all variables

close all, % Close all figure windows

dt = 0.01; % Specify the step of time variable t = -2:dt:2; % Specify the interval of time x = sin(2*pi*t); % Generate the signal

plot(t,x) % Open a figure window and draw the plot of x(t) title('Sinusoidal signal x(t)') xlabel('Time t (sec)')

常用的图形控制函数

axis([xmin,xmax,ymin,ymax])：图型显示区域控制函数，其中xmin为横轴的显示起点，xmax为横轴的显示终点，ymin为纵轴的显示起点，ymax为纵轴的显示终点。

有时，为了使图形具有可读性，需要在所绘制的图形中，加上一些网格线来反映信号的幅度大小。MATLAB中的grid on/grid off可以实现在你的图形中加网格线。

grid on %在图形中加网格线。 grid off %取消图形中的网格线。

hold on %图形显示窗口原来的图像保持。 hold off %关闭图形保持功能。 figure %打开新的显示窗口。

x = input(‘Type in signal x(t) in closed form:’)%接收键盘输入。 在绘制信号的波形图时，有时我们需要将若干个图形绘制在图一个图形窗口中，这就需要使用MATLAB的图形分割函数subplot()，其用法是在绘图函数stem或plot之前，使用图形分割函数subplot(n1,n2,n3)，其中的参数n1，n2和n3的含义是，该函数将把一个图形窗口分割成n1×n2个子图，即将绘制的图形将绘制在第n3个子图中。

在《信号与系统》课程中，单位阶跃信号u(t) 和单位冲激信号δ(t) 是二个非常有用的信号。它们的定义如下

?t?????(t)dt?1t?0 1.1(a)

?(t)?0,?1,u(t)???0,t?0 1.1(b) t?0这里分别给出相应的简单的产生单位冲激信号和单位阶跃信号的扩展函数。产生单位冲激信号的扩展函数为：

function y = delta(t) dt = 0.01;

y = (u(t)-u(t-dt))/dt;

产生单位阶跃信号的扩展函数为： function y = u(t)

y = (t>=0); % y = 1 for t > 0, else y = 0

请将这二个MATLAB函数分别以delta 和u为文件名保存在work文件夹中，以后，就可以像教材中的方法使用单位冲激信号δ(t) 和单位阶跃信号u(t)。

单位冲激信号： x=-10:0.1:10;

y=(x==0); %就是个条件判断，只有x=0的时候，y才为“1” plot(x,y);

如果是单位脉冲信号（离散的），改成stem(x,y);就行了 t = -5:0.01:5; y = (t==0); subplot(121); plot(t, y, 'r'); n = -5:5; x = (n==0); subplot(122); >> stem(n, x);

x=-2:0.01:2; y=100*(x==0); stairs(x,y)

练习产生正弦信号（sin、cos）、指数信号（exp）、符号函数（sign）、抽样函数（sinc）、矩形脉冲（rectpuls(t,width)）、三角脉冲（tripuls(t,width,skew)）

2.2离散时间信号的仿真

程序Program1_2用来产生离散时间信号x[n]=sin(0.2πn)。

% Program1_2

% This program is used to generate a discrete-time sinusoidal signal and draw its plot clear, % Clear all variables

close all, % Close all figure windows n = -10:10; % Specify the interval of time x = sin(0.2*pi*n); % Generate the signal

stem (n,x) % Open a figure window and draw the plot of x[n] title ('Sinusoidal signal x[n]') xlabel ('Time index n')

请仔细阅读该程序，比较程序Program1_1和Program1_2中的不同之处，以便自己编程时能够正确使用这种方法方针连续时间信号和离散时间信号。 程序Program1_3用来仿真下面形式的离散时间信号： x[n]={...., 0.1, 1.1, -1.2, 0, 1.3, ….} ↑n=0

% Program1_3

% This program is used to generate a discrete-time sequence % and draw its plot

clear, % Clear all variables

close all, % Close all figure windows

n = -5:5; % Specify the interval of time, the number of points of n is

11.

x = [0, 0, 0, 0, 0.1, 1.1, -1.2, 0, 1.3, 0, 0]; % Generate the signal

stem(n,x,'.') % Open a figure window and draw the plot of x[n] grid on,

title ('A discrete-time sequence x[n]') xlabel ('Time index n')

由于在程序的stem(n,x,'.') 语句中加有'.'选项，因此绘制的图形中每根棒条线的顶端是一个实心点。

如果需要在序列的前后补较多的零的话，可以利用函数zeros()，其语法为： zeros(1, N)：圆括号中的1和N表示该函数将产生一个一行N列的矩阵，矩阵中的所有元素均为零。利用这个矩阵与序列x[n]进行组合，从而得到一个长度与n相等的向量。

例如，当 x[n]={ 0.1, 1.1, -1.2, 0, 1.3} 时，为了得到程序Program1_3中的序列， ↑n=0

可以用这个MATLAB语句x = [zeros(1,4) x zeros(1, 2)] 来实现。用这种方法编写的程序如下：

% Program1_4

% This program is used to generate a discrete-time sinusoidal signal % and draw its plot

clear, % Clear all variables

close all, % Close all figure windows n = -5:5; % Specify the interval of time

x = [zeros(1,4), 0.1, 1.1, -1.2, 0, 1.3, zeros(1,2)]; % Generate the sequence stem (n,x,'.') % Open a figure window and draw the plot of x[n] grid on,

title ('A discrete-time sequence x[n]') xlabel ('Time index n')

离散时间单位阶跃信号u[n]定义为 u[n]???1,?0,n?0 1.2 n?0离散时间单位阶跃信号u[n]除了也可以直接用前面给出的扩展函数来产生，还可以利用MATLAB内部函数ones(1,N) 来实现。这个函数类似于zeros(1,N)，所不同的是它产生的矩阵的所有元素都为1。

单位冲激序列： n=-10:10;

y=[zeros(1,10),1,zeros(1,10)]; stem(n,y)

单位阶跃序列： n=-10:10;

y=[zeros(1,10),ones(1,11)];

stem(n,y)

值得注意的是，利用ones(1,N) 来实现的单位阶跃序列并不是真正的单位阶跃序列，而是一个长度为N单位门(Gate)序列，也就是u[n]-u[n-N]。但是在一个有限的图形窗口中，我们看到的还是一个单位阶跃序列。

三、实验内容及步骤

Q1-1：修改程序Program1_1，将dt改为0.2，再执行该程序，保存图形，看看所得图形

的效果如何？

dt = 0.01时的信号波形 dt = 0.2时的信号波形

此处粘贴图形 此处粘贴图形

这两幅图形有什么区别，哪一幅图形看起来与实际信号波形更像？ 答：

Q1-2：修改程序Program1_1，并以Q1_2为文件名存盘，产生实指数信号x(t)=e-2t。 要求

在图形中加上网格线，并使用函数axis()控制图形的时间范围在0~2秒之间。然后执行该程序，保存所的图形。

修改Program1_1后得到的程序Q1_2如下： 信号x(t)=e-2t的波形图

Q1-3：将实验原理中所给的单位冲激信号和单位阶跃信号的函数文件在MATLAB文件编

辑器中编写好，并分别以以文件名delta和u存入work文件夹中以便于使用。

抄写函数文件delta如下： 抄写函数文件u如下：

Q1-4：修改程序Program1_4，并以Q1_4为文件名存盘，利用axis()函数，将图形窗口的

横坐标范围改为-2≤n≤5，纵坐标范围改为-1.5≤ x ≤1.5。

修改Program1_4后得到的程序Q1_4如下： 信号的波形图

此处粘贴图形

实验二 连续系统的时域分析

一、实验目的

1.掌握信号的时域变换。 2.信号的时域运算。 3.掌握卷积的计算。 4.掌握微分方程的求解。

二、实验原理

1信号的时域变换 1.1 信号的时移

信号的时移可用下面的数学表达式来描述：

设一个连续时间信号为x(t)，它的时移y(t) 表示为：

y(t) = x(t - t0) 2.1

其中，t0为位移量。若t0为正数，则y(t)等于将x(t)右移t0秒之后的结果。反之，若t0为负数，则y(t)等于将x(t)左移t0秒之后的结果。

在MATLAB中，时移运算与数学上习惯表达方法完全相同。

程序Program2_1对给定一个连续时间信号x(t) = e-0.5tu(t)，对它分别左移2秒钟和右移2秒钟得到信号x1(t) = e-0.5(t+2)u(t+2)和x2(t) = e-0.5(t-2)u(t-2)。

% Program2_1

% This program is used to implement the time-shift operation

% on a continuous-time signal and to obtain its time-shifted versions % and to draw their plots. clear,close all, t = -5:0.01:5;

x = exp(-0.5*t).*u(t); % Generate the original signal x(t)

x1 = exp(-0.5*(t+2)).*u(t+2); % Shift x(t) to the left by 2 second to get x1(t) x2 = exp(-0.5*(t-2)).*u(t-2); % Shift x(t) to the right by 2 second to get x2(t) subplot(311)

plot(t,x) % Plot x(t) grid on,

title ('Original signal x(t)') subplot (312)

plot (t,x1) % Plot x1(t) grid on,

title ('Left shifted version of x(t)') subplot (313)

plot (t,x2) % Plot x2(t) grid on,

title ('Right shifted version of x(t)') xlabel ('Time t (sec)')

1.2 信号的时域反褶

对一个信号x[n]的反褶运算在数学上表示为

y[n] = x[-n] 2.2

这种反褶运算，用MATLAB实现起来也是非常简单的。有多种方法可以实现信号的反褶运算。

方法一，修改绘图函数plot(t,x)和stem(n,x)中的时间变量t和n，即用-t和-n替代原来的t和n，这样绘制出来的图形，看起来就是原信号经时域反褶后的版本。

方法二，直接利用原信号与其反褶信号的数学关系式来实现。这种方法最符合信号反褶运算的实际意义。

1.3 信号的时域尺度变换

信号x(t)的时域尺度变换在数学描述为

y(t) = x(at)， 2.3

其中a为任意常数。根据a的不同取值，这种时域尺度变换对信号x(t)具有非常不同的影响。 当a = 1时，y(t) = x(t)；

当a = -1时，y(t) = x(-t)，即y(t)可以通过将x(t)反褶运算而得到； 当a > 1时，y(t) = x(at)，y(t)是将x(t)在时间轴上的压缩而得到； 当0 < a < 1时，y(t) = x(at)，y(t)是将x(t)在时间轴上的扩展而得到；

当 -1 < a < 0时，y(t) = x(at)，y(t)是将x(t)在时间轴上的扩展同时翻转而得到； 当 a < -1时，y(t) = x(at)，y(t)是将x(t)在时间轴上的压缩同时翻转而得到；

由此可见，信号的时域尺度变换，除了对信号进行时域压缩或扩展外，还可能包括对信号的时域反褶运算。实际上，MATLAB完成式2.3的运算，并不需要特殊的处理，按照数学上的常规方法即能完成。

1.4周期信号

在《信号与系统》课程中，周期信号是一类非常重要的信号。给定一个信号x(t)或x[n]，如果满足

x(t) = x(t+kT) 2.4 x[n] = x[n+kN] 2.5

则该信号叫做周期信号。其中，k为任意整数，T和N为常数，通常称为信号的基本周期或最小周期。

周期信号可以看作是一个时限的非周期信号经过周期延拓之后形成的。在数字信号处理中，周期延拓这一信号处理方法非常重要。

下面的程序段，就是将一个非周期信号x1(t) = e-2t[u(t)-u(t-2)]经过周期延拓之后而得到一个周期信号：

clear, close all; t = -4:0.001:4; T = 2; x = 0; for k = -2:2;

x = x+exp(-2*(t-k*T)).*(u(t-k*T)-u(t-(k+1)*T)); end plot(t,x)

仔细阅读该程序，可以发现其算法就是：

x(t)?

k????x(t?kT) 2.6

1?

由于k无法计算到无穷，而是以有限值加以替代，反映到有限宽度图形窗口中得到的效果完全符合要求。

2 信号的时域运算 2.1相加

syms x f1=sin(x) f2=x^2 f= f1+f2 ezplot(f)

?、f2是符号表达式。

2.2 相乘

syms x f1=sin(x) f2=x^2 f= f1*f2 ezplot(f)

?、f2是符号表达式。

2.3微分

syms x f=sin(x^2) h=diff(f) ezplot(h)

2.4积分

syms x f=sin(x)

h=int(f,x) %h=int(f,x,a,b)为定积分 ezplot(h)

3卷积的计算

卷积的计算通常可按下面的五个步骤进行（以卷积积分为例）：

1) 该换两个信号波形图中的横坐标，由t改为τ，τ变成函数的自变量； 2) 把其中一个信号反褶，如把h(τ)变成h(-τ)；

3) 把反褶后的信号做移位，移位量是t，这样t是一个参变量。 在τ坐标系中，t > 0时图形右移， t < 0时图形左移。

4) 计算两个信号重叠部分的乘积x(τ)h(t-τ)； 5) 完成相乘后图形的积分。

对于两个时限信号（Time-limited signal），按照上述的五个步骤，作卷积积分运算时，关键是正确确定不同情况下的积分限。只要正确地确定了积分限都能得到正确定积分结果。

尽管如此，在时域中计算卷积积分，总体上来说是一项比较困难的工作。

借助MATLAB的内部函数conv()可以很容易地完成两个信号的卷积积分运算。其语法为：y = conv(x,h)。其中x和h分别是两个作卷积运算的信号，y为卷积结果。

为了正确地运用这个函数计算卷积，这里有必要对conv(x,h)做一个详细说明。conv(x,h)函数实际上是完成两个多项式的乘法运算。例如，两个多项式p1和p2分别为：

p1?s3?2s2?3s?4 和 p2?4s3?3s2?2s?1

这两个多项式在MATLAB中是用它们的系数构成一个行向量来表示的，如果用x来表示多项式p1，h表示多项式p2，则x和h分别为 x = [1 2 3 4] h = [4 3 2 1] 在MATLAB命令窗口依次键入

>> x = [1 2 3 4]; >> h = [4 3 2 1]; >> y=conv(x,h)

在屏幕上得到显示结果：

y = 4 11 20 30 20 11 4 这表明，多项式p1和p2的乘积为：

p3?4s6?11s5?20s4?30s3?20s2?11s?4

正如前所述，用MATLAB处理连续时间信号时，独立时间变量t的变化步长应该是很小的，假定用符号dt表示时间变化步长，那么，用函数conv()作两个信号的卷积积分时，应该在这个函数之前乘以时间步长方能得到正确的结果。也就是说，正确的语句形式应为：y = dt*conv(x,h)。

对于定义在不同时间段的两个时限信号x(t)，t0 ≤ t ≤ t1，和h(t)，t2 ≤ t ≤ t3。 如果用y(t)来表示它们的卷积结果，则y(t)的持续时间范围要比x(t)或h(t)要长，其时间范围为t0+t2 ≤ t ≤ t1+t3。这个特点很重要，利用这个特点，在处理信号在时间上的位置时，可以很容易地将信号的函数值与时间轴的位置和长度关系保持一致性。

根据给定的两个连续时间信号x(t) = t[u(t)-u(t-1)]和h(t) = u(t)-u(t-1)，编写程序，完成这两个信号的卷积运算，并绘制它们的波形图。范例程序如下：

% Program2_2

% This program computes the convolution of two continuou-time signals clear;close all;

t0 = -2; t1 = 4; dt = 0.01; t = t0:dt:t1; x = u(t)-u(t-1); h = t.*(u(t)-u(t-1));

y = dt*conv(x,h); % Compute the convolution of x(t) and h(t) subplot(221)

plot(t,x), grid on, title('Signal x(t)'), axis([t0,t1,-0.2,1.2]) subplot(222)

plot(t,h), grid on, title('Signal h(t)'), axis([t0,t1,-0.2,1.2]) subplot(212)

t = 2*t0:dt:2*t1; % Again specify the time range to be suitable to the

% convolution of x and h.

plot(t,y), grid on, title('The convolution of x(t) and h(t)'), axis([2*t0,2*t1,-0.1,0.6]), xlabel('Time t sec')

在有些时候，做卷积和运算的两个序列中，可能有一个序列或者两个序列都非常长，甚至是无限长，MATLAB处理这样的序列时，总是把它看作是一个有限长序列，具体长度由编程者确定。实际上，在信号与系统分析中所遇到的无限长序列，通常都是满足绝对可和或绝对可积条件的信号。因此，对信号采取这种截短处理尽管存在误差，但是通过选择合理的信号长度，这种误差是能够减小到可以接受的程度的。若这样的一个无限长序列可以用一个数学表达式表示的话，那么，它的长度可以由编程者通过指定时间变量n的范围来确定。

例如，对于一个单边实指数序列x[n] = 0.5nu[n]，通过指定n的范围为0 ≤n ≤ 100，则对应的x[n]的长度为101点，虽然指定更宽的n的范围，x[n]将与实际情况更相符合，但是，注意到，当n大于某一数时，x[n]之值已经非常接近于0了。对于序列x[n] = 0.5nu[n]，当n = 7时，x[7] = 0.0078，这已经是非常小了。所以，对于这个单边实指数序列，指定更长的n的范围是没有必要的。当然，不同的无限长序列具有不同的特殊性，在指定n的范围时，只要能够反映序列的主要特征就可以了。

4 用线性常系数微分方程描述LTI系统

线性常系数微分方程或差分方程是描述LTI系统的另一个时域模型。一个连续时间LTI系统，它的输入信号x(t)输出信号y(t)关系可以用下面的微分方程来表达

dky(t)Mdkx(t) ?ak 2.7 ??bkkkdtdtk?0k?0N式2.7中，max (N, M)定义为系统的阶。式2.7描述了LTI系统输入信号和输出信号的一种

隐性关系（Implicit relationship）。为了求得系统响应信号的显式表达式（Explicit expression），必须对微分方程和差分方程求解。

在MATLAB中，一个LTI系统也可以用系统微分方程的系数来描述，例如，一个LTI连续时间系统的微分方程为

d2y(t)dy(t)?3?2y(t)?x(t)

dtdt2MATLAB则用两个系数向量num = [1]和den = [1 3 2]来描述该系统，其中num和den

分别表示系统微分方程右边和左边的系数，按照微分运算的降阶排列。

MATLAB的内部函数impulse()，step()，initial()，lsim() 可以用来计算并绘制连续时间LTI系统的单位冲激响应，单位阶跃响应，零输入响应和任意信号作用于系统的零状态响应。这些函数的用法描述如下：

h= impulse(num, den, T) 和 impulse(num, den, T) s = step(num, den, T) 和 step(num, den, T)

y = lsim(num, den, x, t) 和 lsim(num, den, x, t)

函数impulse()，step()用来计算由num和den表示的LTI系统的单位冲激响应和单位阶跃响应，响应的时间范围为0~T，其中den和num分别为系统微分方程左右两边的系数向量，T为指定的响应的终点时间。h和s的点数默认值为101点。由此可以计算时间步长为

dt = T/(101-1)。不带返回值的函数如impulse(num, den, T)和step(num, den, T)将直接在屏幕上绘制系统的单位冲激响和单位阶跃响应曲线。带返回值的函数如y = lmis(num, den, x, t)和y = lmis(num, den, x, t)，用来计算由num和den表示的LTI系统在输入信号x作用下的零状态响应。其中t为指定的时间变化范围，x为输入信号，它们的长度应该是相同的。如带返回参数y，则将计算的响应信号保存在y中，若不带返回参数y，则直接在屏幕上绘制输入信号x和响应信号y的波形图。

微分方程（课本42页，例题2.1-2）的全解计算： y=dsolve('D2y+5*Dy+6*y=10*cos(t)','y(0)=2,Dy(0)=0') 所得结果为：

y =2/exp(2*t) - 1/exp(3*t) + cos(t) + sin(t)

yht=dsolve('D2y+5*Dy+6*y=0') % 求齐次通解

yt=dsolve('D2y+5*Dy+6*y=10*cos(t)') % 求非齐次通解

yp=yt-yht % 求特解，即强迫响应 yh=y-yp % 求齐次解，即自由响应 求零输入响应和零状态响应

>> yzi=dsolve(' D2y+5*Dy+6*y =0',' y(0)=2,Dy(0)=0')

>> yzs=dsolve(' D2y+5*Dy+6*y=10*cos(t)','y(0)=0,Dy(0)=0')

例如，编写程序，计算并绘制由下面的微分方程表示的系统的单位冲激响应h(t)，单位阶跃响应s(t)。

d2y(t)dy(t)?3?2y(t)?8x(t) 2dtdtMATLAB范例程序如下：

% Program2_3

% This program is used to compute the impulse response h(t) and the step response s(t) of a

% continuous-time LTI system clear, close all;

num = input('Type in the right coefficient vector of differential equation：'); den = input('Type in the left coefficient vector of differential equation：'); t = 0:0.01:8;

subplot(211), impulse(num,den,8); subplot(212), step(num,den,8) 命令窗口为：

Type in the right coefficient vector of differential equation：[8] Type in the left coefficient vector of differential equation：[1,3,2] 三、实验内容及步骤

Q2-1：编写一个MATLAB程序，以Q2_1为文件名存盘，使之能够在同一个图形窗口中

的两个子图中分别绘制信号x(t)=0.5|t| 和x(t)=cos(2πt)[u(t)-u(t-3)]。要求选择的时间窗能够表

现出信号的主要部分（或特征）。

编写的程序Q2_1如下：

信号x(t)=0.5|t| 的波形图和信号x(t)=cos(2πt)[u(t)-u(t-3)]的波形图

Q2-2：根据示例程序的编程方法，编写一个MATLAB程序，以Q2_2为文件名存盘，由

给定信号

x(t) = e-0.5tu(t) 求信号y(t) = x(1.5t+3)，并绘制出x(t) 和y(t)的图形。

编写的程序Q2_2如下：

信号x(t)的波形图 和 信号y(t) = x(1.5t+3) 的波形图

此处粘贴图形 此处粘贴图形

Q2-3：给定一个离散时间信号x[n] = u[n] – u[n-8]，仿照示例程序Program2_1，编写程序

Q2_3，产生x[n]的左移序列x1[n] = x[n+6]和右移序列x2[n] = x[n-6]，并在同一个图形窗口的三个子图中分别绘制这三个序列的图形。

编写的程序Q2_3如下：

n=-10:15; x=u(n)-u(n-8); x1=u(n+6)-u(n-2); x2=u(n-6)-u(n-14); subplot(311) stem(n,x) subplot(312) stem(n,x1) subplot(313) stem(n,x2)

信号波形图

此处粘贴图形

Q2-4：编写程序Q2_4，使之能够接受以键盘方式输入的定义在不同时间段的两个不同连

续时间信号并完成卷积运算，分别绘制这两个信号及其卷积的结果的图形，图形按照2?2分割成四个子图。

编写的程序Q2_4如下：

t0=input('信号x1的开始时间：'); t1=input('信号x1的结束时间：'); t2=input('信号x2的结束时间：'); t3=input('信号x2的结束时间：'); dt=0.01; tx1=t0:dt:t1;

x1=exp(-0.5*tx1); tx2=t2:dt:t3;

x2=sin(2*pi*tx2); x=dt*conv(x1,x2); subplot(221) plot(tx1,x1) subplot(222) plot(tx2,x2)

tx=(t0+t2):dt:(t1+t3); subplot(212) plot(tx,x)

信号x (t)、h(t)和x (t)*h(t)的波形图

此处粘贴图形

Q2-5 仿照范例程序Program2_3，编写程序Q2_5，计算并绘制由如下微分方程表示的系统

在输入信号为x(t) = (e-2t - e-3t)u(t)时的零状态响应和你手工计算得到的系统零状态响应曲线。

d2y(t)dy(t)?3?2y(t)?8x(t) 2dtdt编写的程序Q2_5如下： 执行程序Q2_5得到的系统响应

此处粘帖执行程序Q2_5得到的系统响应

clear, close all; num = [8] den = [1 3 2]; t = 0:0.01:10;

x=(exp(-2*t)-exp(-3*t)).*u(t); y=lsim(num,den,x,t) plot(t,y)

实验三、连续系统的频域分析

一、实验目的

1、掌握连续时间周期信号的傅里叶级数的物理意义和分析方法； 2、观察截短傅里叶级数而产生的“Gibbs现象”，了解其特点以及产生的原因； 3、掌握连续时间傅里叶变换的分析方法及其物理意义；

4、掌握各种典型的连续时间非周期信号的频谱特征以及傅里叶变换的主要性质； 5、学习掌握利用MATLAB语言编写计算CTFS和CTFT的仿真程序，并能利用这些程序对一些典型信号进行频谱分析，验证CTFT的若干重要性质。

基本要求：掌握并深刻理傅里叶变换的物理意义，掌握信号的傅里叶变换的计算方法，掌握利用MATLAB编程完成相关的傅里叶变换的计算。

二、实验原理及方法

1、连续时间周期信号的傅里叶级数CTFS分析

任何一个周期为T1的正弦周期信号，只要满足狄利克利条件，就可以展开成傅里叶级数。

其中三角傅里叶级数为：

x(t)?a0??[akcos(k?0t)?bksin(k?0t)] 3.1

k?1??或： x(t)?a0??ck?1kcos(k?0t??k) 3.2

其中?0?2?，称为信号的基本频率（Fundamental frequency），a0，ak，和bk分别是信T1号x(t)的直流分量、余弦分量幅度和正弦分量幅度，ck、?k为合并同频率项之后各正弦谐波分量的幅度和初相位，它们都是频率k?0的函数，绘制出它们与k?0之间的图像，称为信号的频谱图（简称“频谱”），ck－k?0图像为幅度谱，?k－k?0图像为相位谱。 三角形式傅里叶级数表明，如果一个周期信号x(t)，满足狄里克利条件，那么，它就可以被看作是由很多不同频率的互为谐波关系（harmonically related）的正弦信号所组成，其中每一个不同频率的正弦信号称为正弦谐波分量 (Sinusoid component)，其幅度（amplitude）为ck。也可以反过来理解三角傅里叶级数：用无限多个正弦谐波分量可以合成一个任意的非正弦周期信号。

指数形式的傅里叶级数为：

x(t)?

k????aek?jk?0t 3.3

其中，ak为指数形式的傅里叶级数的系数，按如下公式计算：

1 ak?T1T1/2?T1/2?x(t)e?jk?0tdt 3.4

指数形式的傅里叶级数告诉我们，如果一个周期信号x(t)，满足狄里克利条件，那么，它就可以被看作是由很多不同频率的互为谐波关系（harmonically related）的周期复指数信号所组成，其中每一个不同频率的周期复指数信号称为基本频率分量，其复幅度（complex amplitude）为ak。这里“复幅度（complex amplitude）”指的是ak通常是复数。

上面的傅里叶级数的合成式说明，我们可以用无穷多个不同频率的周期复指数信号来合成任意一个周期信号。然而，用计算机（或任何其它设备）合成一个周期信号，显然不可能做到用无限多个谐波来合成，只能取这些有限个谐波分量来近似合成。

假设谐波项数为N，则上面的和成式为：

x(t)?k??N?aekNjk?0t 3.5

显然，N越大，所选项数越多，有限项级数合成的结果越逼近原信号x(t)。本实验可以比较直观地了解傅里叶级数的物理意义，并观察到级数中各频率分量对波形的影响包括“Gibbs”现象：即信号在不连续点附近存在一个幅度大约为9%的过冲，且所选谐波次数越多，过冲点越向不连续点靠近。这一现象在观察周期矩形波信号和周期锯齿波信号时可以看得很清楚。

2、连续时间信号傅里叶变换----CTFT

傅里叶变换在信号分析中具有非常重要的意义，它主要是用来进行信号的频谱分析的。

傅里叶变换和其逆变换定义如下：

? X(j?)????x(t)e??j?tdt 3.6

1x(t)?2???j?tX(j?)ed? 3.7 ? 连续时间傅里叶变换主要用来描述连续时间非周期信号的频谱。按照教材中的说法，任

意非周期信号，如果满足狄里克利条件，那么，它可以被看作是由无穷多个不同频率（这些频率都是非常的接近）的周期复指数信号ej?t的线性组合构成的，每个频率所对应的周期复指数信号ej?t称为频率分量（frequency component），其相对幅度为对应频率的|X(j?)|之值，其相位为对应频率的X(j?)的相位。

X(j?)通常为关于?的复函数，可以按照复数的极坐标表示方法表示为：

X(j?)=| X(j?)|ej? X(j?)

其中，| X(j?)|称为x(t)的幅度谱，而?X(j?)则称为x(t)的相位谱。

给定一个连续时间非周期信号x(t)，它的频谱也是连续且非周期的。对于连续时间周期信号，也可以用傅里变换来表示其频谱，其特点是，连续时间周期信号的傅里叶变换时有冲激序列构成的，是离散的——这是连续时间周期信号的傅里叶变换的基本特征。

3、连续周期信号的傅里叶级数CTFS的MATLAB实现

3.1 傅里叶级数的MATLAB计算

设周期信号x(t)的基本周期为T1，且满足狄里克利条件，则其傅里叶级数的系数可由式3.4计算得到。式3.4重写如下：

1ak?T1基本频率为： ?0?T1/2?T1/2?jk?0tx(t)edt ?2? T1对周期信号进行分析时，我们往往只需对其在一个周期内进行分析即可，通常选择主周期（Principle period）。假定x1(t)是x(t)中的主周期，则

1ak?T1T1/2?T1/2?jk?0tx(t)edt 1?计算机不能计算无穷多个系数，所以我们假设需要计算的谐波次数为N，则总的系数个

数为2N+1个。在确定了时间范围和时间变化的步长即T1和dt之后，对某一个系数，上述系数的积分公式可以近似为：

11ak?x1(t)e?jk?0tdt??x(tn)e?jk?0tdt/T1 ?T1?T1/2n ?[x(t1),x(t2),?x(tM)]?[e?jk?0t1T/2,e?jk?0t2,?e?jk?0tM]?dt/T1

对于全部需要的2N+1个系数，上面的计算可以按照矩阵运算实现。MATLAB实现系

数计算的程序如下：

dt = 0.01;

T = 2; t = -T/2:dt:T/2; w0 = 2*pi/T;

x1 = input('Type in the periodic signal x(t) over one period x1(t)='); %u(t)-u(t-0.3) N = input('Type in the number N='); 0 for k = -N:N;

ak(N+1+k) = x1*exp(-j*k*w0*t')*dt/T;%1×n矩阵与n×1矩阵相乘=一个数。 end k=-N:N; stem(k,ak)

需要强调的是，时间变量的变化步长dt的大小对傅里叶级数系数的计算精度的影响

非常大，dt越小，精度越高，但是，计算机计算所花的时间越长。

例题3-1：给定一个周期为T1 = 2s的连续时间周期方波信号，如图所示，其一个周期内

的数学表达式为：

x1(t)???1,?0,0?t?1

1?t?2-2 -1 1 x(t) 0 1 2 t 解：首先，我们根据前面所给出的公式，计算该信

图3.1 周期方波信号

号的傅里叶级数的系数。

?jk?t1e1?jk?t1100?jk?t?jk?t?ed(?jk?t)?edt?ak?x(t)edt01?j2k?0?2??j2k?0T1?T?/200T1/2100011?e?1?e?j2k?0?jk?0k?j?02ek?j?02?e?j2k?0kj?02ksin(?0)?jk?02?e2 k?0sin(因为：?0 = 2π/T1 = π，代入上式得到：ak?(?j)kk?k?)2

在MATLAB命令窗口，依次键入：

>> k = -10:10;

>> ak = ((-j).^k).* (sin((k+eps)*pi/2)./((k+eps)*pi)) % The expression of ak ak =

Columns 1 through 4

-0.0000 0 + 0.0354i -0.0000 0 + 0.0455i Columns 5 through 8

-0.0000 0 + 0.0637i -0.0000 0 + 0.1061i Columns 9 through 12

-0.0000 0 + 0.3183i 0.5000 0 - 0.3183i Columns 13 through 16

-0.0000 0 - 0.1061i -0.0000 0 - 0.0637i Columns 17 through 20

-0.0000 0 - 0.0455i -0.0000 0 - 0.0354i Column 21

-0.0000 从MATLAB命令窗口，我们得到了该周期信号从a?10到a10共21个系数。 紧接着再键入以下命令：

>> stem(k,abs(ak),'k.')

>> title('The Fourier series coefficients') >> xlabel('Frequency index k')

就得到一幅如右图所示的描述ak与k之间的关系的

图形。

以上是我们通过手工计算得到的这个周期信号的傅里叶级数表达式及其频谱图，下面给出完成傅里叶级数系数计算的相应MATLAB范例程序。

% Program3_1

% This program is used to evaluate the Fourier series coefficients ak of a periodic square

wave

clear, close all

T = 2; dt = 0.001; t = -2:dt:2; x1 = u(t) - u(t-1-dt); x = 0;

for m = -1:1 % Periodically extend x1(t) to form a periodic signal x = x + u(t-m*T) - u(t-1-m*T-dt); end

w0 = 2*pi/T;

N = 10; % The number of the harmonic components L = 2*N+1;

for k = -N: N; % Evaluate the Fourier series coefficients ak ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt; end

phi = angle(ak); % Evaluate the phase of ak

执行程序Program2_1后，就完成了信号的傅里叶级数的系数的计算，在命令窗口键入

>> ak

命令窗口就可以显示傅里叶级数的21个系数：

ak =

Columns 1 through 4

0.0000 + 0.0000i 0.0000 + 0.0354i 0.0000 - 0.0000i 0.0000 + 0.0455i Columns 5 through 8

0.0000 - 0.0000i 0.0000 + 0.0637i 0.0000 - 0.0000i 0.0000 + 0.1061i Columns 9 through 12

0.0000 - 0.0000i 0.0000 + 0.3183i 0.5000 0.0000 - 0.3183i Columns 13 through 16

0.0000 + 0.0000i 0.0000 - 0.1061i 0.0000 + 0.0000i 0.0000 - 0.0637i Columns 17 through 20

0.0000 + 0.0000i 0.0000 - 0.0455i 0.0000 + 0.0000i 0.0000 - 0.0354i Column 21

0.0000 - 0.0000i

将这里的ak之值同前面手工计算得到的ak比较，可见两者是完全相同的。

再次特别提示：程序中，时间变量的变化步长dt的大小对傅里叶级数系数的计算精度

的影响非常大，dt越小，精度越高，本程序中的dt之所以选择0.001就是为了提高计算精度。但是，计算机所花的计算时间越长。

在程序Program3_1中添加相应的计算| ak |和绘图语句，就可以绘制出信号的幅度谱和相位谱的谱线图。在程序最后加上：

m=-N:N;

stem(m,abs(ak)) figure

stem(m,phi)

3.2 周期信号的合成以及Gibbs现象

从傅里叶级数的合成式（Synthesis equation）

?x(t)?

k????aekjk?0t

可以看出，用无穷多个不同频率和不同振幅的周期复指数信号可以合成一个周期信号。然而，我们无法用计算机实现对无穷多个周期复指数信号的合成。但是，用有限项来合成却是可行的，在实际应用中，多半也就是这么做的。然而，这样做的一个必然结果，就是引入了误差。

如果一个周期信号在一个周期有内断点存在，那么，引入的误差将除了产生纹波之外，还将在断点处产生幅度大约为9%的过冲（Overshot），这种现象被称为吉伯斯现象（Gibbs phenomenon）。

为了能够观察到合成信号与原信号的不同以及Gibbs现象，我们可以利用前面已经计算出的傅里叶级数的系数，计算出截短的傅里叶级数：

x(t)?k??N?aekNjk?0t

这个计算可用L = 2N+1次循环来完成：

x2?x2?ak(r)?ej(r?1?N)?0t

其中r作为循环次数，x2在循环之前应先清零。完成这一计算的MATLAB程序为：

x2 = 0; L = 2*N+1; for r = 1:L;

x2 = x2+ak(r)*exp(j*(r-1-N)*w0*t); end;

完成了所有的计算之后，就可以用绘图函数：plot()和stem()将计算结果包括x, x2, abs(ak)和angle(ak)以图形的形式给出，便于我们观察。

plot(t,x)

axis([-2 2 -0.2 1.2]) figure plot(t,x2)

观察吉伯斯现象的最好的周期信号就是图3-1所示的周期方波信号，这种信号在一个周期内有两个断点，用有限项级数合成这个信号时，吉伯斯现象的特征非常明显，便于观察。

例题3-2：修改程序Program3_1，使之能够用有限项级数合成例题3-1所给的周期方波信

号，并绘制出原始周期信号、合成的周期信号、信号的幅度谱和相位谱。

为此，只要将前述的for循环程序段和绘图程序段添加到程序Program3_1中即可，范例程序如下：

% Program3_2

% This program is used to compute the Fourier series coefficients ak of a periodic square wave

clear,close all

T = 2; dt = 0.001; t = -2:dt:2; x1 = u(t)-u(t-1-dt); x = 0; for m = -1:1

x = x + u(t-m*T) - u(t-1-m*T-dt); % Periodically extend x1(t) to form a periodic signal

end

w0 = 2*pi/T;

N = input('Type in the number of the harmonic components N = ');

L = 2*N+1; for k = -N:1:N;

ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt; end

phi = angle(ak); y=0;

for q = 1:L; % Synthesis the periodic signal y(t) from the finite Fourier series

y = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T);% y = y+ak(q)*exp(j*(q-1-N)*2*pi*t/T); end;

subplot(221),

plot(t,x), title('The original signal x(t)'), axis([-2,2,-0.2,1.2]), subplot(223),

plot(t,y), title('The synthesis signal y(t)'), axis([-2,2,-0.2,1.2]), xlabel('Time t'), subplot(222)

k=-N:N; stem(k,abs(ak),'k.'), title('The amplitude |ak| of x(t)'), axis([-N,N,-0.1,0.6])

subplot(224)

stem(k,phi,'r.'), title('The phase phi(k) of x(t)'), axis([-N,N,-2,2]), xlabel('Index k') 在用这个程序观察吉伯斯现象时，可以反复执行该程序，每次执行时，输入不同之N值，比较所的图形的区别，由此可以观察到吉伯斯现象的特征。

4 用MATLAB实现CTFT及其逆变换的计算 4.1 用MATLAB实现CTFT的计算

MATLAB进行傅里叶变换有两种方法，一种利用符号运算的方法计算，另一种是数值计算，本实验要求采用数值计算的方法来进行傅里叶变换的计算。严格来说，用数值计算的方法计算连续时间信号的傅里叶变换需要有个限定条件，即信号是时限信号（Time limited signal），也就是当时间|t|大于某个给定时间时其值衰减为零或接近于零，这个条件与前面提到的为什么不能用无限多个谐波分量来合成周期信号的道理是一样的。计算机只能处理有限大小和有限数量的数。

符号运算的方法： syms t v w x

fourier(1/t) % returns i*pi*(Heaviside(-w)-Heaviside(w)) fourier(exp(-x^2),x,t) % returns pi^(1/2)*exp(-1/4*t^2) fourier(exp(-t)*sym('Heaviside(t)'),v) % returns 1/(1+i*v) fourier(diff(sym('F(x)')),x,w) % returns i*w*fourier(F(x),x,w) 其中：

function f=Heaviside(t) f=(t>0);

采用数值计算算法的理论依据是：

? X(j?)????x(t)e?j?tdt?lim?x(kT)e?jk?TT

T?0k????若信号为时限信号，当时间间隔T取得足够小时，上式可演变为：

X(j?)?Tk??N?x(kT)eN?jk?T

?[x(t1),x(t2),?,x(t2N?1)]?[e?j?t1,e?j?t2,?,e?j?t2N?1]T

上式用MATLAB表示为：

X=x*exp(j*t’*w)*T

其中X为信号x(t)的傅里叶变换，w为频率Ω，T为时间步长。 相应的MATLAB程序：

T = 0.01; dw = 0.1; %时间和频率变化的步长 t = -10:T:10;

w = -4*pi:dw:4*pi;

X(j?)可以按照下面的矩阵运算来进行：

x=u(t)-u(t-2);

X=x*exp(-j*t'*w)*T; %傅里叶变换 X1=abs(X); %计算幅度谱 phai=angle(X); %计算相位谱

为了使计算结果能够直观地表现出来，还需要用绘图函数将时间信号x(t)，信号的幅度谱|X(j?)|和相位谱? X(j?)分别以图形的方式表现出来，并对图形加以适当的标注。

plot(w,X) figure plot(w,X1) figure

plot(w,phai)

4.2 用MATLAB实现傅里叶逆变换

连续时间傅里叶逆变换可用式2.7进行计算。式2.7重写如下：

1x(t)?2?????X(j?)ej?td?

从定义式可看出，其计算方法与傅里叶变换是一样的，因此可以采用同样的矩阵运算的方法来计算，即

x(t)=X(j?)*exp(j?’*t)*d?/(2*pi) 具体的MATLAB函数如下：

t = -4*pi:0.01:4*pi; % 指定信号的时间范围，此范围应根据信号的持续时间确定。 dw = 0.01; w = -4*pi:dw:4*pi;

X = input('Type in the expression of X(jw)');%u(w)-u(w-1) x = X* exp(j*w'*t)*dw/(2*pi); plot(t,x)

然后用绘图函数就可以绘制出逆变换得到的时域信号波形图。

三、实验内容和要求

实验前，必须首先阅读本实验原理，读懂所给出的全部范例程序。实验开始时，先在计算机上运行这些范例程序，观察所得到的信号的波形图。并结合范例程序应该完成的工作，

进一步分析程序中各个语句的作用，从而真正理解这些程序。 实验前，一定要针对下面的实验项目做好相应的实验准备工作，包括事先编写好相应的实验程序等事项。

Q3-1 编写程序Q3_1，绘制下面的信号的波形图：

?111n?(0t)?cos3(?0t)?cos5(?0t)????sin()cosn( x(t)?cos??0t)

35n2n?1其中，?0 = 0.5π，要求将一个图形窗口分割成四个子图，分别绘制cos(?0t)、cos(3?0t)、cos(5?0t)

和x(t) 的波形图，给图形加title，网格线和x坐标标签，并且程序能够接受从键盘输入的和式中的项数。

抄写程序Q3_1如下：

t = -pi:0.001:pi; w0=0.5*pi; subplot(221) plot(t,cos(w0*t)) subplot(222)

plot(t,cos(3*w0*t)) subplot(223)

plot(t,cos(5*w0*t))

N= input('Type in the number of the harmonic components N = '); x=0;

for n=1:N

x=x+1/n*sin(n*pi/2)*cos(n*w0*t); end

subplot(224) plot(t,x)

执行程序Q3_1所得到的图形如下：

Q3-2 给程序Program3_1增加适当的语句，并以Q3_2存盘，使之能够计算例题3-1中的

周期方波信号的傅里叶级数的系数，并绘制出信号的幅度谱和相位谱的谱线图。

通过增加适当的语句修改Program3_1而成的程序Q3_2抄写如下：

执行程序Q3_2得到的图形

此处粘帖执行程序Q3_2所得到的图形

Q3-3 反复执行程序Program3_2，每次执行该程序时，输入不同的N值，并观察所合成的

周期方波信号。通过观察，你了解的吉伯斯现象的特点是：

Q3-4给定如下周期信号：

x(t)1t?2?0.20.22

仿照程序Program3_1，编写程序Q3_4，以计算x(t)的傅里叶级数的系数。

程序Q3_4如下：

clear, close all

T = 2; dt = 0.001; t = -4:dt:4; x1 = u(t+0.2) - u(t-0.2-dt); x = 0;

for m = -1:1 % Periodically extend x1(t) to form a periodic signal x = x + u(t+0.2-m*T) - u(t-0.2-m*T-dt); end

w0 = 2*pi/T;

N = 10; % The number of the harmonic components L = 2*N+1;

for k = -N: N; % Evaluate the Fourier series coefficients ak ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt; end

phi = angle(ak);

m=-N:N;

stem(m,abs(ak)) figure

stem(m,phi) figure plot(t,x)

axis([-4 4 -0.2 1.2])

执行程序Q3_4所得到的x(t)的傅里叶级数的ak从-10到10共21个系数如下：

Q3-5 仿照程序Program3_2，编写程序Q3_5，计算并绘制出原始信号x(t) 的波形图，用

有限项级数合成的y(t) 的波形图，以及x(t) 的幅度频谱和相位频谱的谱线图。

编写程序Q3_5如下：

clear,close all

T = 2; dt = 0.001; t = -2:dt:2; x1 = t.*(u(t)-u(t-1-dt)); x = 0; for m = -1:1

x = x + (t-m*T).*(u(t-m*T) - u(t-1-m*T-dt)); % Periodically extend x1(t) to form a periodic signal end

w0 = 2*pi/T;

N = input('Type in the number of the harmonic components N = '); L = 2*N+1; for k = -N:1:N;

ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt; end

phi = angle(ak); y=0;

for q = 1:L; % Synthesis the periodic signal y(t) from the finite Fourier series

y = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T);% y = y+ak(q)*exp(j*(q-1-N)*2*pi*t/T); end;

subplot(221),

plot(t,x), title('The original signal x(t)'), axis([-2,2,-0.2,1.2]), subplot(223),

plot(t,y), title('The synthesis signal y(t)'), axis([-2,2,-0.2,1.2]), xlabel('Time t'), subplot(222)

k=-N:N; stem(k,abs(ak),'k.'), title('The amplitude |ak| of x(t)'), axis([-N,N,-0.1,0.6]) subplot(224)

stem(k,phi,'r.'), title('The phase phi(k) of x(t)'), axis([-N,N,-2,2]), xlabel('Index k')

执行程序Q3_5，输入N = 10所得到的图形如下：

Q3-6给定时限信号:

x(t)?cos(?2t)[u(t?1)?u(t?1)]

编写MATLAB程序Q3_6，计算并绘制出信号x(t)的时域波形、幅度谱和相位谱。

程序Q3_6抄写如下

T = 0.01; dw = 0.1; %时间和频率变化的步长 t = -10:T:10;

w = -4*pi:dw:4*pi;

x=cos(pi*t/2).*(u(t+1)-u(t-1));

X=x*exp(-j*t'*w)*T; %傅里叶变换 X1=abs(X); %计算幅度谱 phai=angle(X); %计算相位谱 plot(t,x) figure plot(w,X1) figure

plot(w,phai)

执行程序Q3_6，得到的信号时域波形、幅度谱和相位谱如下：

实验四 连续系统的复频域分析

一、实验目的

1、掌握拉普拉斯变换的物理意义、基本性质及应用；

2、掌握用拉普拉斯变换求解连续时间LTI系统的时域响应；

3、掌握系统函数的概念，掌握系统函数的零、极点分布（零、极点图）与系统的稳定性、时域特性等之间的相互关系；

4、掌握用MATLAB对系统进行变换域分析的常用函数及编程方法。

基本要求：掌握拉普拉斯变换及其基本性质，掌握应用拉普拉斯变换求解系统的微分方程，能够自己编写程序完成对系统时域响应的求解。掌握并理解系统函数的概念，掌握系统函数零极点与系统时域和频域特性之间的关系，能够编写程序完成对系统的一些主要特性如稳定性、因果性等的分析。

二、实验原理及方法

1、连续时间LTI系统的复频域描述

拉普拉斯变换（The Laplace transform）主要用于系统分析。描述系统的另一种数学模型就是建立在拉普拉斯变换基础上的“系统函数（System Function）”——H(s)：

H(s)?Y(s)?系统响应的拉氏变换L?y(t)?X(s)?系统激励信号的拉氏变换L?x(t)? 4.1

系统函数H(s)的实质就是系统单位冲激响应（Impulse Response）h(t)的拉普拉斯变换。因此，系统函数也可以定义为：

? H(s)????h(t)e?stdt 4.2

所以，系统函数H(s)的一些特点是和系统的时域响应h(t)的特点相对应的。在教材中，我们求系统函数的方法，除了按照拉氏变换的定义式的方法之外，更常用的是根据描述系统的线性常系数微分方程（Linear Constant-Coefficient Deferential Equation），经过拉氏变换之后得到系统函数H(s)。

假设描述一个连续时间LTI系统的线性常系数微分方程为：

dky(t)Mdkx(t) 4.3 ak??bk?kkdtdtk?0k?0N对式4.3两边做拉普拉斯变换，则有

?asY(s)??bskkkk?0k?0NMkX(s)

即 H(s)?Y(s)?X(s)?bskMk?ak?0k?0N 4.4

ksk式4.4告诉我们，对于一个能够用线性常系数微分方程描述的连续时间LTI系统，它的系统函数是一个关于复变量s的有理多项式的分式，其分子和分母的多项式系数与系统微分方程左右两端的系数是对应的。根据这一特点，可以很容易的根据微分方程写出系统函数表达式，或者根据系统函数表达式写出系统的微分方程。

系统函数H(s)大多数情况下是复变函数，因此，H(s)可以有多种表示形式：

1、直角坐标形式： H(s)?Re(s)?jIm(s)

k?(s?zj)2、零极点形式： H(s)?j?1M?(s?p)ii?1N

3、部分分式和形式： H(s)?Ak（假设系统的N>M，且无重极点） ?s?sk?0kN根据我们所要分析的问题的不同，可以采用不同形式的系统函数H(s)表达式。

在MATLAB中，表达系统函数H(s)的方法是给出系统函数的分子多项式和分母多项式的系数向量。由于系统函数的分子和分母的多项式系数与系统微分方程左右两端的系数是

对应的，因此，用MATLAB表示系统函数，就是用系统函数的两个系数向量来表示。

应用拉普拉斯变换分析系统的主要内容有：

1、分析系统的稳定性； 2、分析系统的频率响应。

分析方法主要是通过绘制出系统函数的零极点分布图，根据零极点分布情况，判断系统的稳定性。

MATLAB中有相应的复频域分析函数，下面简要介绍如下：

[z,p,k] = tf2zp(num,den)：求系统函数的零极点，返回值z为零点行向量，p为极点

行向量，k为系统传递函数的零极点形式的增益。num为系统函数分子（numerator）多项式的系数向量，den为系统函数分母（denominator）多项式系数向量。

H = freqs(num,den,w)：计算由num,den描述的系统的频率响应特性曲线。返回值H

为频率向量规定的范围内的频率响应向量值。如果不带返回值H，则执行此函数后，将直接在屏幕上给出系统的对数频率响应曲线（包括幅频特性曲线和相频特性曲线）。

clear a=[1];

b=[1 0.4 0.08];

w=linspace(-2*pi,2*pi,512); h=freqs(a,b,w);

subplot(2,1,1),plot(w/pi,abs(h)) subplot(2,1,2),plot(w/pi,angle(h));

[x,y] = meshgrid(x1,y1)：用来产生绘制平面图的区域，由x1，y1来确定具体的区

域范围，由此产生s平面区域。

[X,Y] = meshgrid(-2:.2:2, -2:.2:2); Z = X .* exp(-X.^2 - Y.^2); surf(X,Y,Z)

mesh(x,y,fs)：绘制系统函数的零极点曲面图。 [X,Y] = meshgrid(-3:0.1:3,-3:0.1:3);

Z=sin(sqrt(X.^2+Y.^2)+eps)./(sqrt(X.^2+Y.^2)+eps); mesh(X,Y,Z)

mesh和surf命令都可以绘出某一区间内的完整曲面。它们的调用方法类似，不同

的是，mesh命令绘制的图形是一个一排排的彩色曲线组成的网格图，而surf命令绘制得到的是着色的三维曲面。

H = impulse(num,den)：求系统的单位冲激响应，不带返回值，则直接绘制响应曲

线，带返回值则将冲激响应值存于向量H之中。

2、系统函数的零极点分布图

系统函数的零极点图（Zero-pole diagram）能够直观地表示系统的零点和极点在s平面上的位置，从而比较容易分析系统函数的收敛域（Regin of convergence）和稳定性（stablity）。

下面给出一个用于绘制连续时间LTI系统的零极点图的扩展函数splane(num,den)：

% splane

% This function is used to draw the zero-pole plot in the s-plane function splane(num,den)

p = roots(den); % Determine the poles q = roots(num); % Determine the zeros p = p'; q = q'; %原来的p、q为一列矩阵 x = max(abs([p q])); % Determine the range of real-axis x = x+1;

y = x; % Determine the range of imaginary-axis plot([-x x],[0 0],':');hold on; % Draw the real-axis plot([0 0],[-y y],':');hold on; % Draw the imaginary-axis

plot(real(p),imag(p),'x');hold on; % Draw the poles plot(real(q),imag(q),'o');hold on; % Draw the zeros title('zero-pole plot');

xlabel('Real Part');ylabel('Imaginal Part') axis([-x x -y y]); % Determine the display-range

对于一个连续时间LTI系统，它的全部特性包括稳定性、因果性（Causality）和它具有何种滤波特性（Frequency-domain aspect）等完全由它的零极点在s平面上的位置所决定。

3、拉普拉斯变换与傅里叶变换之间的关系

根据课堂上所学的知识可知，拉普拉斯变换与傅里叶变换之间的关系可表述为：傅里叶变换是信号在虚轴上的拉普拉斯变换，也可用下面的数学表达式表示 H(j?)?H(s)s?j? 4.5

上式表明，给定一个信号h(t)，如果它的拉普拉斯变换存在的话，它的傅里叶变换不一定存在，只有当它的拉普拉斯变换的收敛域包括了整个虚轴，则表明其傅里叶变换是存在的。下面的程序可以以图形的方式，表现拉普拉斯变换与傅里叶变换的这种关系。

% Relation_ft_lt

% This program is used to observe the relationship between the Fourier transform % and the Laplace transform of a rectangular pulse. %f(t)=u(t)-u(t-2) clear, close all, a = -0:0.1:5; b = -20:0.1:20;

[a, b] = meshgrid (a, b); c = a+i*b; %确定绘图区域

1?e?2sc = (1-exp (-2* (c+eps)))./ (c+eps); %计算拉普拉斯变换F(s)?

sc = abs (c); subplot (211) mesh (a,b,c); %绘制曲面图 surf (a,b,c); view (-60,20) %调整观察视角 axis ([-0,5,-20,20,0,2]);

title ('The Laplace transform of the rectangular pulse'); w = -20:0.1:20;

Fw = (2*sin(w+eps).*exp(i*(w+eps)))./(w+eps);%计算傅里叶变换

F(j?)?2Sa(?)ej??2sin(?)?ej?

subplot (212); plot(w,abs(Fw))

title ('The Fourier transform of the rectangular pulse') xlabel ('frequence w')

上面的程序不要求完全读懂，重点是能够从所得到的图形中，观察拉和理解普拉斯变换与傅里叶变换之间的相互关系就行。

4、系统函数的极点分布与系统的稳定性和因果性之间的关系

一个稳定的LTI系统，它的单位冲激响应h(t)满足绝对可积条件，即

?

???h(t)dt?? 4.6

同时，我们还应该记得，一个信号的傅里叶变换的存在条件就是这个信号满足绝对可积条件，所以，如果系统是稳定的话，那么，该系统的频率响应也必然是存在的。又根据傅里叶变换与拉普拉斯变换之间的关系，可进一步推理出，稳定的系统，其系统函数的收敛域必然包括虚轴。稳定的因果系统，其系统函数的全部极点一定位于s平面的左半平面。

所以，对于一个给定的LTI系统，它的稳定性、因果性完全能够从它的零极点分布图上直观地看出。

例题4-1：已知一个因果的LTI系统的微分方程为

d6y(t)d5y(t)d4y(t)d3y(t)d2y(t)dy(t)?10?48?148?306?401?262y(t)?262x(t)dtdt6dt5dt4dt3dt2

编写程序，绘制出系统的零极点分布图，并说明它的稳定性如何。

解：这是一个高阶系统，显然手工计算它的极点是很困难的。可以利用前面给出的扩展函

数splane()，来绘制系统的零极点分布图。范例程序如下：

% Program4_1

% This program plots the zero-pole diagram of an LTI system described % by the linear constant-coefficient differential equation clear, close all, b = 262;

a = [1 10 48 148 306 401 262]; subplot (221) splane (b,a)

title ('The zero-pole diagram')

执行该程序后，得到系统的零极点分布图如图4.1所示。由于已知该系统是因果系统，从零极点分布图上看，它的全部极点都位于s平面的左半平面上，所以系统是稳定的。

图4.1

然后，直接在命令窗口键入

>> roots(a)

回车后，就得到系统的极点为：

ans =

-0.5707 + 2.4716i -0.5707 - 2.4716i -2.7378 + 0.0956i -2.7378 - 0.0956i -1.6915 + 1.6014i -1.6915 - 1.6014i

若题目中没有说明该系统是否是因果的，则需要做详细的分析。从零极点分布图上可以看出，该系统可能的收敛域共有四种可能，另外三种可能如下：

（a） 收敛域为Re{s}< -2.7378，此种情况说明，该系统是一个反因果系统（Anticausal system），由于收敛域不包含虚轴，故此系统是不稳定的。

（b）、（c）收敛域为 -2.7378 < Re{s} < -1.6915和-1.6915 < Re{s} < -0.5707，此两种情况说明该系统是一个单位冲激响应为双边信号的非因果系统，收敛域仍不包含虚轴，所以，系统是不稳定的。

总之，系统的稳定性主要取决于系统函数的收敛域是否包含整个虚轴，而系统的因果性则取决于系统极点位置的分布。

需要特别强调的是，MATLAB总是把由分子和分母多项式表示任何系统都当作是因果

系统。所以，利用impulse ()函数求得的单位冲激响应总是因果信号。

5、系统函数的零极点分布与系统的滤波特性

系统具有何种滤波特性，主要取决于系统的零极点所处的位置。没有零点的系统，通常是一个低通滤波器。

例题4-2 已知一个系统的系统函数为

H(s)?1 s?1图4.2

显然，这是一个一阶系统，无零点。为了确定该系统具有何种滤波特性，需要把系统的频率响应特性曲线绘制出来加以判断。借助实验三中的范例程序Program3_2，可以绘制系统的频率特性曲线如图4.2所示。

% Program4_2

% This Program is used to compute and draw the plots of the frequency response % of a continuous-time system

b = [1]; % The coefficient vector of the right side of the differential equation a = [1 1]; % The coefficient vector of the left side of the differential equation [H,w] = freqs(b,a); % Compute the frequency response H Hm = abs(H); % Compute the magnitude response Hm

phai = angle(H); % Compute the phase response phai

Hr = real(H); % Compute the real part of the frequency response

Hi = imag(H); % Compute the imaginary part of the frequency response subplot(221)

plot(w,Hm), grid on, title('Magnitude response'), xlabel('Frequency in rad/sec') subplot(223)

plot(w,phai), grid on, title('Phase response'), xlabel('Frequency in rad/sec') subplot(222)

plot(w,Hr), grid on, title('Real part of frequency response'), xlabel('Frequency in rad/sec') subplot(224)

plot(w,Hi), grid on, title('Imaginary part of frequency response'), xlabel('Frequency in rad/sec')

图4.3

通过编程，可以将系统的零极点分布图和系统的频率响应特性以及系统的单位冲激响应特性绘制在一个图形窗口的各个子图中，这样便于观察系统的零极点分布情况与系统的时域和频域之间的关系。如图4.4所示。

clear, close all, b = 262;

a = [1 10 48 148 306 401 262]; subplot (221) splane (b,a)

title ('The zero-pole diagram') [H,w] = freqs(b,a); Hm = abs(H); phai = angle(H); subplot(222) impulse(b,a) subplot(223)

plot(w,Hm), title('Frequency response'), xlabel('Frequency in rad/sec')

subplot(224)

plot(w,phai), title('Phase response'), xlabel('Frequency in rad/sec')

图4.4 系统的零极点分布与系统的单位冲激响应、频率响应

6、拉普拉斯逆变换的计算

我们已经知道，直接用拉普拉斯逆变换（Inverse transform）的定义公式计算逆变换是很困难的，通常的计算拉普拉斯逆变换的方法是长除法（Long division）和部分分式分解法（Partial fraction decomposition）。MATLAB的内部函数residue()可以帮助我们完成拉普拉斯逆变换的计算。

例题4-3 已知某信号的拉普拉斯变换表达式为

X(s)?求该信号的时域表达式。

1 2s?3s?2解：由于题目没有指定收敛域，所以必须考虑所有可能的情况。为此，可以先计算出该信

号的拉普拉斯变换表达式的极点。很显然，X(s)有两个极点，分别为 s = -1，s = -2。零极点分布图如例题图4.5所示。

在MATLAB命令窗口键入：

>> b = 1;

>> a = [1 3 2];

>> [r, p, k] = residue (b, a)

命令窗口立即给出计算结果为：

r = -1 1

p = -2 -1 k = 例题图4.5 [ ]

根据r（系数）、p（极点）、k（余数）之值，可以写出X(s)的部分分式和的表达式为： X(s)??11? s?2s?1然后根据不同的ROC，可写出X(s)的时域表达式x(t)。

第一种情况，ROC为 Re{s} < -2，则x(t)为反因果信号，其数学表达式为 x(t)?e?2tu(?t)?e?tu(?t)

第二种情况，ROC为 -2 < Re{s} < -1，则x(t)为双边非因果信号，其数学表达式为 x(t)??e?2tu(t)?e?tu(?t)

第三种情况，ROC为 Re{s} > -1，则x(t)为因果信号，其数学表达式为 x(t)??e?2tu(t)?e?tu(t)

在这个例题中，函数residue()仅仅完成了部分分式分解的任务，至于逆变换的数学表达式的结果是什么，还得结合收敛域的不同才能写出。

如果X(s)的分子的阶不小于分母的阶，则k将不等于一个空矩阵，例如，当

s3X(s)?2时，我们在命令窗口中键入：

s?3s?2>> b = [1 0 0 0]; >> a = [1 3 2];

>> [r,p,k]=residue(b,a)

则：

r = 8 -1 p = -2 -1 k =

1 -3

这里的k = [1 3]，实际上是将X(s)做了一个长除法后，得到的商的多项式。所以，根据上面的r、p、k之值，可写出X(s)的部分分式和的表达式为：

X(s)?s?3?81? s?2s?1有关函数residue()的详细用法，可通过在线帮助加以了解。

三、实验内容及步骤

实验前，必须首先阅读本实验原理，了解所给的MATLAB相关函数，读懂所给出的全部范例程序。实验开始时，先在计算机上运行这些范例程序，观察所得到的信号的波形图。并结合范例程序所完成的工作，进一步分析程序中各个语句的作用，从而真正理解这些程序。

实验前，一定要针对下面的实验项目做好相应的实验准备工作，包括事先编写好相应的实验程序等事项。

Q4-1 将绘制零极点图的扩展函数文件splane以splane为文件名存盘。

function splane(num,den)

p = roots(den); % Determine the poles q = roots(num); % Determine the zeros p = p'; q = q'; x = max(abs([p q])); % Determine the range of real-axis x = x+1;

y = x; % Determine the range of imaginary-axis plot([-x x],[0 0],':');hold on; % Draw the real-axis plot([0 0],[-y y],':');hold on; % Draw the imaginary-axis plot(real(p),imag(p),'x');hold on; % Draw the poles plot(real(q),imag(q),'o');hold on; % Draw the zeros title('zero-pole plot');

xlabel('Real Part');ylabel('Imaginal Part') axis([-x x -y y]); % Determine the display-range

Q4-2 运行程序

Relation_ft_lt，观察拉普拉斯变换与傅里叶变换之间的关系。在点击工具

条上的旋转按钮，再将鼠标放在曲面图上拖动图形旋转，从各个角度观察拉普拉斯曲面图形，

并同傅立叶变换的曲线图比较，加深对拉普拉斯变换与傅里叶变换之间关系的理解与记忆。

%f(t)=u(t)-u(t-2) clear, close all, a = -0:0.1:5; b = -20:0.1:20;

[a, b] = meshgrid (a, b); c = a+i*b; %确定绘图区域

c = (1-exp (-2* (c+eps)))./ (c+eps); %计算拉普拉斯变换 c = abs (c); subplot (211) mesh (a,b,c); %绘制曲面图 surf (a,b,c); view (-60,20) %调整观察视角 axis ([-0,5,-20,20,0,2]);

title ('The Laplace transform of the rectangular pulse'); w = -20:0.1:20;

Fw = (2*sin(w+eps).*exp(i*(w+eps)))./(w+eps); subplot (212); plot(w,abs(Fw))

title ('The Fourier transform of the rectangular pulse')

xlabel ('frequence w')

Q4-3 编写程序Q4_3，能够接受从键盘输入的系统函数的分子分母多项式系数向量，并绘

制出系统的零极点图、系统的单位冲激响应、系统的幅度频率响应和相位频率相应的图形。各个子图要求按照图4.4布置。

程序Q4_3抄写如下：

clear, close all,

num=input('系统函数分子多项式系数向量=') ; den=input('系统函数分母多项式系数向量='); subplot (221) splane (num,den)

title ('The zero-pole diagram') [H,w] = freqs(num,den); Hm = abs(H); phai = angle(H); subplot(222)

impulse(num,den) subplot(223)

plot(w,Hm), title('Frequency response'), xlabel('Frequency in rad/sec') subplot(224)

plot(w,phai), title('Phase response'), xlabel('Frequency in rad/sec')

Q4-4 执行程序编写Q4_3，输入因果的系统函数H(s)?量，绘制所得到的图形如下：

执行Q4_3所得到的图形

s的分子分母系数向

(s?1)(s?2)从上面的图形中可以看出，该系统的零点和极点分别位于：

从时域和零极点分布特征两个方面说明该系统是否是稳定的系统？

答：

从频率响应特性上看，该系统具有何种滤波特性？ 答：

Q4-5 执行程序编写Q4_3，输入因果的系统函数

12s?1a H(s)?32s?2s?2s?1此处a取1，执行程序Q4_3，输入该系统的分子分母系数向量，得到的图形如下：

从上面的图形中可以看出，该系统的零点和极点分别位于：

从时域和零极点分布特征两个方面说明该系统是否是稳定的系统？ 答：

从频率响应特性上看，该系统具有何种滤波特性？ 答：

改变系统函数中的a值，分别取0.6、0.8、4、16等不同的值，反复执行程序Q4-3，观察

系统的幅度频率响应特性曲线（带宽、过渡带宽和阻带衰减等），贴一张a = 4时的图形如下：

观察a取不同的值时系统的幅度频率响应特性曲线的变化（带宽、过渡带宽和阻带衰减等），

请用一段文字说明零点位置对系统滤波特性的这些影响。

答：

本文来源：https://www.bwwdw.com/article/ekoo.html