MATLAB与信号与系统实验指导书（李敏 教改版）

《MATLAB&信号与系统》实验指导书

《MATLAB&信号与系统》

实验指导书

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

前 言

长期以来，《信号与线性系统》课程一直采用单一理论教学方式，同学们依靠做习题来巩固和理解教学内容，虽然手工演算训练了计算能力和思维方法，但是由于本课程数学公式推导较多，概念抽象，常需画各种波形，做题时难免花费很多时间，现在，我们给同学们介绍一种国际上公认的优秀科技应用软件MATLAB，借助它我们可以在电脑上轻松地完成许多习题的演算和波形的绘制。

MATLAB的功能非常强大，我们此处仅用到它的一部分，在后续课程中我们还会用到它，在未来地科学研究和工程设计中有可能继续用它，所以有兴趣的同学，可以对MATLAB再多了解一些。

MATLAB究竟有那些特点呢？

1．高效的数值计算和符号计算功能，使我们从繁杂的数学运算分析中解脱出来； 2．完备的图形处理功能，实现计算结果和编程的可视化；

3．友好的用户界面及接近数学表达式的自然化语言，易于学习和掌握； 4．功能丰富的应用工具箱，为我们提供了大量方便实用的处理工具；

MATLAB的这些特点，深受大家欢迎，由于个人电脑地普及，目前许多学校已将它做为本科生必须掌握的一种软件。正是基于这些背景，我们编写了这本《信号与系统分析的MATLAB实现》实验指导书，内容包括信号的MATLAB表示、基本运算、系统的时域分析、频域分析、S域分析、Z域分析等。通过这些练习，同学们在学习《信号与线性系统》的同时，掌握MATLAB的基本应用，学会应用MATLAB的数值计算和符号计算功能，摆脱烦琐的数学运算，从而更注重于信号与系统的基本分析方法和应用的理解与思考，将课程的重点、难点及部分习题用MATLAB进行形象、直观的可视化计算机模拟与仿真实现，加深对信号与线性系统的基本原理、方法及应用的理解，为学习后续课程打好基础。另外同学们在进行实验时，最好事先预习一些MATLAB的有关知识，以便更好地完成实验，同时实验中也可利用MATLAB的help命令了解具体语句以及指令的使用方法。

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

实验一 常用基本信号的MATLAB表示和运算

一、实验目的

1． 学会用MATLAB表示常用连续信号的方法； 2． 学会用MATLAB进行信号基本运算的方法； 二、实验原理

1． 信号的MATLAB表示 (1).连续信号的MATLAB表示

MATLAB提供了大量的生成基本信号的函数，例如指数信号、正余弦信号。

表示连续时间信号有两种方法，一是数值法，二是符号法。数值法是定义某一时间范围和取样时间间隔，然后调用该函数计算这些点的函数值，得到两组数值矢量，可用绘图语句画出其波形；符号法是利用MATLAB的符号运算功能，需定义符号变量和符号函数，运算结果是符号表达的解析式，也可用绘图语句画出其波形图。

例1-1指数信号 指数信号在MATLAB中用exp函数表示。

如f(t)?Aeat，调用格式为 ft=A*exp(a*t) 程序是

A=1; a=-0.4;

t=0:0.01:10; %定义时间点

ft=A*exp(a*t); %计算这些点的函数值

plot(t,ft); %画图命令，用直线段连接函数值表示曲线 grid on; %在图上画方格

例1-2 正弦信号 正弦信号在MATLAB中用 sin 函数表示。

调用格式为 ft=A*sin(w*t+phi) A=1; w=2*pi; phi=pi/6;

t=0:0.01:8; %定义时间点

ft=A*sin(w*t+phi); %计算这些点的函数值 plot(t,ft); %画图命令 grid on; %在图上画方格

例1-3 抽样信号 抽样信号Sa(t)=sin(t)/t在MATLAB中用 sinc 函数表示。

定义为 Sa(t)?sinc(t/?)

t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on;

axis([-10,10,-0.5,1.2]); %定义画图范围，横轴，纵轴 title('抽样信号') %定义图的标题名字

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

例1-4 三角信号 三角信号在MATLAB中用 tripuls 函数表示。

调用格式为 ft=tripuls(t,width,skew)，产生幅度为1，宽度为width，且以0为中心左右各展开width/2大小，斜度为skew的三角波。width的默认值是1，skew的取值范围是-1~+1之间。一般最大幅度1出现在t=(width/2)*skew的横坐标位置。 t=-3:0.01:3;

ft=tripuls(t,4,0.5); plot(t,ft); grid on; axis([-3,3,-0.5,1.5]);

例1-5虚指数信号 调用格式是f=exp((j*w)*t) t=0:0.01:15;

w=pi/4;

X=exp(j*w*t);

Xr=real(X); %取实部 Xi=imag(X); %取虚部 Xa=abs(X); %取模 Xn=angle(X); %取相位

subplot(2,2,1),plot(t,Xr),axis([0,15,-(max(Xa)+0.5),max(Xa)+0.5]), title('实部');

subplot(2,2,3),plot(t,Xi),axis([0,15,-(max(Xa)+0.5),max(Xa)+0.5]), title('虚部');

subplot(2,2,2), plot(t,Xa),axis([0,15,0,max(Xa)+1]),title('模');

subplot(2,2,4),plot(t,Xn),axis([0,15,-(max(Xn)+1),max(Xn)+1]),title('相角'); %subplot(m,n,i) 命令是建立m行n列画图窗口，并指定画图位置i 例1-6复指数信号 调用格式是f=exp((a+j*b)*t) t=0:0.01:3;

a=-1;b=10;

f=exp((a+j*b)*t);

subplot(2,2,1),plot(t,real(f)),title('实部') subplot(2,2,3),plot(t,imag(f)),title('虚部') subplot(2,2,2),plot(t,abs(f)),title('模') subplot(2,2,4),plot(t,angle(f)),title('相角')

