MATLAB在信号与系统中的应用举例（课件PPT） 49页

'第9章在信号与系统中的应用 9.1连续信号和系统严格说来，MATLAB（就基本部分而言）是不能表示连续信号的，因为它给出的是各个样本点的数据。只是当样本点取得很密时可看成连续信号，什么叫密，要相对于信号变化的快慢而言，形象地说，在所有相邻样本点之间的数据变化必须非常小才能看成‘密’，其严格的数学定义此处不予讨论。以下均假定相对于采样点密度而言,信号变化足够慢。例9-1-1连续信号的MATLAB描述列出单位脉冲、单位阶跃、复指数函数等连续信号的MATLAB表达式。 程序exn911(1)clear,t0=0;tf=5;dt=0.05;t1=1;t=[t0:dt:tf];%(1)单位脉冲信号,%在t1(t0≤t1≤tf)处有面积为1的脉冲信号。t=[t0:dt:tf];st=length(t);n1=floor((t1-t0)/dt);%求t1对应的样本序号x1=zeros(1,st);%把信号先初始化为零x1(n1)=1/dt;%给出t1处的脉冲信号subplot(2,2,1),stairs(t,x1)%绘图，用stairs命令axis([0,5,0,22])%为了使脉冲顶部避开图框 程序exn911(2)%(2)单位阶跃信号，%信号从t0到tf，在t1前为0，到t1处跃变为1.%程序前几句即求t,st,n1的语句与上同%产生阶跃信号x2=[zeros(1,n1-1),ones(1,st-n1+1)];subplot(2,2,3),stairs(t,x2)%绘图axis([0,5,0,1.1])%改变图框坐标%(3)复数指数信号u=-0.5;w=10;x4=exp((u+j*w)*t);subplot(2,2,2),plot(t,real(x4))%绘图,subplot(2,2,4),plot(t,imag(x4))%绘图, 程序exn911运行结果x1,x2,x3,x4的波形见右图.注意若要显示连续信号波形中的不连续点，用stairs命令；而要使波形光滑些，则用plot命令较好。复数指数信号可以分解为余弦和正弦信号，它们分别是复数信号的实部和虚部。右图中的两个衰减振荡信号就代表了这两个相位差90度的分量。 例9-1-2线性系统冲击响应编写求任意高阶连续常系数线性系统冲击响应的程序。解:◆建模这个问题在第四章介绍多项式函数库时已打下了基础，在第七章例7-3-1中又给出了解法，读者可先看懂那些内容，然后再看本题。任意阶次的连续线性系统可用下列微分方程表述：写成传递函数形式为其特性可用系统传递函数的分子分母系数向量b和a来表示。 传递函数反变换的求法如果分母系数多项式没有重根，则可以把两个多项式之比分解成n个一阶部分分式之和。即：其中p1,p2,…,pn是分母多项式的n个根，而r1,r2,…,rn是则对应于这n个根的留数。一阶分式的反变换可以查表得到，容易写出冲击响应的公式如下：可见只要求出根向量p和留数向量r，线性方程的解就得到了。求根是代数问题，当阶次很高时，代数方程没有解析解。可喜的是MATLAB提供了用数值方法求根和留数的函数residue.m，它的调用格式为：[r,p]=residue(b,a) 程序exn912a=input("分母系数向量a=(书上取poly([0,-1,-2,-5])）")b=input("分子系数向量b=(书上取[1,7,1]）")[r,p]=residue(b,a),%求留数k=input("是否要求波形？是,键入1;否,键入0");ifk==1dt=input("dt=(书上取0.05)");tf=input("tf=(书上取5");%设定时间数组t=0:dt:tf;h=zeros(1,length(t));%h的初始化fori=1:length(a)-1%根数为a的长度减1h=h+r(i)*exp(p(i)*t);%叠加各根分量end,plot(t,h),gridelse,end 程序exn912运行结果给出系统的传递函数为运行程序exn912，依次输入：（注意用poly函数把极点向量p=[0,-1,-2,-5]转换成系数向量a）a=poly([0,-1,-2,-5])B=[1,7,1]，dt=0.05，tf=5得出的h(t)如右图所示。 9-1-3线性系统零输入响应的计算线性时不变连续系统的特性可用常微分方程表示为：求其零输入响应。解：在零输入条件u=0时，等式右端为零。系统响应的通解为：其中，p是特征方程的n个根组成的向量[p1,p2,…,pn]，其每个分量的系数Cn则由y及其各阶导数的初始条件y0,Dy0,…,D(n-1)y0来确定。 代入初始条件得到的矩阵方程初始条件数应该和常数数相等，由此构成一个确定C1,…,Cn的线性代数方程组，写成:矩阵V由特征根向量p确定，这种矩阵称为范德蒙特矩阵。在MATLAB中，有生成它的函数vander(p)。 求零输入响应程序exn913它产生的矩阵与上述矩阵排列转了90度，故用V=rot90(vander(p))，按此思路编成程序exn913:a=input("输入分母系数向量a=[a1,a2,...]=");n=length(a)-1;Y0=input("初始条件Y0=[y0,Dy0,D2y0,...]=");p=roots(a);V=rot90(vander(p))；%生成系数矩阵Vc=VY0‘;%用左除求出系数向量C%以下是计算并画出时间响应的程序段dt=input("dt=");tf=input("tf=")%生成t向量和y初始向量t=0:dt:tf;y=zeros(1,length(t));fork=1:ny=y+c(k)*exp(p(k)*t);endplot(t,y),grid%绘图 数字实例用这一个普遍程序来解一个三阶系统，设其微分方程为：初始条件为：，求零输入响应。解：运行exn913，按提示输入a为[1,2,9,3];Y0为[1,0,0];dt为0.1;tf为5；即可得到t=0~5秒的零输入响应曲线。再分别取Y0为[0,1,0];[0,0,1]，用holdon语句使三次运行生成的图形叠画在一幅图上，得到右图。 例9-1-4重极点求反变换命题：n级放大器，每级的转移函数均为，求阶跃输入下的过渡过程，画出n不同时的波形及频率特性。解：◆建模系统的转移函数为H(s)=，阶跃输入的拉普拉斯变换为1/s，因此输出为两者的乘积，即求Y(s)的拉普拉斯反变换，即可得到输出过渡过程y(t)。这里我们遇到了一个有多重极点-wn的H(s)求拉普拉斯反变换的问题，数学上比较麻烦。为了避开重极点问题，可以有意把极点拉开一些，例如设n个极点散布在-0.98wn到1.02wn之间，由此就可用下面的程序来求： 程序exn914(1)clear,clf,N=input("输入放大器级数N=");wn=1000;dt=1e-4;tf=0.01;t=0:dt:tf;y=zeros(N,length(t));%输出初始化forn=1:Np0=-linspace(.95,1.05,n)*wn;%将H(s)极点分散布置ay=poly([p0,0]);%由Y(s)的极点求分母系数by=prod(abs(p0));%求Y(s)的分子系数[r,p]=residue(by,ay);%求Y(s)的留数极点fork=1:n+1%把各时域分量相加y(n,:)=y(n,:)+r(k)*exp(p(k)*t);endfigure(1),plot(t,y(n,:));grid,holdon%绘制过渡过程end 程序exn914(2)-画波德图%下面的语句用来绘制波德图，如果用bode函数,只要一句%figure(2),bode(prod(abs(p0)),poly(p0));holdonbh=by;ah=poly(p0);%求H(s)的分子分母系数w=logspace(2,4);%给出频率范围和分度H=polyval(bh,j*w)./polyval(ah,j*w);%求H(jw)aH=unwrap(angle(H))*180/pi;%求出相角fH=20*log10(abs(H));%求出振幅figure(2),subplot(2,1,1),semilogx(w,fH),gridon,holdon%绘幅频图subplot(2,1,2),semilogx(w,aH),gridon,holdon%相频图end 程序exn914运行结果(时域)运行此程序，设N=4,可得过渡过程如右图,从中看出输出信号达到0.6处所需的时间约为单级时常数乘以级数。由于极点互相接近，此程序在N>4时又会出现很大误差。 程序exn914运行结果(频域)右图绘出了多级放大器的频率特性，其幅特性（图上为分贝数），显示了低通特性，随级数的增加,通带减小,从相特性看出，随级数的增加,负相移成比例地增加。 例9-1-5方波通过滤波器设方波信号的宽度为5秒，信号持续期为10秒，试求其在0～20【1/秒】频段间的频谱特性。如只取从0～10【1/秒】的频谱分量作反变换（相当于通过了一个低通滤波器），求其输出波形。解：◆建模设信号的时域波形f(t),在0到10秒的区间外信号为零,则其付利叶变换为:按MATLAB作数值计算的要求,必须把t分成N份,用相加来代替积分,对于任一ω,可写成: 积分转化为求和运算这说明求和的问题可以用f(t)行向量乘以ejωt列向量来实现.此处的Δt是t的增量,在程序中,将用dt来代替.由于要求出一系列不同的ω处的F值,而都可用同一公式,这就可以利用MATLAB中的元素群运算能力,把ω设成一个行数组,分别代入本公式左右端的ω中去,写成:F=f*exp(-j*t’*w)*Δt其中,F是与ω同长的行向量,exp中的t’是列向量,w是行向量,t’*w是一个矩阵,其行数与t相同,列数与w相同.这个式子就完成了付利叶变换,类比地可以得出付利叶逆变换表示式.由此得到下面的付利叶变换程序。 程序exn915clear,tf=10;N=256;t=linspace(0,tf,N);%给出时间分割w1=linspace(eps,20,N);dw=20/(N-1);%dw=1/4/tf;w1=[eps:dw:(N-1)/4/tf];%给出频率分割f=[ones(1,N/2),zeros(1,N/2)];%给出信号（此处是方波）F1=f*exp(-j*t"*w1)*tf/(N-1);%求付利叶变换w=[-fliplr(w1),w1(2:N)];%补上负频率F=[fliplr(F1),F1(2:N)];%补上负频率区的频谱w2=w(N/2:3*N/2);%取出中段频率F2=F(N/2:3*N/2);%取出中段频谱subplot(1,2,1),plot(w,abs(F),"linewidth",1.5),gridf1=F2*exp(j*w2"*t)/pi*dw;%对中段频谱求付利叶逆变换subplot(1,2,2),plot(t,f,t,f1,"linewidth",1.5),grid 程序exn915运行结果◆执行这个程序的结果见下图,因为方波含有很丰富的高频分量，要充分恢复其原来波形需要很宽的频带，实践中不太可能做到。 9.2离散信号和系统信号可以粗略地分为模拟信号和数字信号.模拟信号将用x(t)表示,其中变量t可以表示任何物理量,但我们假定它代表以秒为单位的时间.离散信号用x(n)表示.其中变量n为整数并代表时间的离散时刻.因此它也称为离散时间信号.他是一个数字的序列并可用以下符号之一来表述:x(n)={x(n)}={...,x(-1),x(0),x(1),...}↑其中,向上的箭头表示在n=0处的取样. 离散信号的MATLAB表示在MATLAB中,我们可以用一个列向量来表示一个有限长度的序列.然而这样一个向量并没有包含采样位置的信息.因此,完全地表示x(n)要用x和n两个向量，例如序列x(n)={2,1,-1,5,1,4,3,7}（下面的↑箭头为第0个采样点），在MATLAB中表示为:n=[-3,-2,-1,0,1,2,3,4];x=[2,1,-1,0,1,4,3,7];当不需要采样位置信息或这个信息是多余的时候(例如该序列从n=0开始),我们可以只用x向量来表示.由于有限的内存,MATLAB无法表示无限序列. 例9-2-1基本脉冲序列程序1.单位脉冲序列起点n0,终点nf,在ns处有一单位脉冲(n0≤ns≤nf),2.单位阶跃序列起点n0,终点nf,在ns前为0,在ns后为1(n0≤ns≤nf),3.复数指数序列:4.复数指数序列: 例9-2-1的程序解:◆建模这些基本序列的表达式比较简明,编写程序也不难,对单位脉冲序列,此处提供了两种方法,其中用逻辑关系的编法比较简洁.读者可从中看到MATLAB编程的灵活性和技巧性,要多用多想才能编出高明简洁的程序。绘制脉冲序列,通常用stem语句.◆MATLAB程序clear,no=0;nf=10;ns=3;n1=n0:nf;x1=[zeros(1,ns-n0),1,zeros(1,nf-ns)];%n1=n0:nf;x1=[(n1-ns)==0];%显然,用逻辑式是比较高明的方法n2=n0:nf;x2=[zeros(1,ns-n0),ones(1,nf-ns+1)];%也有类似的用逻辑比较语句的方法,留给读者思考 例9-2-1的程序(续)n3=n0:nf;x3=(0.9).^n3;%实数指数序列n4=n0:nf;x4=exp((-0.2+0.3j)*n3);%复数指数序列subplot(2,2,1),stem(n1,x1);subplot(2,2,2),stem(n2,x2);subplot(2,2,3),stem(n3,x3);subplot(4,2,6),stem(n4,real(x4));%注意sunplot的输入变元subplot(4,2,8),stem(n4,imag(x4));line([0,10],[0,0]),%画横坐标 例9-2-1的程序运行结果 例9-2-2离散付利叶变换的计算解：◆建模一个时间序列x(n)的离散时间付利叶变换的定义为：如果序列的长度是有限的，可以把它看作是周期性无限序列中的一个周期，其长度为N,对这个周期性序列可以用离散付利叶变换（注意少了‘时间’两字）进行研究。它的定义为：其中 离散傅里叶变换化为矩阵乘积用例9-1-5中的方法,引入矩阵乘法来实现求和运算,用元素群算法来求不同k时的X,把n和k都设成1×N的行数组，则nk=n’*k(以及WN.^nk)成为一个N×N的方阵则有X=x*WN.^nk即MATLAB只能处理有限长度的序列，因此，适合于计算离散付利叶变换及其逆变换。 离散傅里叶变换程序exn922设有限信号序列xn(n)的长度为Nx，则按定义，其N点付利叶变换Xk(k)的的程序为：xn=input(‘x=‘);N=(length(xn));%取N为Nxn=[0:1:N-1];k=[0:1:N-1];%设定n和k行向量WN=exp(-j*2*pi/N);%Wn因子nk=n"*k;%产生一个含nk值的N乘N维矩阵WNnk=WN.^nk;%DFT矩阵Xk=xn*WNnk;%DFT系数的行向量plot(abs(Xk)),grid%绘幅频特性图 MATLAB中的fft子程序这个程序的运算速度比较低，实际上MATLAB已提供了快速离散付利叶变换的函数fft可直接调用。其调用格式为：X=fft(x,N)x是输入的时间序列，N则是付利叶变换取的点数，若省略N，则它自动把x的长度作为N。当N取2的幂时，变换速度最快，所以要提高fft函数的运行速度，程序应编写如下：xn=input(‘x=‘);N=pow2(nextpow2((length(xn))));%取N为大于Nx而最接近的2的幂tic,X=fft(xn,N);toc要注意X是一个长度为N的复数数组，可以分解出它的振幅和相位频率分别绘图。 程序exn922运行结果令Nx=700;n=1:Nx;x=sin(0.1*n)+randn(1,Nx);然后调用上述程序之一,输入此x,所得其fft的幅度特性如右图所示。在第一个程序的末两行之前加上tic，在其后加上toc,则屏幕上会运行这两语句句所需的时间，在作者的计算机上约为8秒；而同样的fft只需0.05秒。计算逆离散付利叶变换的程序与此相仿，留待读者自行编写。信号的fft的振幅频率特性 9.3系统函数的计算机求法例9-3-1简单信号流图模型的矩阵解法[11]信号流图是用来表示和分析复杂系统内的信号变换关系的工具。其基本概念如下：（1）系统中每个信号用图上的一个节点表示。如图中的u,x1,x2。（一般物流图中是把物流标在箭杆上的。）（2）系统部件对信号实施的变换关系用有向线段表示，箭尾为输入信号，箭头为输出信号，箭身标注对此信号进行变换的乘子。如图上的G1，G2。如果乘子为1，可以不必标注。（3）每个节点信号的值等于所有指向此节点的箭头信号之和，每个节点信号可以向外输出给多个部件，其值不变。 简单信号流图的数学模型根据这几个概念，可以列出右图的方程如下。写成矩阵方程或x=QxPu带反馈的简单信号流图 矩阵方程的化简(1)移项整理，可以得到信号向量x的公式(I–Q)x=Pux=inv(I–Q)*Pu定义系统的传递函数W为输出信号与输入信号之比x/u，则W可按下式求得：W=x/u=inv(I–Q)*P 矩阵方程的化简(2)因为求二阶矩阵的逆可以直接用下面的公式：所以即对x1的传递函数为，对x2的传递函数为。 用MATLAB解信号流图对于阶次高的情况，求逆就必须用软件工具了。如果信号流图中有G1那样的符号变量，那么它的求解要用符号运算工具箱，对于本题，其MATLAB程序exn931为：symsG1G2Q=[0,-G2;G1,0],P=[1;0]W=inv(eye(2)Q)*P程序运行的结果是与前面的结果相同。用经典的“梅森公式”求信号流图的解非常繁琐，而且无法机械化。用矩阵代数方法的最大好处是可以向任意高的阶次、任意复杂的信号流图推广，实现复杂系统传递函数推导的自动化。到现在为止，我们还没有见过任何一本书籍用矩阵方法来推导这个公式。 例9-3-2较复杂的信号流图图9-11就是一个较复杂些的信号流图。照上述方法列出它的方程如下：x1=G4x3ux2=G1x1G5x4x3=G2x2x4=G3x3列为矩阵方程，得到带双重反馈的信号流图 程序exn932公式W=x/u=inv(I–Q)*P同样是正确的，不过这里的Q和P分别为44和41矩阵，用手工求逆是办不到了。我们可采用类似的MATLAB程序exn931：symsG1G2G3G4G5Q=[0,0,-G4,0;G1,0,0,-G5;0,G2,0,0;0,0,G3,0],P=[1;0;0;0]W=inv(eye(4)-Q)*PPretty(W(4))程序运行的结果为： exn932的运行结果当系统内各个环节都是线性集总参数，因而它们的传递函数Gi都可以表为s的多项式有理分式时，不管其阶次有多高，传递函数W可以很容易由计算机直接自动算出。这个方法还可以推广到离散系统，用来计算任意复杂的数字滤波器系统函数。 9.4频谱及其几何意义频谱分析是信号与系统课程中最重要的内容之一，许多读者在学习中感到抽象，往往只能从数学上承认时域信号与它的频谱之间的变换关系，而没有理解它的物理意义。用MATLAB可以帮助读者建立形象的几何概念，真正掌握它。首先来看欧拉公式，它是以最简明的方式建立了信号频域与时域的关系：它说明一个最简单的实余弦信号可以由正、负两个Ω0频率分量合成。在复平面上，正的Ω0对应于反时针旋转的向量，负的Ω0对应于顺时针旋转的向量，当这两个向量的幅度相同，而相角符号相反时，就合成为一个在实轴上的向量。它的相角为零，大小按正弦变化，形成了实信号cosΩ0t。 实周期信号是频谱向量的合成推而广之，任何实周期信号必然具有正、负两组频频率的频谱成分，正、负频率频谱的幅度对称而相位反对称，或者说，是共軛的。如果频谱不止这两项，而是有四项或更多，它们的合成仍然可以用几何动画来表示。可以把每个频谱看作一根长度等于频谱幅度、按频率Ω旋转的杆件，频谱的相加等价于多节杆件首尾相接，杆件末端的轨迹就描述了生成的时域波形。因为这个端点是在平面上运动，所以它将产生复信号，在实轴和虚轴上的投影分别为实信号和虚信号。 多个频率向量的合成程序【例9-4-1】设计一个演示程序，它能把四个用户任意给定集总频谱合成并生成对应的时域信号。解：建模按上述多节杆合成模型程序设计包括三个主要部分：（1）各频谱分量的输入，包括其幅度和频率（有正负号）；（2）将各分量当作转动的杆件首尾相接；（3）记录多节杆系末端的轨迹画出图形。 程序exn941%（1）给频谱向量赋值N=input("N(输入向量个数，限定N不大于4)=");fori=1:Ni,a(i)=input("振幅a(i)=");w(i)=input("角频率w(i)=");end%（2）将各个频谱向量相加合成并画图%此处应该把各时刻的图形转为动画，语句省略）t=0:0.1:20;lt=length(t);%给出时间数组%各频谱分量随时间变化的复数值p=a"*ones(1,lt).*exp(j*w"*t);q=cumsum(p);%各频谱分量的累加 exn941(绘图语句)%画合成复信号的端点轨迹figure(1),plot(real(q(4,:)),imag(q(4,:))),gridon%（3）将此轨迹的在实轴与虚轴两个方向的投影画成时间信号figure(2),subplot(2,2,1),plot(real(q(4,:)),imag(q(4,:))),gridon%画出实信号的时间波形subplot(2,2,3),plot(real(q(4,:)),t),gridon%画出虚信号的时间波形subplot(2,2,2),plot(t,imag(q(4,:))),gridon 程序exn941运行的结果运行程序并按提示，输入四个频谱向量组成的多节杆及其端点轨迹a(1)=1，w(1)=-1；a(2)=1，w(2)=--1；a(3)=0.5，w(3)=3；a(4)=0.5，w(1)=-4；在此处，为了显示复信号，我们有意把输入频谱设成不对称的。于是读者将看到四节杆的运动动画，并得到杆系及其端点在复平面上的轨迹，如图所示。四个频谱向量组成的多节杆及其端点轨迹 不对称频谱生成的实、虚信号改变了比例尺的轨迹见图9-13子图(a)。将它在x,y两方向的投影与时间轴的关系画在子图(b)和(c)中，我们就得到信号与系统课程中常见的实信号曲线。图9-13四个集总频谱生成的复信号及其实信号分量。四个集总频谱生成的复信号及其实信号分量 进一步理解负频率的意义频谱的频率则只能是可正可负的实数，正频率和负频率以及在该频率上频谱的意义在此不言自明。读者可以做各种各样的试验，例如当两组频率具有倍频关系时，得到的是周期信号，如果频率比是任意小数，那将得出非周期的信号；可见，只有在复信号平面上，才能看到频率的正负。许多人往往看不到负频率频谱的意义，这是因为他们老是停留在实信号的局部范畴来思考问题。要想知道象有没有鼻子（负频率），必须去摸象头（复信号），只摸象腿（实信号）是得不到信号理论普遍而科学的规则的。'