第五讲信号系统分析的MATLAB实现

第五讲信号系统分析的MATLAB实现基于MATLAB的离散时间系统分析基于MATLAB的连续时间系统分析一、离散时间信号的表示表示离散序列并将其可视化时：用向量表示，无法用符号运算表示无法表示无限序列绘制图形时，用stem指令X(n)={1,2,-1,3,2,4,-1}0kk=-3:3;f=[1,2,-1,3,2,4,-1];stem(k,f),axis([-4,4,-1.5,4.5])单位脉冲序列)(n)0(nn=1n+n0=0=0其它functiondwxulie(k1,k2,k0)k=k1:k2;n=length(k);f=zeros(1,n);f(1,-k0-k1+1)=1;stem(k,f)axis([k1,k2,0,1.5])title(‘单位序列’）dwxulie(-5,5,0)单位阶跃序列)0(kkfunctionjyxulie(k1,k2,k0)k=k1:-k0-1kk=-k0:k2;n=length(k);nn=length(kk)u=zeros(1,n);uu=ones(1,nn)stem(kk,uu)holdonstem(k,u)holdofftitle(‘单位阶跃序列’）axis([k1,k2,0,1.5])jyxulie(-3,8,0)二、序列的运算1、信号加数学描述：MATLAB描述：x=x1+x2注意：x1和x2应该具有相同的长度，位置对应，才能相加，否则，需要先通过函数补零后再相加)(2)(1)(nxnxnxfunction[f,k]=lsxj(f1,f2,k1,k2)k=min(min(k1),min(k2)):max(max(k1),max(k2));s1=zeros(1,length(k));s2=s1;s1(find((k=min(k1))&(k=max(k1))==1))=f1;s2(find((k=min(k2))&(k=max(k2))==1))=f2;f=s1+s2;stem(k,f,’filled’)axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)])例：有两离散序列，用MATLAB绘出它们的波形及f1(k)+f2(k).解：f1=-2:2;k1=-2:2;f2=[1,1,1];k2=-1:1;stem(k1,f1),axis([-3,3,-2.5,2.5])f1(k)={-2,-1,0,1,2}f2(k)={1,1,1}stem(k2,f2),axis([-3,3,-2.5,2.5])[f,k]=lsxj(f1,f2,k1,k2)f=-20122k=-2-10122、离散序列的卷积MATLAB实现01)(1kf0=k=2其它}0,3,2,1{)(2kfk=1其它f1=ones(1,3);f2=0:3;f=conv(f1,f2)运行结果：f=013653function[f,k]=deconv(f1,f2,k1,k2)f=conv(f1,f2)k0=k1(1)+k2(1);k3=length(f1)+length(f2)-2;k=k0:k0+k3subplot(2,2,1)stem(k1,f1)title('f1(k)')subplot(2,2,2)stem(k2,f2)title('f2(k)')suplot(2,2,3)stem(k,f)f1=[1,1,1];k1=0:2;f2=[1,2,3];k2=1:3;[f,k]=deconv(f1,f2,k1,k2)f=13653k=12345三、差分方程filtery=filter(b,a,x)a(1)*y(n)=b(1)*x(n)+b(2)*x(n-1)+...+b(nb+1)*x(n-nb)-a(2)*y(n-1)-...-a(na+1)*y(n-na)例：y(n)-y(n-1)+0.9y(n-2)=x(n)a.计算并画出n=-20，……,100的脉冲响应b.计算并画出n=-20…,100的单位阶跃响应c.这个系统是稳定的吗？b=[1];a=[1,-1,0.9];s=dwxulie(-20,120,0);n=[-20:120];s=dwxulie(-20,120,0);n=[-20:120];h=filter(b,a,s);subplot(2,1,1);stem(n,h);title('impulseresponse');xlabel('n');ylabel('h(n)')单位阶跃响应x=jyxulie(-20,120,0);s=filter(b,a,x);subplot(2,1,2);stem(n,s)title('stepresponse');xlabel('n');ylabel('s(n)')四、Z变换Z反变换residuez[R,P,C]=residuez(b,a)nmkkknkKKnnmmzCZPRZAZBzazaazbzbbzx0111101101)()(......)(例：求211243143)(zzzzzzzX的z反变换b=[0,1];a=[3,-4,1];[R,P,C]=residuez(b,a)R=0.5000-0.5000P=1.00000.3333C=[]113333.015.015.0)(zzzX五、从差分方程表示系统函数roots－可以对分子分母求零点、极点zplane(b,a)－在已知分子行向量b和分母行向量a下画出零点和极点[H,W]=freqz(b,a,N)－b是由分子多项式系数构成的数组，a是由分母多项式系数构成的数组，h是返回的复数形式的频率响应数组，w是由对应于h的各个数字角频率构成的数组。1z1z例：给定一因果系统)()1(9.0)(nxnynya、求H(Z)并大致画出它的零极点图。b、画出和c、求脉冲响应h(n)解：)(jeH)(jeH19.011)(ZZHb=[1,0];a=[1,-0.9];zplane(b,a)[H,W]=freqz(b,a,100);magh=abs(H);phah=angle(H);subplot(2,1,1);plot(W/pi,magh);gridsubplot(2,1,2);plot(W/pi,phah/pi);grid六、Fourier分析abs－模angle－相角freqz－知道z变换，求频率特性。freqs－知道s变换，求频率特性。fft－快速离散Fourier变换fft2－二维快速离散Fourier变换fftn－高维快速离散Fourier变换ifft－快速离散Fourier逆变换ifft2－二维快速离散Fourier逆变换abs功能：求绝对值（幅值）调用格式：y=abs(x)angle功能：求相角调用格式：p=angle(h)freqs功能：模拟滤波器的频率响应调用格式h=freqs(b,a,w)[h,w]=freqs(b,a)[h,w]=freqs(b,a,n))1(...)2()1()1(...)2()1()()()(11naasasanbbsbsbsAsBsHnananbnbfreqspace功能：设置频率响应中的频率间隔调用格式：f=freqspace(n)f=freqspace(n,’whole’)freqzfreqzplot功能：利用频率响应数据绘图调用格式：freqzplot(h,w)freqz(h,w,s)impz功能：计算数字滤波器的冲激响应调用格式：[h,t]=impz(b,a)[h,t]=impz(b,a,n)impz(b,a)例：一个四阶的低通椭圆滤波器，绘制其冲激响应的前50个采样点。[b,a]=ellip(4,0.5,20,0.4);impz(b,a,50);gridonSignalProcessingToolboxFilterDesign&AnalysisTool(FDATool)SignalProcessingTool(SPTool)