例1-7 矩形脉冲信号 矩形脉冲信号可用rectpuls函数产生，

调用格式为y=rectpuls(t,width),幅度是1，宽度是width，以t=0为对称中心。 t=-2:0.01:2; width=1;

ft=2*rectpuls(t,width); plot(t,ft) grid on;

例1-8 单位阶跃信号 单位阶跃信号u(t)用“t>=0”产生，调用格式为ft=(t>=0) t=-1:0.01:5; ft=(t>=0);

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

plot(t,ft); grid on; axis([-1,5,-0.5,1.5]); 例1-9 正弦信号符号算法

syms t %定义符号变量t y=sin(pi/4*t) %符号函数表达式 ezplot(y,[-16,16]) %符号函数画图命令 或者

f=sym('sin(pi/4*t)') %定义符号函数表达式 ezplot(f,[-16,16])

例1-10单位阶跃信号 MATTLAB符号数学函数Heaviside表示阶跃信号，但要画图需在工作目录创建Heaviside的M文件

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

保存，文件名是Heaviside ，调用该函数即可画图，例 t=-1:0.01:3; f=heaviside(t); plot(t,f)

axis([-1,3,-0.2,1.2]) 或者

y=sym('Heaviside(t)'); ezplot(y,[-1,5]);grid on (2)离散信号的MATLAB表示 例1-11 单位脉冲序列

单位脉冲序列的表达式: 延迟ks的单位脉冲序列表达式:

?(k)???1,?0,k?0?1, x1(k)??(k?ks)??k?其余?0,k?ksk?其余

％单位脉冲序列m文件

clear,k0=0;kf=10;ks=3; %本例取ks=3, k1=k0:kf;

x1=［zeros(1,ks-k0),1,zeros(1,kf-ks)］;%单位脉冲序列的产生 stem(k1,x1,′.′);title(′单位脉冲序列′)%绘图 例1-12 单位阶跃序列

单位阶跃序列的表达式: 延迟ks的单位阶跃序列表达式: ?1k?0?1,k?ks??k???x2 ?k????k－ks????0k?0?0,k?ks

％本例取ks=3。单位阶跃序列m文件 clear,k0=0;kf=10;ks=3;

k2=k0:kf;x2=[zeros(1,ks-k0),ones(1,kf-ks+1)]; %单位阶跃序列的产生 stem(k2,x2,'.');title('单位阶跃序列') %绘图

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

例1-13 复指数序列

复指数序列的表达式: ?????k?ek?0 ??x3k??k?0 ?0若ω=0,它是实指数序列;若α=0,则为虚指数序列,其实部为余弦序列,虚部为正弦序列。本例取α=-0.2,ω=0.5,该复指数序列的实部和虚部如图所绘，其m文件如下 ％复指数序列m文件 clear,k0=0;kf=20;ks=3;

k3=k0:kf;x3=exp((-0.2+0.5j)*k3);%复指数序列的产生

subplot(1,2,1),stem(k3,real(x3),′.′);line(［0,10］,［0,0］)%绘图 xlabel(′实部′)

subplot(1,2,2),stem(k3,imag(x3),′.′);line(［0,10］,［0,0］)%绘图 xlabel(′虚部′)

2． 信号基本运算的MATLAB实现

信号基本运算是乘法、加法、尺度、反转、平移、微分、积分，实现方法有数值法和符号法.

例1-14 以f(t)为三角信号为例，求f(2t) , f(2-2t)

t=-3:0.001:3;

ft=tripuls(t,4,0.5);

subplot(3,1,1); plot(t,ft); grid on; title ('f(t)');

ft1= tripuls(2*t,4,0.5);

subplot(3,1,2); plot(t,ft1); grid on; title ('f(2t)');

ft2= tripuls(2-2*t,4,0.5);

subplot(3,1,3); plot(t,ft2); grid on; title ('f(2-2t)');

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

f(t)10.50-310.50-310.50-3-2-10f(2t)123-2-10f(2-2t)123-2-10123

例1-15 已知f1(t)=sinwt , f2(t)=sin8wt , w=2pi , 求f1(t)+f2(t)和f1(t)f2(t) 的波形图 w=2*pi;

t=0:0.01:3; f1=sin(w*t); f2=sin(8*w*t); subplot(211)

plot(t,f1+1,':',t,f1-1,':',t,f1+f2) grid on,title('f1(t)+f2(t))') subplot(212)

plot(t,f1,':',t,-f1,':',t,f1.*f2) grid on,title('f1(t)*f2(t)')

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

f1(t)+f2(t))210-1-200.511.5f1(t)*f2(t)22.5310.50-0.5-100.511.522.53

离散序列的差分与求和 离散序列的差分 ? f ? k ? ? f ? 1 ? ,在MATLAB中用diff函数实现,其调用格式???kf ?k为:

y=diff(f)

离散序列的求和

?f(k) 与信号相加运算不同,求和运算是把k1和k2之间的所有样本

k1k2f［k］加起来,在MATLAB中用sum函数实现,其调用格式为:

y=sum(f(k1:k2))

例1-16 用MATLAB计算指数信号(-1.6)kε［k］的能量。 解:离散信号的能量定义为

N 2E?limf?k?N?? k??N其MATLAB实现如下: ％program10.2－3 k=0:10;A=1;a=-1.6; fk=A*a.^k;

W=sum(abs(fk).^2) 运行结果为

W=1.9838e+004

?三、实验内容

1． 验证实验原理中程序

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

2． 画出信号波形

（1）f(t)?(2?e?2t)u(t) （2）f(t)?(1?cos?t)[u(t)?u(t?2)] 3．信号f(t)?(2?e?2t)u(t)，画出f(2t)、f(2?t)波形 四、 实验报告要求

1、 简述实验目的及实验原理。

2、 计算相应的信号运算的理论值，并与实验结果进行比较。 3、 记录实验结果，分析系统作用。 4、 写出程序清单。 5、 收获与建议。

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

实验二 连续时间信号与系统的时域分析

一、实验目的

