信号与系统及MATLAB实现.doc

信号与系统实验指导书实验一 基本信号在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 delta(k1,k2,k0) % k1 , k2 是画图时间范围，k0是脉冲位置（k1=k0 =k2）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);保存文件名delta.m画图时在命令窗口调用，例：delta (-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(f2(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