MATLAB在信号与系统中的应用

实验4MATLAB在信号与系统中的应用一、连续信号和系统【例6.1】连续信号的MATLAB描述·程序cleart0=0;tf=5;dt=0.05;t1=1;t=[t0:dt:tf];st=length(t);%（1）单位冲激函数n1=floor((t1-t0)/dt);x1=zeros(1,st);x1(n1)=1/dt;subplot(2,2,1),stairs(t,x1),gridonaxis([0,5,0,1.1])%（2）单位阶跃函数x2=[zeros(1,n1-1),ones(1,st-n1+1)];subplot(2,2,3),stairs(t,x2),gridonaxis([0,5,0,1.1])%（3）复指数函数alpha=-0.5;w=10;x3=exp((alpha+j*w)*t);subplot(2,2,2),plot(t,real(x3)),gridonsubplot(2,2,4),plot(t,imag(x3)),gridon·总结1、单位冲激函数无法直接描述，可以把其看成宽度是dt的，幅度是1/dt的矩形脉冲。2、用stairs命令，可以显示连续信号波形中的不连续点；用plot命令，使波形更光滑。3、复指数信号可以分解为余弦和正弦信号，它们分别是复数信号的实部和虚部。【例6.2】LTI系统的零输入响应a=input('输入分母系数向量a=[a1,a2,...]=');n=length(a)-1;Y0=input('输入初始条件向量Y0=[y0,Dy0,D2y0,...]=');p=roots(a);V=rot90(vander(p));c=V\Y0';dt=input('dt=');tf=input('tf=');t=0:dt:tf;y=zeros(1,length(t));fork=1:ny=y+c(k)*exp(p(k)*t);endplot(t,y),gridon·总结1、当LTI系统的输入为0时，其零输入响应为微分方程的齐次解，即令微分方程等号右端为0的解。2、特征方程的根，可以用roots(a)语句求得。【例6.3】n阶LTI系统的冲激响应a=input('多项式分母系数向量a=')多项式分母系数向量a=poly([0,-1+2i,-1-2i,-2,-5]);a=192955500b=input('多项式分子系数向量b=')多项式分子系数向量b=[8,3,1];b=831[r,p]=residue(b,a),r=0.62000.1300-0.3900i0.1300+0.3900i-0.90000.0200p=-5.0000-1.0000+2.0000i-1.0000-2.0000i-2.00000disp('解析式h(t)=r(i)*exp(p(i)*t)')解析式h(t)=r(i)*exp(p(i)*t)disp('给出时间数组t=[0:dt:tf]')给出时间数组t=[0:dt:tf]dt=input('dt=');dt=0.2tf=input('tf=');tf=8t=0:dt:tf;h=zeros(1,length(t));fori=1:length(a)-1h=h+r(i)*exp(p(i)*t);endplot(t,h),grid调用工具箱函数：a=input('多项式分母系数向量a=(书上取poly([0,-1+2j,-1-2j,-2,-5])）');b=input('多项式分子系数向量b=(书上取[8,3,1]）');[r,p]=residue(b,a),%求留数disp('解析式h(t)=Σr(i)*exp(p(i)*t)')disp('给出时间数组t=[0:dt:tf]')dt=input('dt=');%给定时间数组tf=input('tf=');t=0:dt:tf;y=impulse(b,a,t)plot(t,y),grid·总结1、冲激函数的拉式变换等于1，则系统对冲激函数的响应的拉式变换Y(s)=H(s)U(s)=H(s)。2、冲激响应就是H(s)的拉式反变换。【例6.4】卷积的计算clear,closeallu1s=input('输入u数祖u=（例如ones(1,10)）');lu=length(u1s);h1s=input('输入h数祖h=（例如exp(-0.1*[1:15])）');lh=length(h1s);lmax=max(lu,lh);iflulhnu=0;nh=lu-lh;%若u比h长，对h补nh个零elseiflulhnh=0;nu=lh-lu;%若h比u长，对u补nu个零elsenu=0;nh=0;enddt=input('输入时间间隔dt=（例如0.5）');lt=lmax;u=[zeros(1,lt),u1s,zeros(1,nu),zeros(1,lt)];t1=(-lt+1:2*lt)*dt;h=[zeros(1,2*lt),h1s,zeros(1,nh)];hf=fliplr(h);y=zeros(1,3*lt);fork=0:2*ltp=[zeros(1,k),hf(1:end-k)];y1=u.*p*dt;yk=sum(y1);y(k+lt+1)=yk;%绘图，注意如何用axis命令把各子图的横坐标统一起来，使纵坐标随数据自动调整subplot(4,1,1);stairs(t1,u)%用stairs是为了避免plot函数在突跳点形成的斜边axis([-lt*dt,2*lt*dt,min(u),max(u)]),holdonylabel('u(t)')subplot(4,1,2);stairs(t1,p)axis([-lt*dt,2*lt*dt,min(p),max(p)])ylabel('h(k-t)')subplot(4,1,3);stairs(t1,y1)axis([-lt*dt,2*lt*dt,0,1])ylabel('s=u.*h(k-t)')subplot(4,1,4);stem(k*dt,yk,'.')%用stem函数表示每一次卷积积分求和的结果axis([-lt*dt,2*lt*dt,floor(min(y)+eps),ceil(max(y+eps))]),holdonylabel('y(k)=sum(s)*dt')ifk==round(0.8*lt)disp('暂停，按任意键继续'),pauseelsepause(0.4)endend二、傅里叶分析【例6.7】方波分解为多次正弦波之和t=0:.01:2*pi;y=sin(t);plot(t,y),figure(gcf),pausey=sin(t)+sin(3*t)/3;plot(t,y),pause%用1，3，5，7，9次谐波叠加y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;plot(t,y)%为了绘制三维曲面，要把各次波形数据存为一个三维数组，因此必须重新定义y，重编程。y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:19x=x+sin(k*t)/k;y((k+1)/2,:)=x;end%将各波形迭合绘出pause,figure(1),plot(t,y(1:9,:)),gridline([0,pi+0.5],[pi/4,pi/4])text(pi+0.5,pi/4,'pi/4')%将各半波形绘成三维网格图，看出增加谐波阶次对方波逼近程度的影响halft=ceil(length(t)/2);pause,figure(2),mesh(t(1:halft),[1:10],y(:,1:halft)),pause【例6.9】周期信号的滤波clear,formatcompactL=0.4;C=10e-6;R=200;Um=100;w1=100*pi;%输入只有偶次谐波，故只分析偶次谐N=input('需分析的谐波次数2N=（键入偶数）');n=1:N/2;w=[eps,2*n*w1];Us=4*Um/pi*[0.5,-1./(4*n.^2-1)];z1=j*w*L;z2=1./(j*w*C);z3=R;z23=z2.*z3./(z2+z3);Ur=Us.*z23./(z1+z23)disp('谐波次数谐波幅度谐波相移(度)')disp([2*[0,n]',abs(Ur)',angle(Ur)'*180/pi])需分析的谐波次数2N=（键入偶数）10Ur=63.6620-0.0000i12.8383+27.8571i1.3050+0.6169i0.2546+0.0726i0.0799+0.0165i0.0326+0.0053i谐波次数谐波幅度谐波相移(度)063.6620-0.00002.000030.673165.25684.00001.443425.30146.00000.264815.92538.00000.081611.702910.00000.03309.2740【例6.10】调幅信号通过带通滤波器clearN=1000;t=linspace(0,2*pi,N);dt=2*pi/(N+1);%信号周期为2*pi,分成1000份w=[99,100,101];%输入信号的三个频率分量U=[0.5,1,0.5];%三个频率分量对应的向量（虚部为零）b=[2,0];a=[1,2,10001];%滤波器分子分母系数向量u1=U*cos(w'*t+angle(U')*ones(1,N));%输入信号的时间曲线H=polyval(b,j*w)./polyval(a,j*w);%求滤波器在三个频点上的频率响应，也可用H=freqs(b,a,w);%画出滤波器的频率响应曲线，只用三个频点，图形不好看，读者可修改程序得到其完整的频率响应figure(1)subplot(2,1,1),plot(w,abs(H)),grid%幅度subplot(2,1,2),plot(w,angle(H)),grid%相位%U2=U.*H;%ewn=exp(w'*t*i);u2=U2*ewn;%求u2=idtft(U2,w,t);u21=abs(U(1)*H(1))*cos(99*t+angle(U(1)*H(1)));%角频率为99的分量u22=abs(U(2)*H(2))*cos(100*t+angle(U(2)*H(2)));%角频率为100的分量u23=abs(U(3)*H(3))*cos(101*t+angle(U(3)*H(3)));%角频率为101的分量u2=u21+u22+u23;%求和%巧妙地利用元素群运算和矩阵运算相结合可把四条语句合成一条语句如下%u2=abs(U.*H)*cos(w'*t+angle((U.*H).')*ones(1,1001));%注意对复数矩阵(U.*H)，(U.*H)'为其共轭转置，(U.*H).'为转置而不共轭figure(2)%画出原信号和滤波后信号的波形作比较subplot(2,1,1),plot(t,u1)subplot(2,1,2),plot(t,u2)a)clearN=2^16;tf=2*pi;t=linspace(0,tf,N);dt=tf/N;%信号周期为2*pi,分成1000份w1=[99,100,101];%输入信号的三个频率分量U=[0.5,1,0.5];%三个频率分量对应的向量（虚部为零）b=[2,0];a=[1,2,10001];%滤波器分子分母系数向量u1=U*cos(w1'*t+angle(U')*ones(1,N));%输入信号的时间曲线U1=fft([u1,u1],2*N)*dt;w=[0:2*N-1]*2*pi/2/N/dt;H=freqs(b,a,w);%求滤波器在三个频点上的频率响应，也可用H=freqs(b,a,w);%画出滤波器的频率响应曲线，只用三个频点，图形不好看，读者可修改程