1．学会用MATLAB求解连续系统的零状态响应； 2. 学会用MATLAB求解冲激响应及阶跃响应； 3．学会用MATLAB实现连续信号卷积的方法； 二、实验原理

1．连续时间系统零状态响应的数值计算

我们知道，LTI连续系统可用如下所示的线性常系数微分方程来描述，

?ayii?0N(i)(t)??bjf(j)(t)

j?0M在MATLAB中，控制系统工具箱提供了一个用于求解零初始条件微分方程数值解的函数lsim。其调用格式

y=lsim(sys,f,t)

式中，t表示计算系统响应的抽样点向量，f是系统输入信号向量，sys是LTI系统模型，用来表示微分方程，差分方程或状态方程。其调用格式

sys=tf(b,a)

式中，b和a分别是微分方程的右端和左端系数向量。例如，对于以下方程:

a3y'''(t)?a2y''(t)?a1y'(t)?a0y(t)?b3f'''(t)?b2f''(t)?b1f'(t)?b0f(t)

可用a?[a3,a2,a1,a0];b?[b3,b2,b1,b0]; sys?tf(b,a) 获得其LTI模型。

注意，如果微分方程的左端或右端表达式中有缺项，则其向量a或b中的对应元素应为零，不能省略不写，否则出错。

例2-1 已知某LTI系统的微分方程为 y’’(t)+ 2y’(t)+100y(t)=f(t)

其中，y(0)?y(0)?0,f(t)?10sin(2?t)，求系统的输出y(t).

解：显然，这是一个求系统零状态响应的问题。其MATLAB计算程序如下： ts=0;te=5;dt=0.01; sys=tf([1],[1,2,100]); t=ts:dt:te;

f=10*sin(2*pi*t); y=lsim(sys,f,t); plot(t,y);

xlabel('Time(sec)'); ylabel('y(t)');

2．连续时间系统冲激响应和阶跃响应的求解 在MATLAB中，对于连续LTI系统的冲激响应和阶跃响应，可分别用控制系统工具箱提供的函数impluse和step来求解。其调用格式为

乐山师范学院.物理与电子工程学院

'

《MATLAB&信号与系统》实验指导书

y=impluse(sys,t) y=step(sys,t)

式中，t表示计算系统响应的抽样点向量，sys是LTI系统模型。

例2-2已知某LTI系统的微分方程为 y’’(t)+ 2y’(t)+100y(t)=10f(t)

求系统的冲激响应和阶跃响应的波形. 解：ts=0;te=5;dt=0.01; sys=tf([10],[1,2,100]);

t=ts:dt:te;

h=impulse(sys,t); figure; plot(t,h);

xlabel('Time(sec)'); ylabel('h(t)');

g=step(sys,t); figure; plot(t,g);

xlabel('Time(sec)');

ylabel('g(t)');

3. 用MATLAB实现连续时间信号的卷积

信号的卷积运算有符号算法和数值算法，此处采用数值计算法，需调用MATLAB 的conv( )函数近似计算信号的卷积积分。连续信号的卷积积分定义是 f(t)?f1(t)?f2(t)?????f1(?)f2(t??)d?

如果对连续信号f1(t)和f2(t)进行等时间间隔?均匀抽样，则f1(t)和f2(t)分别变为离散时间信号f1(m?)和f2(m?)。其中，m为整数。当?足够小时，f1(m?)和f2(m?)既为连续时间信号f1(t)和f2(t)。因此连续时间信号卷积积分可表示为

f(t)?f1(t)?f2(t)??f1(?)f2(t??)d?????lim??0m????f(m?)?f1?

2(t?m?)??采用数值计算时，只求当t?n?时卷积积分f(t)的值f(n?)，其中，n为整数，既

f(n?)?

m?????f(m?)?f1m????2(n??m?)??

???f1(m?)?f2[(n?m)?]其中，

m????f(m?)?f1?2[(n?m)?]实际就是离散序列f1(m?)和f2(m?)的卷积和。当?足

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

够小时，序列f(n?)就是连续信号f(t)的数值近似，既 f(t)?f(n?)??[f1(n)?f2(n)]

上式表明，连续信号f1(t)和f2(t)的卷积，可用各自抽样后的离散时间序列的卷积再乘以抽样间隔?。抽样间隔?越小，误差越小。

例2-3用数值计算法求f1(t)?u(t)?u(t?2)与f2(t)?e?3tu(t)的卷积积分。

解：因为f2(t)?e?3tu(t)是一个持续时间无限长的信号，而计算机数值计算不可能计算真正的无限长信号，所以在进行f2(t)的抽样离散化时，所取的时间范围让f2(t)衰减到足够小就可以了，本例取t?2.5。程序是 dt=0.01; t=-1:dt:2.5;

f1=Heaviside(t)-Heaviside(t-2); f2=exp(-3*t).*Heaviside(t);

f=conv(f1,f2)*dt; n=length(f); tt=(0:n-1)*dt-2; subplot(221), plot(t,f1), grid on;

axis([-1,2.5,-0.2,1.2]); title('f1(t)'); xlabel('t') subplot(222), plot(t,f2), grid on;

axis([-1,2.5,-0.2,1.2]); title('f2(t)'); xlabel('t') subplot(212), plot(tt,f), grid on; title('f(t)=f1(t)*f2(t)'); xlabel('t')

由于f1(t)和f2(t)的时间范围都是从t=-1开始，所以卷积结果的时间范围从 t=-2开始，增量还是取样间隔?，这就是语句tt=(0:n-1)*dt-2的由来。 4.利用MATLAB中的Simulink进行系统时域特性仿真

已知或求解出系统传输函数H(s)，分子和分母都写出多项式的形式，例如：

H(s)?s?1 2s?1.3s?0.8在continuous的子库中选择传输函数（Transfer Fcn）子库，用鼠标把传输函数模块拖到 Untitled视窗中，置于激励信号源和示波器之间，然后用鼠标拖出的连线将信号源、传输函数模块和示波器等按照系统的要求连接起来即可。

