信号与系统及MATLAB实现

信号与系统及MATLAB实现实验指导书前 言长期以来，信号与系统课程一直采用单一理论教学方式，同学们依靠做习题来巩固和理解教学内容，虽然手工演算训练了计算能力和思维方法，但是由于本课程数学公式推导较多，概念抽象，常需画各种波形，作题时难免花费很多时间，现在，我们给同学们介绍一种国际上公认的优秀科技应用软件MATLAB，借助它我们可以在电脑上轻松地完成许多习题的演算和波形的绘制。 MATLAB的功能非常强大，我们此处仅用到它的一部分，在后续课程中我们还会用到它，在未来地科学研究和工程设计中有可能继续用它，所以有兴趣的同学，可以对MATLAB再多了解一些。MATLAB究竟有那些特点呢？1高效的数值计算和符号计算功能，使我们从繁杂的数学运算分析中解脱出来；2完备的图形处理功能，实现计算结果和编程的可视化；3友好的用户界面及接近数学表达式的自然化语言，易于学习和掌握；4功能丰富的应用工具箱，为我们提供了大量方便实用的处理工具；MATLAB的这些特点，深受大家欢迎，由于个人电脑地普及，目前许多学校已将它做为本科生必须掌握的一种软件。正是基于这些背景，我们编写了这本信号与系统及MATLAB实现指导书，内容包括信号的MATLAB表示、基本运算、系统的时域分析、频域分析、S域分析、状态变量分析等。通过这些练习，同学们在学习信号与系统的同时，掌握MATLAB的基本应用，学会应用MATLAB的数值计算和符号计算功能，摆脱烦琐的数学运算，从而更注重于信号与系统的基本分析方法和应用的理解与思考，将课程的重点、难点及部分习题用MATLAB进行形象、直观的可视化计算机模拟与仿真实现，加深对信号与系统的基本原理、方法及应用的理解，为学习后续课程打好基础。另外同学们在进行实验时，最好事先预习一些MATLAB的有关知识，以便更好地完成实验，同时实验中也可利用MATLAB的help命令了解具体语句以及指令的使用方法。实验一 基本信号在MATLAB中的表示和运算一、实验目的1 学会用MATLAB表示常用连续信号的方法；2 学会用MATLAB进行信号基本运算的方法；二、实验原理1 连续信号的MATLAB表示MATLAB提供了大量的生成基本信号的函数，例如指数信号、正余弦信号。表示连续时间信号有两种方法，一是数值法，二是符号法。数值法是定义某一时间范围和取样时间间隔，然后调用该函数计算这些点的函数值，得到两组数值矢量，可用绘图语句画出其波形；符号法是利用MATLAB的符号运算功能，需定义符号变量和符号函数，运算结果是符号表达的解析式，也可用绘图语句画出其波形图。例1-1指数信号 指数信号在MATLAB中用exp函数表示。如，调用格式为 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 函数表示。定义为 t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on; axis(-10,10,-0.5,1.2); %定义画图范围，横轴，纵轴 title(抽样信号) %定义图的标题名字例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); plot(t,ft); grid on; axis(-1,5,-0.5,1.5);例1-9 正弦信号符号算法syms t %定义符号变量ty=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=(t0); 保存，文件名是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 以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);例1-12 已知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)符号算法也可实现上述运算，以信号的微积分运算为例说明符号算法应用微分的调用格式为 diff(function,variable,n)积分的调用格式为 int(function,variable,a,b)式中function表示要微分或积分的函数，variable表示运算变量，n表示求导阶数，默认值是求一阶导数，a是积分下限，b是积分上限，a b默认是求不定积分。例1-13 求一阶导数的例题，已知，clear syms a x y1 y2 %定义符号变量a， x ，y1， y2 y1=sin(a*x2); %符号函数y1 y2=x*sin(x)*log(x); %符号函数y2 dy1=diff(y1,x) %无分号直接显示结果 dy2=diff(y2) %无分号直接显示结果例1-14 求积分的例题， clearsyms a x y3 y4y3=x5-a*x2+sqrt(x)/2;y4=(x*exp(x)/(1+x)2;iy3=int(y3,x)iy4=int(y4,0,1)三、上机实验内容1 验证实验原理中程序2 画出信号波形（1） （2）3信号，求、波形实验二 离散信号与系统的时域分析一、实验目的1学会用MATLAB表示常用离散信号的方法；2学会用MATLAB实现离散信号卷积的方法；3 学会用MATLAB求解离散系统的单位响应；4 学会用MATLAB求解离散系统的零状态响应；二、实验原理1离散信号的MATLAB表示表示离散时间信号f(k)需要两个行向量，一个是表示序号k= ，一个是表示相应函数值f= ，画图命令是stem。例2-1正弦序列信号 正弦序列信号可直接调用MATLAB函数cos，例，当是整数或分数时，才是周期信号。画，波形程序是：k=0:40;subplot(2,1,1)stem(k,cos(k*pi/8),filled)title(cos(k*pi/8)subplot(2,1,2)stem(k,cos(2*k),filled)title(cos(2*k) 例2-2 单位序列信号 本题先建立一个画单位序列的M函数文件，画图时调用。M文件建立方法：file / new / m-file 在文件编辑窗输入程序，保存文件名用函数名。function dwxulie(k1,k2,k0) % k1 , k2 是画图时间范围，k0是脉冲位置k=k1:k2;n=length(k);f=zeros(1,n);f(1,-k0-k1+1)=1; stem(k,f,filled) axis(k1,k2,0,1.5)title(单位序列(k) 保存文件名dwxulie.m画图时在命令窗口调用，例：dwxulie(-5,5,0)例2-3 单位阶跃序列信号 本题也可先建立一个画单位阶跃序列的M函数文件，画图时调用。 function jyxulie(k1,k2,k0)k=k1:-k0-1;kk=-k0:k2;n=length(k);nn=length(kk)u=zeros(1,n);uu=ones(1,nn); stem(kk,uu,filled)hold onstem(k,u,filled)hold offtitle(单位阶跃序列)axis(k1 k2 0 1.5)保存文件名jyxulie.m画图时在命令窗口调用，例：jyxulie(-3,8,0)例2-4 实指数序列信号 ，c、 a是实数。建立一个画实指数序列的M函数文件，画图时调用。function dszsu(c,a,k1,k2)%c：指数序列的幅度%a：指数序列的底数%k1：绘制序列的起始序号%k2：绘制序列的终止序号k=k1:k2;x=c*(a.k);stem(k,x,filled)hold onplot(k1,k2,0,0)hold off调用该函数画信号： ，波形。dszsu(1,5/4,0,40)dszsu(1,-3/4,0,40)2 .离散信号的卷积和两个有限长序列f1，f2卷积可调用MATLAB函数conv，调用格式是f=conv(f1,f2), f是卷积结果，但不显示时间序号，可自编一个函数dconv给出f和k，并画图。function f,k=dconv(f1,f2,k1,k2) %The function of compute f=f1*f2% f: 卷积和序列f(k)对应的非零样值向量% k： 序列f(k)的对应序号向量% f1: 序列f1(k)非零样值向量% f2: 序列f2(k)的非零样值向量% k1: 序列f1(k)的对应序号向量% k2: 序列f2(k)的对应序号向量f=conv(f1,f2) %计算序列f1与f2的卷积和fk0=k1(1)+k2(1);%计算序列f非零样值的起点位置k3=length(f1)+length(f2)-2;%计算卷积和f的非零样值的宽度k=k0:k0+k3 %确定卷积和f非零样值的序号向量subplot(2,2,1)stem(k1,f1)%在子图1绘序列f1(k)时域波形图title(f1(k)xlabel(k)ylabel(f1(k)subplot(2,2,2)stem(k2,f2)%在图2绘序列f2(k)时波形图title(f1(k)xlabel(k)ylabel(f2(k)subplot(2,2,3)stem(k,f);%在子图3绘序列f(k)的波形图title(f(k)f1(k)与f2(k)的卷积和f(k)xlabel(k)ylabel(f(k)h=get(gca,position);h(3)=2.5*h(3);set(gca,position,h)%将第三个子图的横坐标范围扩为原来的2.5倍例2-5求卷积和， f1=1 2 1;k1=-1 0 1;f2=ones(1,5);k2=-2:2;f, k=dconv(f1,f2,k1,k2) 由运行结果知，f的长度等于f1和f2长度之和减一， f的起点是f1和f2的起点之和，f的终点是f1和f2的终点之和。3 离散系统的单位响应MATLAB提供画系统单位响应函数impz，调用格式是impz(b,a) 式中b和a是表示离散系统的行向量；impz(b,a,n) 式中b和a是表示离散系统的行向量，时间范围是0n；impz(b,a,n1,n2) 时间范围是n1n2 ；y=impz(b,a,n1,n2) 由y给出数值序列；例2-6已知 求单位响应。a=1,-1,0.9;b=1;impz(b,a)impz(b,a,60)impz(b,a,-10:40)4 离散系统的零状态响应MATLAB提供求离散系统零状态响应数值解函数filter，调用格式为filter(b,a,x),式中b和a是表示离散系统的向量，x是输入序列非零样值点行向量，输出向量序号同x一样。例2-7 已知 ， 求零状态响应， 范围020。a=1 -0.25 0.5;b=1 1;t=0:20;x=(1/2).t;y=filter(b,a,x)subplot(2,1,1)stem(t,x)title(输入序列)subplot(2,1,2)stem(t,y)title(响应序列)三、上机实验内容1验证实验原理中程序2已知，画单位响应波形。3已知，输入，画输出波形，范围015。实验三 连续时间LTI系统的时域分析一、实验目的1学会用MATLAB求解连续系统的零状态响应；2. 学会用MATLAB求解冲激响应及阶跃响应； 3学会用MATLAB实现连续信号卷积的方法；二、实验原理1连续时间系统零状态响应的数值计算 我们知道，LTI连续系统可用如下所示的线性常系数微分方程来描述， 在MATLAB中，控制系统工具箱提供了一个用于求解零初始条件微分方程数值解的函数lsim。其调用格式y=lsim(sys,f,t)式中，t表示计算系统响应的抽样点向量，f是系统输入信号向量，sys是LTI系统模型，用来表示微分方程，差分方程或状态方程。其调用格式sys=tf(b,a)式中，b和a分别是微分方程的右端和左端系数向量。例如，对于以下方程: 可用 获得其LTI模型。注意，如果微分方程的左端或右端表达式中有缺项，则其向量a或b中的对应元素应为零，不能省略不写，否则出错。例3-1 已知某LTI系统的微分方程为 y(t)+ 2y(t)+100y(t)=f(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来求解。其调用格式为 y=impluse(sys,t) y=step(sys,t)式中，t表示计算系统响应的抽样点向量，sys是LTI系统模型。例3-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( )函数近似计算信号的卷积积分。连续信号的卷积积分定义是 如果对连续信号和进行等时间间隔均匀抽样，则和分别变为离散时间信号和。其中，为整数。当足够小时，和既为连续时间信号和。因此连续时间信号卷积积分可表示为 采用数值计算时，只求当时卷积积分的值，其中，n为整数，既 其中，实际就是离散序列和的卷积和。当足够小时，序列就是连续信号的数值近似，既 上式表明，连续信号和的卷积，可用各自抽样后的离散时间序列的卷积再乘以抽样间隔。抽样间隔越小，误差越小。例3-3用数值计算法求与的卷积积分。解：因为是一个持续时间无限长的信号，而计算机数值计算不可能计算真正的无限长信号，所以在进行的抽样离散化时，所取的时间范围让衰减到足够小就可以了，本例取。程序是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)由于和的时间范围都是从t=-1开始，所以卷积结果的时间范围从 t=-2开始，增量还是取样间隔，这就是语句tt=(0:n-1)*dt-2的由来。三、上机实验内容1 验证实验原理中所述的相关程序2 已知描述系统的微分方程和激励信号f(t)如下，试用解析法求系统的零状态响应y(t)，并用MATLAB绘出系统零状态响应的时域仿真波形，验证结果是否相同y(t)+ 4y(t)+4y(t)=f(t)+3f(t) f(t)= exp(-t)u(t) 3已知描述系统的微分方程如下，试用MATLAB求系统在010秒范围内冲激响应和阶跃响应的数值解，并用绘出系统冲激响应和阶跃响应的时域波形y(t)+3y(t)+2y(t)=f(t)y(t)+ 2y(t)+2y(t)=f(t) 4画出信号卷积积分的波形，实验四 傅里叶变换、系统的频域分析一、 实验目的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，即 Fourier逆变换的调用格式f=ifourier(F):它是符号函数F的fourier逆变换，默认的独立变量为w，默认返回是关于x的函数。f=ifourier(f,u):它的返回函数f是u的函数，而不是默认的x.注意：在调用函数fourier()及ifourier()之前，要用syms命令对所用到的变量（如t,u,v,w）进行说明,即将这些变量说明成符号变量。例4-1 求的傅立叶变换解: 可用MATLAB解决上述问题：syms tFw=fourier(exp(-2*abs(t)例4-2 求的逆变换f(t)解： 可用MATLAB解决上述问题syms t wft=ifourier(1/(1+w2),t)2连续时间信号的频谱图例4-3 求调制信号的频谱，式中解：MATLAB程序如下所示ft=sym(4*cos(2*pi*6*t)*(Heaviside(t+1/4)-Heaviside(t-1/4);Fw=simplify(fourier(ft)subplot(121)ezplot(ft,-0.5 0.5),grid onsubplot(122)ezplot(abs(Fw),-24*pi 24*pi),grid 用MATLAB符号算法求傅里叶变换有一定局限，当信号不能用解析式表达时，会提示出错，这时用MATLAB的数值计算也可以求连续信号的傅里叶变换，计算原理是当足够小时，近似计算可满足要求。若信号是时限的，或当时间大于某个给定值时，信号已衰减的很厉害，可以近似地看成时限信号时，n的取值就是有限的，设为N，有 是频率取样点时间信号取样间隔应小于奈奎斯特取样时间间隔，若不是带限信号可根据计算精度要求确定一个频率 W0为信号的带宽。例4-4 用数值计算法求信号的傅里叶变换解，信号频谱是，第一个过零点是，一般将此频率视为信号的带宽，若将精度提高到该值的50倍，既W0=50，据此确定取样间隔，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的有理多项式时，有 MATLAB信号处理工具箱提供的freqs函数可直接计算系统的频率响应的数值解。其调用格式如下H=freqs(b,a,w)其中，a和b分别是H(jw)的分母和分子多项式的系数向量，w为形如w1:p:w2的向量，定义系统频率响应的频率范围，w1为频率起始值，w2为频率终止值,p为频率取样间隔。H返回w所定义的频率点上，系统频率响应的样值。例如，运行如下命令，计算02pi频率范围内以间隔0.5取样的系统频率响应的样值a=1 2 1;b=0 1;h=freqs(b,a,0:0.5:2*pi)例 4-5 三阶归一化的butterworth 低通滤波器的频率响应为 试画出该系统的幅度响应和相位响应。解 其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);ylabel(|H(jomega)|);title(H(jw)的幅频特性);subplot(2,1,2);plot(w,angle (H);grid;xlabel(omega(rad/s);ylabel(phi(omega);title(H(jw)的相频特性);4用MATLAB分析LTI系统的输出响应例 4-6已知一RC电路如图所示 系统的输入电压为f(t)，输出信号为电阻两端的电压y(t).当RC=0.04，f(t)=cos5t+cos100t, 试求该系统的响应y(t)-+-+f(t)y(t)RC解 由图可知 ，该电路为一个微分电路，其频率响应为 由此可求出余弦信号通过LTI系统的响应为 计算该系统响应的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);三、 上机实验内容1.验证实验原理中所述的相关程序；2.试用MATLAB求单边指数数信号的傅立叶变换，并画出其波形；3.设，试用MATLAB画出该系统的幅频特性和相频特性，并分析系统具有什么滤波特性。实验五 信号抽样与恢复一、实验目的学会用MATLAB实现连续信号的采样和重建二、实验原理 1抽样定理若是带限信号，带宽为, 经采样后的频谱就是将的频谱 在频率轴上以采样频率为间隔进行周期延拓。因此，当时，不会发生频率混叠；而当 时将发生频率混叠。2信号重建经采样后得到信号经理想低通则可得到重建信号，即：=*其中：=所以：=*=* =上式表明，连续信号可以展开成抽样函数的无穷级数。利用MATLAB中的来表示，有 ，所以可以得到在MATLAB中信号由重建的表达式如下：=我们选取信号=作为被采样信号，当采样频率=2时，称为临界采样。我们取理想低通的截止频率=。下面程序实现对信号=的采样及由该采样信号恢复重建：例5-1 Sa(t)的临界采样及信号重构；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); %信号重构t1=-15:0.5:15;f1=sinc(t1/pi);subplot(211);stem(t1,f1);xlabel(kTs);ylabel(f(kTs);title(sa(t)=sinc(t/pi)的临界采样信号);subplot(212);plot(t,fa)xlabel(t);ylabel(fa(t);title(由sa(t)=sinc(t/pi)的临界采样信号重构sa(t);grid;例5-2 Sa(t)的过采样及信号重构和绝对误差分析程序和例4-1类似，将采样间隔改成Ts=0.7*pi/wm , 滤波器截止频率该成wc=1.1*wm ，添加一个误差函数wm=1;wc=1.1*wm; Ts=0.7*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Tsf=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(311);stem(t1,f1);xlabel(kTs);ylabel(f(kTs);title(sa(t)=sinc(t/pi)的采样信号);subplot(312);plot(t,fa)xlabel(t);ylabel(fa(t);title(由sa(t)=sinc(t/pi)的过采样信号重构sa(t);grid;subplot(313);plot(t,error);xlabel(t);ylabel(error(t);title(过采样信号与原信号的误差error(t);例5-3 Sa(t)的欠采样及信号重构和绝对误差分析程序和例4-2类似，将采样间隔改成Ts=1.5*pi/wm , 滤波器截止频率该成wc=wm=1三、上机实验内容1验证实验原理中所述的相关程序；2设f(t)=0.5*(1+cost)*(u(t+pi)-u(t-pi) ，由于不是严格的频带有限信号，但其频谱大部分集中在0，2之间，带宽wm可根据一定的精度要求做一些近似。试根据以下两种情况用 MATLAB实现由f(t)的抽样信号fs(t)重建f(t) 并求两者误差，分析两种情况下的结果。(1) wm=2 , wc=1.2wm , Ts=1;(2) wm=2 , wc=2 , Ts=2.5 实验六 信号与系统复频域分析一、实验目的1.学会用MATLAB进行部分分式展开；2.学会用MATLAB分析LTI系统的特性；3.学会用MATLAB进行Laplace正、反变换。4.学会用MATLAB画离散系统零极点图；5.学会用MATLAB分析离散系统的频率特性；二、实验原理及内容1用MATLAB进行部分分式展开用MATLAB函数residue可以得到复杂有理分式F(s)的部分分式展开式，其调用格式为 其中，num,den分别为F(s)的分子和分母多项式的系数向量，r为部分分式的系数，p为极点，k为F(s)中整式部分的系数，若F(s)为有理真分式，则k为零。例6-1 用部分分式展开法求F(s)的反变换 解:其MATLAB程序为format rat;num=1,2;den=1,4,3,0;r,p=residue(num,den)程序中format rat是将结果数据以分数形式显示F(s)可展开为 所以，F(s)的反变换为 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）的分子和分母多项式的系数向量。如果已知系统函数H（s），求系统的单位冲激响应h(t)和频率响应可以用以前介绍过的impulse和freqs函数。例6-2 已知系统函数为 试画出其零极点分布图，求系统的单位冲激响应h(t)和频率响应，并判断系统是否稳定。解：其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和F分别为时域表示式和s域表示式的符号表示，可以应用函数sym实现，其调用格式为 S=sym(A)式中，A为待分析表示式的字符串，S为符号数字或变量。例6-3 试分别用Laplace和ilaplace函数求（1）的Laplace变换；（2）的Laplace反变换。解：（1）其程序为f=sym(exp(-t)*sin(a*t);F=laplace(f)或syms a tF=laplace(exp(-t)*sin(a*t)（2）其程序为F=sym(s2/(s2+1);ft=ilaplace(F)或syms sft= ilaplace(s2/(s2+1)4离散系统零极点图离散系统可以用下述差分方程描述：Z变换后可得系统函数：用MATLAB提供的root函数可分别求零点和极点，调用格式是p=a0,a1an,q=b0,b1bm,0,00, 补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 systemp=roots(A)； %求系统极点q=roots(B)； %求系统零点p=p； %将极点列向量转置为行向量q=q； %将零点列向量转置为行向量x=max(abs(p q 1)；%确定纵坐标范围x=x+0.1；y=x；%确定横坐标范围clfhold onaxis(-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)%画纵坐标轴