如果需要对传输函数的参数进行设置，则可通过双击传输函数模块打开它的设置环境窗，可以分别设置分子、分母多项式的系数和阶数等。

在untitled窗口的菜单中选simulation的start功能则可执行仿真，待仿真完毕后，则可通过激活示波器观察到激励信号和响应信号的波形。

由于source库中没有单位冲激信号模块，所以由阶跃信号模块经过微分来产生。仿真冲激响应的系统中，改变输入信号为任意波形，例如周期矩形波等，将得到在任意激励下的零状态响应。 三、实验内容

1． 验证实验原理中所述的相关程序

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

2． 已知描述系统的微分方程和激励信号f(t)如下，试用解析法求系统的零状态响应y(t)，

并用MATLAB绘出系统零状态响应的时域仿真波形，验证结果是否相同

?(t)u(t) y''(t)?4y'(t)?4y(t)?f'(t)?3f(t); f(t)?exp3.已知两连续时间信号如下图所示，试用MATLAB求f(t)?f1(t)*f2(t)，并绘出f(t)的时域波形图。（设定取样时间间隔为p）

4、对如下连续时间系统

H(s)?s?1

s2?1.3s?0.8 通过Simulink仿真分别观察其单位冲激响应波形和周期矩形信号作用下的零状态响应波形。

四、 实验报告要求

1、 简述实验目的及实验原理。

2、 计算相应的冲激响应、零状态响应及卷积积分的理论值，并与实验结果进行比较。 3、 记录连续系统的冲激响应、零状态响应特性曲线，分析系统作用。 4、 写出程序清单。 5、 收获与建议。

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

实验三 信号的幅度调制及MATLAB实现

一. 实验目的

熟悉使用MATLAB软件来分析信号的调制问题并可视化相关结果，同时分析和对比有关结果。

二. 实验原理

设信号f(t)的频谱为F(j?)，现将f(t)乘以载波信号cos(?0t)，得到高频的已调信号y(t)，即：

y(t)?f(t)cos(?0t)

其中，f(t)称为调制信号。实现信号调制的原理图如下所示：

从频域上看，已调制信号y(t)的频谱为原调制信号f(t)的频谱搬移到??0处，幅度将为原F(j?)的12，即：

1f(t)cos(?0t)?{F[j(???0)]?F[j(???0)]}

2上式即为调制原理，也是傅里叶变换性质中“频移特性”的一种特别情形。这里采用的调制方法为抑制载波方式，即y(t)的频谱中不含有cos(?0t)的频率分量。

MATLAB提供了专门的函数modulate()用于实现信号的调制。调用格式为：

y=modulate(x,Fc,Fs,’method’) [y,t]=modulate(x,Fc,Fs)

其中，x为被调信号，Fc为载波频率，Fs为信号x的采样频率，method为所采用的调制方式，若采用幅度调制、双边带调制、抑制载波调制，则’method’为’am’或’amdsd-sc’。其执行算法为：

y=x*cos(2*pi*Fc*t)

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

其中，y为已调信号，t为函数计算时间间隔向量。

在MATLAB的实现程序中，为了观察f(t)及y(t)的频谱，可使用“信号处理工具箱函数”中估计信号的功率谱密度函数psd()，其格式是：

[Px,f]=psd(x,Nfft,Fs,window,noverlap,dflag)

其中，x是被调信号（即本例中的f(t)），Nfft指定快速傅氏变换FFT的长度，Fs为对信号的采样频率。后面三个参数的意义将涉及信号处理的更深的知识，在此暂不介绍。

三. 实验内容

1、 设信号f(t)?sin(100?t)，载波y(t)为频率为400Hz的余弦信号。试用MATLAB

实现调幅信号y(t)，并观察y(t)的频谱和f(t)的频谱，以及两者在频域上的关系。

2、 设f(t)?u(t?1)?u(t?1),f1(t)?f(t)cos(10?t),，试用MATLAB画出f(t),f1(t)的时域波形及其频谱，并观察傅里叶变换的频移特性。

四. 实验报告要求

1、 简述实验目的及实验原理。

2、 写出程序清单,列出实验结果图片，并进行必要的文字解释或

原理阐述。

3、 写出收获与建议。

附录：

1、内容1的MATLAB参考程序如下：

Fs=1000; Fc=400; N=1000; n=0:N-2; t=n/Fs;

x=sin(2*pi*50*t); subplot(221) plot(t,x); xlabel('t(s)'); ylabel('x');

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

title('被调信号'); axis([0 0.1 -1 1]) Nfft=1024;

window=hamming(512); noverlap=256; dflag='none';

[Pxx,f]=psd(x,Nfft,Fs,window,noverlap,dflag); subplot(222) plot(f,Pxx) xlabel('频率(Hz)');

ylabel('功率谱(X)');

title('被调信号的功率谱') grid

y=modulate(x,Fc,Fs,'am'); subplot(223) plot(t,y) xlabel('t(s)'); ylabel('y'); axis([0 0.1 -1 1]) title('已调信号')

[Pxx,f]=psd(y,1024,Fs,window,noverlap,dflag); subplot(224) plot(f,Pxx) xlabel('频率(Hz)');

ylabel('功率谱(Y)');

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

title('已调信号的功率谱'); grid

2、内容2的MATLAB参考程序如下：

R=0.005;t=-1.2:R:1.2;

f=Heaviside(t+1)-Heaviside(t-1); f1=f.*cos(10*pi*t); subplot(221) plot(t,f) xlabel('t'); ylabel('f(t)'); subplot(222); plot(t,f1); xlabel('t');

ylabel('f1(t)=f(t)*cos(10*pi*t)'); W1=40; N=1000; k=-N:N; W=k*W1/N; F=f*exp(-j*t'*W)*R; F=real(F);

F1=f1*exp(-j*t'*W)*R; F1=real(F1); subplot(223); plot(W,F); xlabel('w'); ylabel('F(jw)'); subplot(224); plot(W,F1); xlabel('w');

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

ylabel('F1(jw)');

实验四 连续时间信号与系统的频域分析

一、 实验目的

1、学会用MATLAB实现连续时间信号傅里叶变换 2、学会用MATLAB分析LTI系统的频域特性 3、学会用MATLAB分析LTI系统的输出响应 二、实验原理

1．傅里叶变换的MATLAB求解

MATLAB的symbolic Math Toolbox 提供了直接求解傅里叶变换及逆变换的函数fourier()及ifourier()两者的调用格式如下。 Fourier 变换的调用格式

F=fourier(f)：它是符号函数f的fourier变换默认返回是关于w的函数。

F=fourier(f，v):它返回函数F是关于符号对象v的函数，而不是默认的w，即

F(v)??????f(x)e?jvxdx

Fourier逆变换的调用格式

f=ifourier(F):它是符号函数F的fourier逆变换，默认的独立变量为w，默认返回是

关于x的函数。

f=ifourier(f,u):它的返回函数f是u的函数，而不是默认的x.

注意：在调用函数fourier()及ifourier()之前，要用syms命令对所用到的变量（如t,u,v,w）进行说明,即将这些变量说明成符号变量。 例4-1 求f(t)?e?2t的傅立叶变换 解: 可用MATLAB解决上述问题： syms t

Fw=fourier(exp(-2*abs(t))) 例4-2 求F(jw)?1的逆变换f(t) 1??2解： 可用MATLAB解决上述问题 syms t w

ft=ifourier(1/(1+w^2),t)

2．连续时间信号的频谱图

例4-3 求调制信号f(t)?AG?(t)cos?0t的频谱，式中

A?4,?0?12?,??1??,G?(t)?u(t?)?u(t?) 222解：MATLAB程序如下所示

ft=sym('4*cos(2*pi*6*t)*(Heaviside(t+1/4)-Heaviside(t-1/4))');

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

Fw=simplify(fourier(ft)) subplot(121)

ezplot(ft,[-0.5 0.5]),grid on subplot(122)

ezplot(abs(Fw),[-24*pi 24*pi]),grid

4 cos(2 ? 6 t) (Heaviside(t+1/4)-Heaviside(t-1/4))43210-1-20.2-3-4-0.50t0.50.10-500w5010.90.80.70.60.50.40.38 abs(w sin(1/4 w)/(w2-144 ?2))

用MATLAB符号算法求傅里叶变换有一定局限，当信号不能用解析式表达时，会提示出错，这时用MATLAB的数值计算也可以求连续信号的傅里叶变换，计算原理是

F(j?)??f(t)e????j?tdt?lim?f(n?)e?j?n??

??0n????当?足够小时，近似计算可满足要求。若信号是时限的，或当时间大于某个给定值时，信号已衰减的很厉害，可以近似地看成时限信号时，n的取值就是有限的，设为N，有

F(k)??f(n?)e?j?kn??n?0N?1,0?k?N,?k?2?k 是频率取样点 N?时间信号取样间隔?应小于奈奎斯特取样时间间隔，若不是带限信号可根据计算精度要求确定一个频率 W0为信号的带宽。

例4-4 用数值计算法求信号f(t)?u(t?1)?u(t?1)的傅里叶变换

解，信号频谱是F(j?)?2Sa(?)，第一个过零点是?，一般将此频率视为信号的带宽，若将精度提高到该值的50倍，既W0=50?，据此确定取样间隔，??1?0.02 2F0乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

R=0.02;t=-2:R:2;

f=Heaviside(t+1)-Heaviside(t-1); W1=2*pi*5;

N=500;k=0:N;W=k*W1/N; F=f*exp(-j*t'*W)*R; F=real(F);

W=[-fliplr(W),W(2:501)]; F=[fliplr(F),F(2:501)]; subplot(2,1,1);plot(t,f); xlabel('t');ylabel('f(t)'); title('f(t)=u(t+1)-u(t-1)'); subplot(2,1,2);plot(W,F); xlabel('w');ylabel('F(w)'); title('f(t)的付氏变换F(w)');

3．用MATLAB分析LTI系统的频率特性

当系统的频率响应H（jw）是jw的有理多项式时，有

jwM)?1?L?b1jw(?)b0B(w)bM(jw)M?bM?1( H(jw)? ?A(w)aN(jw)N?aN?1(jw)N?1?L?a1(jw)?a0MATLAB信号处理工具箱提供的freqs函数可直接计算系统的频率响应的数值解。其调用格式如下

H=freqs(b,a,w)

其中，a和b分别是H(jw)的分母和分子多项式的系数向量，w为形如w1:p:w2的向量，定义系统频率响应的频率范围，w1为频率起始值，w2为频率终止值,p为频率取样间隔。H返回w所定义的频率点上，系统频率响应的样值。

例如，运行如下命令，计算0~2pi频率范围内以间隔0.5取样的系统频率响应的样值 a=[1 2 1]; b=[0 1];

h=freqs(b,a,0:0.5:2*pi)

例 4-5 三阶归一化的butterworth 低通滤波器的频率响应为

H(jw)?1 32(jw)?2(jw)?2(jw)?1 试画出该系统的幅度响应H(jw)和相位响应?(?)。 解 其MATLAB程序及响应的波形如下 w=0:0.025:5; b=[1];a=[1,2,2,1]; H=freqs(b,a,w); subplot(2,1,1);

plot(w,abs(H));grid; xlabel('\\omega(rad/s)');

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

ylabel('|H(j\\omega)|'); title('H(jw)的幅频特性'); subplot(2,1,2);

plot(w,angle (H));grid; xlabel('\\omega(rad/s)'); ylabel('\\phi(\\omega)'); title('H(jw)的相频特性');

H(jw)的幅频特性1|H(j?)|0.5000.511.52.53?(rad/s)H(jw)的相频特性23.544.5542?(?)0-2-400.511.522.5?(rad/s)33.544.55

4．用MATLAB分析LTI系统的输出响应

例 4-6已知一RC电路如图所示 系统的输入电压为f(t)，输出信号为电阻两端的电压y(t).当RC=0.04，f(t)=cos5t+cos100t, ???t??? 试求该系统的响应y(t)

+ f(t) - C R + y(t) -

解 由图可知 ，该电路为一个微分电路，其频率响应为 H(jw)?Rjw?

R?1jwCjw?1RC乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

由此可求出余弦信号cos?0t通过LTI系统的响应为

y(t)?H(j0w)co?s0(?t??0( ))计算该系统响应的MATLAB程序及响应波形如下

RC=0.04;

t=linspace(-2,2,1024); w1=5;w2=100;

H1=j*w1/(j*w1+1/RC); H2=j*w2/(j*w2+1/RC); f=cos(5*t)+cos(100*t);

y=abs(H1)*cos(w1*t+angle(H1))+ abs(H2)*cos(w2*t+angle(H2)); subplot(2,1,1); plot(t,f); ylabel('f(t)'); xlabel('Time(s)'); subplot(2,1,2); plot(t,y); ylabel('y(t)'); xlabel('Time(s)');

21f(t)0-1-2-2-1.5-1-0.50Time(s)0.511.5221y(t)0-1-2-2-1.5-1-0.50Time(s)0.511.52

三、实验内容

1.验证实验原理中所述的相关程序；

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

2.设H(jw)?1，试用MATLAB画出该系统的幅频特性H(jw)和相20.08(jw)?0.4jw?1频特性?(?)。

四、 实验报告要求

1、 简述实验目的及实验原理

2、 记录连续系统的幅频特性和相频特性曲线，分析系统作用 3、 写出程序清单 4、 收获与建议

实验五 连续信号的采样与恢复

一． 实验目的

通过MATLAB仿真验证抽样定理，进一步加深同学们对抽样定理的理解。 二． 实验原理

1、 连续信号的采样

对某一连续时间信号f(t)的采样原理图为：

由图可知，fs(t)?f(t)??TS(t)，其中，单位冲激采样信号?TS(t)的表达式为：

?T(t)?Sn?????(t?nT)

s?其傅里叶变换为?sn?????(??n??s)，其中?s?2?T,设F(j?)为f(t)的傅里叶变换，

fs(t)的频谱为Fs(j?)，由傅里叶变换的频域卷积定理，有：

?11fs(t)?f(t)??TS(t)?FS(j?)?F(j?)??s??(??n?s)?2?TSn???n????F[j(??n?)]s? 若设f(t)是带限信号，带宽为?m，即当???m时，f(t)的频谱F(j?)?0，则f(t)乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

经过采样后的频谱FS(j?)就是FS(j?)在频率轴上搬移至0,??s,?2?s,...,?n?s,...处（幅度为原频谱的1Ts倍）。因此，当?s?2?m时，频谱不会发生混叠；而当?s?2?m时，频谱发生混叠。

2、 连续信号的恢复

设信号f(t)被采样后形成的采样信号为fs(t)，信号的重构是指由fs(t)经内插处理后，恢复出原来的信号f(t)的过程，因此又称为信号恢复。设f(t)为带限信号，带宽为?m，经采样后的频谱为FS(j?)。设采样频谱?s?2?m，则Fs(j?)是以?s为周期的谱线。现取一个频率特性为

??TS,???CH(j?)??(其中，截止频率?c满足?m??c??s2)

??0,???C的理想低通滤波器与FS(j?)相乘，得到的频谱即为原信号的频谱F(j?)。 根据时域卷积定理，有：

f(t)?h(t)?fs(t)

fs(t)?f(t)?n?????(t?nT)??f(nT)??(t?nT)sssn?????其中，h(t)?F?1[H(j?)]?Ts?c?Sa(?ct)

?c为H(j?)的截止角频率。因此，得到：

f(t)?fs(t)?Ts?c?Sa(?ct)?Ts?c?n????f(nT)Sa[?s?c(t?nTs)]

上式即为用f(nTs)表达f(t)的表达式，其中的抽样函数Sa(?ct)为内插函数。 3、 应用举例

例5-1：下面我们选取信号f(t)?Sa(t)?sin(t)/t作为采样的信号，其：

???,??1F(j?)??

??0,??1 即信号的带宽

?m?1。当采样频率?s?2?m时，被称为临界采样（取?c??m）

。在

临界采样状态下实现对信号Sa(t)的采样及由该采样信号恢复Sa(t)的参考程序如下：

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

clear;

wm=1; ％信号带宽 wc=wm; ％滤波器截止频率 Ts=pi/wm; ％采样间隔 ws=2*pi/Ts; ％采样角频率 n=-100:100; %时域采样点数 nTs=n*Ts; ％时域采样点 f=sinc(nTs/pi); Dt=0.005; t=-15:Dt:15;

fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); ％信号重构 error=abs(fa-sinc(t/pi)); ％求重构信号与原信号的误差 t1=-15:0.5:15; f1=sinc(t1/pi); subplot(3,1,1); stem(t1,f1); xlabel('kTs'); ylabel('f(kTs)');

title('sa(t)=sinc(t/pi)临界采样信号'); subplot(3,1,2); plot(t,fa); xlabel('t'); ylabel('fa(t)');

title('由sa(t)=sinc(t/pi)的临界采样信号重构sa(t)'); grid;

subplot(3,1,3); plot(t,error); xlabel('t');

ylabel('error(t)');

title('临界采样信号与原信号的误差error(t)');

MATLAB程序代码运行结果如下：

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

三． 实验内容

设信号f(t)?Sa(t)?sin(t)/t，在取样间隔分别为Ts?0.7?（令

?m?1,?c?1.1?m）和Ts?1.5?（令?m?1,?c??m）的两种情况下，对信号f(t)进

行采样，试编写MATLAB程序代码，并绘制出采样信号波形、由采样信号得到的重构信号波形以及两信号的绝对误差波形。 四． 实验报告要求

1、 简述实验目的及实验原理 2、 写出程序清单

3、 记录实验结果，并进行原理阐述 4、 收获与建议

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

实验六

一、实验目的

连续时间信号与系统的复频域分析

1.学会用MATLAB进行部分分式展开； 2.学会用MATLAB分析LTI系统的特性； 3.学会用MATLAB进行Laplace正、反变换。

二、实验原理

1．用MATLAB进行部分分式展开

用MATLAB函数residue可以得到复杂有理分式F(s)的部分分式展开式，其调用格式为 ?r,p,k??residue(num,den)

其中，num,den分别为F(s)的分子和分母多项式的系数向量，r为部分分式的系数，p为极点，k为F(s)中整式部分的系数，若F(s)为有理真分式，则k为零。

例6-1 用部分分式展开法求F(s)的反变换 F(s)?s?2

s3?4s2?3s解:其MATLAB程序为 format rat; num=[1,2]; den=[1,4,3,0];

[r,p]=residue(num,den)

程序中format rat是将结果数据以分数形式显示

?1?0.56 ?F(s)可展开为 F(s)?3?ss?1s?3所以，F(s)的反变换为 f(t)???e?e6?322?21?t1?3t?ut( )??2．用MATLAB分析LTI系统的特性

系统函数H（s）通常是一个有理分式，其分子和分母均为多项式。计算H（s）的零极点可以应用MATLAB中的roots函数，求出分子和分母多项式的根，然后用plot命令画图。

在MATLAB中还有一种更简便的方法画系统函数H（s）的零极点分布图，即用pzmap函数画图。其调用格式为

pzmap(sys)

sys表示LTI系统的模型，要借助tf函数获得，其调用格式为

sys=tf(b,a)

式中，b和a分别为系统函数H（s）的分子和分母多项式的系数向量。

（j?）如果已知系统函数H（s），求系统的单位冲激响应h(t)和频率响应H可以用以前

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

介绍过的impulse和freqs函数。 例6-2 已知系统函数为 H(s)=1 s3?2s2?2s?1（j?）试画出其零极点分布图，求系统的单位冲激响应h(t)和频率响应H，并判断系统

是否稳定。

解：其MATLAB程序如下： num=[1]; den=[1,2,2,1]; sys=tf(num,den); figure(1);pzmap(sys); t=0:0.02:10;

h=impulse(num,den,t); figure(2);plot(t,h)

title('Impulse Response') [H,w]=freqs(num,den); figure(3);plot(w,abs(H)) xlabel('\\omega')

title('Magnitude Response')

3．用MATLAB进行Laplace正、反变换

MATLAB的符号数学工具箱提供了计算Laplace正、反变换的函数Laplace和ilaplace,其调用格式为

F?laplace(f)

f?ilaplace(F)上述两式右端的f和F分别为时域表示式和s域表示式的符号表示，可以应用函数sym实现，其调用格式为

S=sym(A)

式中，A为待分析表示式的字符串，S为符号数字或变量。 例6-3 试分别用Laplace和ilaplace函数求 （1）f(t)?esin(at)u(t)的Laplace变换；

?ts2（2）F(s)?2的Laplace反变换。

s?1解：（1）其程序为

f=sym('exp(-t)*sin(a*t)'); F=laplace(f) 或

syms a t

F=laplace(exp(-t)*sin(a*t)) （2）其程序为

F=sym('s^2/(s^2+1)'); ft=ilaplace(F)

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

或 syms s

ft= ilaplace(s^2/(s^2+1)) 三、实验内容

1．验证实验原理中所述的相关程序； 2．求信号f(t)?te?3tu(t)的拉普拉斯变换

s3?5s2?9s?73．求函数F(s)?的反变换

s2?3s?24．已知连续系统的系统函数如下，试用MATLAB绘制系统的零极点图，并根据零极点图判断系统的稳定性

s2?s?2 H(s)=3

3s?5s2?4s?6四、 实验报告要求

1、 简述实验目的及实验原理

2、 计算相应S变换或反S变换的理论值，并与实验结果进行比较 3、 记录连续系统的频率响应特性曲线，分析系统作用 4、 写出程序清单 5、 收获与建议

实验七 离散时间信号与系统的Z域分析

一、实验目的

1.学会用MATLAB画离散系统零极点图； 2.学会用MATLAB分析离散系统的频率特性；

二、实验原理

1．离散系统零极点图

离散系统可以用下述差分方程描述：

?ay(k?i)??bii?0m?0NMmf(k?m)

Y(z)b0?b1z?1?...?bMz?MZ变换后可得系统函数：H(z)? ?F(z)a0?a1z?1?...?aNz?N用MATLAB提供的root函数可分别求零点和极点，调用格式是

p=[a0,a1…an],q=[b0,b1…bm,0,0…0], 补0使二者维数一样。画零极点图的方法有多种，可以用MATLAB函数[z,p,k]=tf2zp(b,a)和zplane(q,p)，也可用plot命令自编一函数ljdt.m,画图时调用。

function ljdt(A,B)

% The function to draw the pole-zero diagram for discrete system

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

p=roots(A)； %求系统极点 q=roots(B)； %求系统零点 p=p'； %将极点列向量转置为行向量 q=q'； %将零点列向量转置为行向量 x=max(abs([p q 1]))； %确定纵坐标范围 x=x+0.1； y=x； %确定横坐标范围 clf hold on

axis([-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,'jIm[z]') text(y,1/10,'Re[z]')

plot(real(p),imag(p),'x') %画极点 plot(real(q),imag(q),'o') %画零点 title('pole-zero diagram for discrete system') %标注标题 hold off 例7-4 求系统函数零极点图H(z)? a=[3 -1 0 0 0 1]; b=[1 1]; ljdt(a,b) p=roots(a) q=roots(b) pa=abs(p)

2．离散系统的频率特性

离散系统的频率特性可由系统函数求出，既令z?ej?z?1 543z?z?1,MATLAB函数freqz可计算频率

特性，调用格式是：

[H，W]=freqz(b,a,n)，b和a是系统函数分子分母系数，n是0-?范围 n个等份点，默认

值512，H是频率响应函数值，W是相应频率点；

[H，W]=freqz(b,a,n,’whole’), n是0-2?范围 n个等份点； freqz(b,a,n)，直接画频率响应幅频和相频曲线； 例7-5 系统函数H(z)?z?0.5 运行如下语句，可得10个频率点的计算结果

乐山师范学院.物理与电子工程学院

z

《MATLAB&信号与系统》实验指导书

A=[1 0]; B=[1 -0.5];

[H,W]=freqz(B,A,10)

继续运行如下语句，可将400个频率点的计算结果用plot语句画幅频和相频曲线 B=[1 -0.5]; A =[1 0];

[H,w]=freqz(B,A,400,'whole'); Hf=abs(H); Hx=angle(H); clf

figure(1) plot(w,Hf)

title('离散系统幅频特性曲线') figure(2) plot(w,Hx)

title('离散系统相频特性曲线')

还可用freqz语句直接画图，注意区别 A=[1 0]; B=[1 -0.5]; freqz(B,A,400)

例7-6 用几何矢量法，自编程序画频率响应.

j?(e??qj)M原理：频率响应H(ej?)?j?1N?(e??p)jii?1

编程流程：定义Z平面单位圆上k个频率等分点；求出系统函数所有零点和极点到这些等

分点的距离；求出系统函数所有零点和极点到这些等分点的矢量的相角；求出单位圆上各 频率等分点的

H(ej?)和?(?)画指定范围内的幅频与相频。若要画零极点图，可调用ljdt.m函数。 function dplxy(k,r,A,B)

%The function to draw the frequency response of discrete system p=roots(A); %求极点 q=roots(B); %求零点 figure(1) ljdt(A,B) %画零极点图 w=0:l*pi/k:r*pi; y=exp(i*w); %定义单位圆上的k个频率等分点 N=length(p); %求极点个数 M=length(q); %求零点个数 yp=ones(N,1)*y; %定义行数为极点个数的单位圆向量

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

yq=ones(M,1)*y; %定义行数为零点个数的单位圆向量 vp=yp-p*ones(1,r*k+1); %定义极点到单位圆上各点的向量 vq=yq-q*ones(1,r*k+1); %定义零点到单位圆上各点的向量 Ai=abs(vp); %求出极点到单位圆上各点的向量的模 Bj=abs(vq); %求出零点到单位圆上各点的向量的模 Ci=angle(vp); %求出极点到单位圆上各点的向量的相角 Dj=angle(vq); %求出零点到单位圆上各点的向量的相角 fai=sum(Dj,1)-sum(Ci,1); %求系统相频响应 H=prod(Bj,1)./prod(Ai,1); %求系统幅频响应 figure(2) plot(w,H); %绘制幅频特性曲线 title('离散系统幅频特性曲线') xlabel('角频率') ylabel('幅度') figure(3) plot(w,fai) title('离散系统的相频特性曲线') xlabel('角频率') ylabel('相位') 5/4(1?z?1)例7-7已知系统函数H(z)?，画频率响应和零极点图。 ?11?1/4zA=[1 -1/4]; B=[5/4 -5/4];

dplxy(500,2,A,B) %绘制系统2π频率范围内500个频率点的幅频和相频特性曲线及零极

点图

三、实验内容

1、已知离散系统的系统函数如下所示：

z2?2H(z)?3 2z?2z?4z?1试用MATLAB实现下列分析过程： （1） 求出系统的零极点位置；

（2） 绘出系统的零极点图，根据零极点图判断系统的稳定性；

（3） 绘出系统单位响应的时域波形，并分析系统稳定性与系统单位响应时域特性的

关系。

2、已知描述离散系统的差分方程为：

y(k)?y(k?1)?y(k?2)?4f(k)?f(k?1)?f(k?2)

试用MATLAB绘出系统的零极点分布图，并绘出系统的幅频和相频特性曲线，分析该系统的作用。

3、已知因果（单边）离散序列的Z变换如下，试用MATLAB求出其逆Z变换。

乐山师范学院.物理与电子工程学院

《MATLAB&信号与系统》实验指导书

2z2?z?1F(z)?

1z3?z2?z2

四、 实验报告要求

1、 简述实验目的及实验原理。 2、 计算相应z变换或反z变换的理论值，并与实验结果进行比较。 3、 记录离散系统的频率响应特性曲线，分析系统作用。 4、 写出程序清单。 5、 收获与建议。

实验八 线性系统的特性分析和输出响应求解

（综合设计性）

一、实验目的

该实验在常规实验的基础上，综合运用“信号与系统”的理论知识，结合MATLAB软件仿真、分析线性系统的特性和系统的响应。通过该实验，同学们初步掌握综合运用理论知识、软件仿真进行简单系统的设计与分析的基本方法，为同学们创造应用实践契机，培养同学们的自学能力和钻研能力。

二、实验内容

设系统的频率特性H(j?)为：

H(j?)?24

(j??1)(j??4)(j??6)（1） 用MATLAB画出该系统的单位冲激响应和单位阶跃响应曲线；

（2） 设输入信号为f(t)?sin(3t)?sin(90t),???t??,，画出f(t)通过该系统的响应

y(t)的波形；

（3） 用MATLAB画出该系统的幅频特性和相频特性；

（4） 若系统输入信号为f(t)?e，试用MATLAB求输出信号Y(s)的部分分式展开式，

并写出其时域表达式y(t)。

三、基本要求

（1） 阅读教师给出的参考资料；

（2） 通过查找纸质资料和网络资料，独立地提出研究问题的设计方案； （3） 编程实现系统特性的仿真分析和响应的求解； （4） 综合理论基础、软件仿真进行分析与讨论； （5） 最后撰写出实验报告。

乐山师范学院.物理与电子工程学院

?t

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