九讲MATLAB在电路中的应用课件

利用利用MATLABMATLAB进行进行信号与系统分析信号与系统分析 信号的信号的MATLABMATLAB表示表示 基本信号的基本信号的MATLABMATLAB表示表示 指数信号指数信号指数信号指数信号AeAeAeAeatatatat 、指数序列、指数序列、指数序列、指数序列a a a ak k k k 、抽样函数、抽样函数、抽样函数、抽样函数Sa(tSa(tSa(tSa(t)、正弦型信号、矩形脉冲信号、三角、正弦型信号、矩形脉冲信号、三角、正弦型信号、矩形脉冲信号、三角、正弦型信号、矩形脉冲信号、三角脉冲信号脉冲信号脉冲信号脉冲信号 信号基本运算的信号基本运算的MATLABMATLAB实现实现尺度变换、翻转、时移、尺度变换、翻转、时移、尺度变换、翻转、时移、尺度变换、翻转、时移、相加、相乘、差分相加、相乘、差分相加、相乘、差分相加、相乘、差分与求和、微分与积分与求和、微分与积分与求和、微分与积分与求和、微分与积分 一、基本信号的一、基本信号的MATLABMATLAB表示表示 指数信号指数信号AeAeataty=A*y=A*expexp(a(a*t);*t);指数序列指数序列akak幂运算幂运算a.ka.k实现实现正弦型信号正弦型信号内部函数内部函数coscos()()和和sin()sin()抽样函数抽样函数Sa(tSa(t)sinc(tsinc(t)矩形脉冲信号矩形脉冲信号 y=y=rectpuls(t,widthrectpuls(t,width)三角波脉冲信号三角波脉冲信号y=y=tripuls(ttripuls(t,width,skewwidth,skew)%decaying exponentialt=0:001:10;A=1;a=-0.4;ft=A*exp(a*t);plot(t,ft)t=0:0.1:10;A=1;a=-0.4;ft=A*exp(a*t);stem(t,ft)一、基本信号的一、基本信号的MATLABMATLAB表示表示%rectpulst=0:0.001:4;T=1;ft=rectpuls(t-2*T,T);plot(t,ft)axis(0,4,-0.5,1.5)一、基本信号的一、基本信号的MATLABMATLAB表示表示%tripulst=-3:0.001:3;ft=tripuls(t,4,0.5);plot(t,ft)ft=tripuls(t,4,1);一、基本信号的一、基本信号的MATLABMATLAB表示表示%unit%unit impulsimpuls sequence sequence k=-50:50;k=-50:50;delta=zeros(1,50),1,zeros(1,50);delta=zeros(1,50),1,zeros(1,50);stem(k,deltastem(k,delta)function function f,kf,k=impseqimpseq(k0,k1,k2)(k0,k1,k2)%产生产生产生产生 fkfk=delta(k=delta(k-k0)k0)；k1=k=k2 k1=k=k2 k=k1:k2;f=(kk=k1:k2;f=(k-k0)=0;k0)=0;k0=0;k1=k0=0;k1=-50;k2=50;50;k2=50;f,kf,k=impseqimpseq(k0,k1,k2);(k0,k1,k2);stem(k,fstem(k,f)一、基本信号的一、基本信号的MATLABMATLAB表示表示%unit step sequence%unit step sequencek=k=-50:50;50:50;ukuk=zeros(1,50),ones(1,51);=zeros(1,50),ones(1,51);stem(k,ukstem(k,uk)function function f,kf,k=stepseqstepseq(k0,k1,k2)(k0,k1,k2)%产生产生产生产生 fkfk=u(k=u(k-k0);k1=k=k2 k0);k1=k=0;k0)=0;k0=0;k1=k0=0;k1=-50;k2=50;50;k2=50;f,kf,k=stepseqstepseq(k0,k1,k2);(k0,k1,k2);stem(k,fstem(k,f)function function f,kf,k=stepseqstepseq(k0,k1,k2)(k0,k1,k2)%产生产生产生产生 fkfk=u(k=u(k-k0);k1=k=k2 k0);k1=k=0;k0)=0;k0=0;k1=k0=0;k1=-50;k2=50;50;k2=50;f,kf,k=stepseqstepseq(k0,k1,k2);(k0,k1,k2);stem(k,fstem(k,f)一、基本信号的一、基本信号的MATLABMATLAB表示表示 二、信号基本运算的二、信号基本运算的MATLABMATLAB实现实现t=t=-3:0.001:3;3:0.001:3;ft1=tripuls(ft1=tripuls(2*t2*t,4,0.5);,4,0.5);subplot(2,1,1)subplot(2,1,1)plot(t,ft1)plot(t,ft1)title(title(f(2t)f(2t)ft2=tripuls(ft2=tripuls(2 2-2*t2*t),4,0.5);),4,0.5);subplot(2,1,2)subplot(2,1,2)plot(t,ft2)plot(t,ft2)title(title(f(2f(2-2t)2t)1.信号的尺度信号的尺度变换、翻、翻转、时移（平移）移（平移）已知三角波已知三角波已知三角波已知三角波f f(t t)，用，用，用，用MATLABMATLAB画出的画出的画出的画出的f f(2(2t t)和和和和f f(2(2-2 2t t)波形波形波形波形二、信号基本运算的二、信号基本运算的MATLABMATLAB实现实现2.信号的相加与相乘信号的相加与相乘相加相加相加相加用算术运算符用算术运算符用算术运算符用算术运算符“+”实现实现实现实现 相乘相乘相乘相乘用数组运算符用数组运算符用数组运算符用数组运算符“.*.*”实现实现实现实现 oo 例例例例:画信号画信号画信号画信号AeAeatatcos(cos(0 0t t+f f f f)的波形的波形的波形的波形 t=0:0.001:8;t=0:0.001:8;A=1;a=A=1;a=-0.4;0.4;w0=2*w0=2*pi;phipi;phi=0;=0;ft1=A*ft1=A*exp(aexp(a*t)*t).*.*sin(w0*sin(w0*t+phit+phi););plot(t,ft1)plot(t,ft1)3.离散序列的差分与求和离散序列的差分与求和 连续信号的微分与信号的微分与积分分差分差分差分差分y=y=diffdiff(f(f););求和求和求和求和y=y=sumsum(f(k1:k2);(f(k1:k2);微分微分微分微分 y=y=diff(f)/hdiff(f)/h;h;h为数值计算所取时间间隔为数值计算所取时间间隔为数值计算所取时间间隔为数值计算所取时间间隔 定积分定积分定积分定积分 quadquad(function_name,a,b(function_name,a,b););function_namefunction_name为被积函数名，为被积函数名，为被积函数名，为被积函数名，a a和和和和b b指定积分区间。指定积分区间。指定积分区间。指定积分区间。二、信号基本运算的二、信号基本运算的MATLAB实现实现3.离散序列的差分与求和离散序列的差分与求和 连续信号的微分与信号的微分与积分分 例例例例:已知三角波已知三角波f(t)，画出其微分与积分的波形，画出其微分与积分的波形%differentiation%differentiation h=0.001;t=h=0.001;t=-3:h:3;3:h:3;y1=diff(ft2)*1/h;y1=diff(ft2)*1/h;plot(t(1:length(t)plot(t(1:length(t)-1),y1)1),y1)%integration%integration t=t=-3:0.1:3;3:0.1:3;forfor x=1:length(t)x=1:length(t)y2(x)=quad(y2(x)=quad(ft2ft2,-3,t(x);3,t(x);endend plot(t,y2)plot(t,y2)三角波三角波f(t)微分与积分的波形微分与积分的波形利用利用MATLAB进行系统的时域分析进行系统的时域分析连续时间系统零状态响应的求解连续时间系统零状态响应的求解 连续时间系统冲激响应和阶跃响应的求解连续时间系统冲激响应和阶跃响应的求解 离散时间系统零状态响应的求解离散时间系统零状态响应的求解 离散时间系统单位脉冲响应的求解离散时间系统单位脉冲响应的求解 离散卷积的计算离散卷积的计算一、连续时间系统零状态响应的求解一、连续时间系统零状态响应的求解t 表示计算系统响应的抽样点向量表示计算系统响应的抽样点向量f是系统输入信号向量是系统输入信号向量sys是是LTI系统模型，借助系统模型，借助tf函数获得函数获得a=a3,a2,a1,a0;a=a3,a2,a1,a0;b=b3,b2,b1,b0;b=b3,b2,b1,b0;sys=sys=tf(b,atf(b,a)y=y=lsimlsim(sys,f,t(sys,f,t)sys=sys=tf tf(b,a(b,a)b和和a分别为微分方程右端和左端各项的系分别为微分方程右端和左端各项的系二、连续时间系统冲激响应和阶跃响应的求解二、连续时间系统冲激响应和阶跃响应的求解 连续时间系统冲激响应可用连续时间系统冲激响应可用连续时间系统冲激响应可用连续时间系统冲激响应可用impulseimpulse函数直接求出，函数直接求出，函数直接求出，函数直接求出，其调用形式为其调用形式为其调用形式为其调用形式为y=y=impulseimpulse(sys,t(sys,t)连续时间系统阶跃响应可用连续时间系统阶跃响应可用连续时间系统阶跃响应可用连续时间系统阶跃响应可用stepstep函数直接求出，函数直接求出，函数直接求出，函数直接求出，其调用形式为其调用形式为其调用形式为其调用形式为y=y=stepstep(sys,t(sys,t)t t 表示计算系统响应的抽样点向量表示计算系统响应的抽样点向量表示计算系统响应的抽样点向量表示计算系统响应的抽样点向量 sys sys 是是是是LTILTI系统模型系统模型系统模型系统模型三、离散时间系统零状态响应的求解三、离散时间系统零状态响应的求解b,a b,a 分别是差分方程左、右端的系数向量分别是差分方程左、右端的系数向量分别是差分方程左、右端的系数向量分别是差分方程左、右端的系数向量b=b0,b1,b2,b=b0,b1,b2,bM;,bM;a=a0,a1,a2,a=a0,a1,a2,aNaN;可用可用可用可用MATLABMATLAB表示为表示为表示为表示为y y=filterfilter(b,a,f(b,a,f)f f 表示输入序列表示输入序列表示输入序列表示输入序列,y,y 表示输出序列表示输出序列表示输出序列表示输出序列四、离散时间系统单位脉冲响应的求解四、离散时间系统单位脉冲响应的求解b,a b,a 分别是差分方程左、右端的系数向量分别是差分方程左、右端的系数向量分别是差分方程左、右端的系数向量分别是差分方程左、右端的系数向量 k k 表示输出序列的取值范围表示输出序列的取值范围表示输出序列的取值范围表示输出序列的取值范围 h h 就是单位脉冲响应就是单位脉冲响应就是单位脉冲响应就是单位脉冲响应h=h=impzimpz(b,a,k(b,a,k)五、离散卷积的计算五、离散卷积的计算例：例：例：例：(s(s3 3+2s+3)(s+2s+3)(s2 2+3s+2)+3s+2)可用下面可用下面可用下面可用下面MATLABMATLAB语句求出语句求出语句求出语句求出 a=1,0,2,3;a=1,0,2,3;b=1,3,2;b=1,3,2;c=c=convconv(a,b(a,b)c=c=convconv(a,b(a,b)式中式中式中式中a,ba,b为待卷积两序列的向量表示，为待卷积两序列的向量表示，为待卷积两序列的向量表示，为待卷积两序列的向量表示，c c是卷积结果。是卷积结果。是卷积结果。是卷积结果。convconv函数也可用于计算两个多项式的积函数也可用于计算两个多项式的积函数也可用于计算两个多项式的积函数也可用于计算两个多项式的积例例1 求系统求系统 y(t)+2y(t)+100y(t)=10f(t)的零状态的零状态响应，已知响应，已知f(t)=sin(2p pt)u(t)。%program3_1%program3_1微分方程求解微分方程求解微分方程求解微分方程求解 ts ts=0;te=5;dt=0.01;=0;te=5;dt=0.01;sys=sys=tf tf(1,1 2 100);(1,1 2 100);t=t=ts:dt:tets:dt:te;f=10*sin(2*pi*t);f=10*sin(2*pi*t);y=y=lsimlsim(sys,f,t(sys,f,t););plot(t,yplot(t,y););xlabel(xlabel(Time(secTime(sec)ylabel(ylabel(y(ty(t)例例2 求系统求系统 y(t)+2y (t)+100y(t)=10f(t)的零状态响应，已知的零状态响应，已知f(t)=d d(t)。%program3_2%program3_2连续时间系统的冲激响应连续时间系统的冲激响应连续时间系统的冲激响应连续时间系统的冲激响应 ts ts=0;te=5;dt=0.01;=0;te=5;dt=0.01;sys=sys=tf tf(10,1 2 100);(10,1 2 100);t=t=ts:dt:tets:dt:te;y=y=impulseimpulse(sys,t(sys,t););plot(t,yplot(t,y););xlabel(xlabel(Time(secTime(sec)ylabel(ylabel(h(th(t)例例例例3 3 分析噪声干扰的信号分析噪声干扰的信号分析噪声干扰的信号分析噪声干扰的信号f f k k=s s k k+d d k k 通过通过通过通过MM点滑动平点滑动平点滑动平点滑动平均系统的响应，均系统的响应，均系统的响应，均系统的响应，其中其中其中其中s sk k=(2=(2k k)0.9)0.9k k是原始信号，是原始信号，是原始信号，是原始信号，d d k k 是噪声。是噪声。是噪声。是噪声。R=51;d=rand(1,R)R=51;d=rand(1,R)-0.5;0.5;k=0:Rk=0:R-1;1;s=2*k.*(0.9.k);f=s=2*k.*(0.9.k);f=s+ds+d;figure(1);figure(1);plot(k,d,plot(k,d,r-.r-.,k,s,k,s,b-b-,k,f,k,f,gg-););M=5;b=ones(M,1)/M;a=1;M=5;b=ones(M,1)/M;a=1;y=y=filterfilter(b,a,f(b,a,f););figure(2);figure(2);plot(k,s,plot(k,s,b-b-,k,y,k,y,r r-););例例例例3 3 分析噪声干扰的信号分析噪声干扰的信号分析噪声干扰的信号分析噪声干扰的信号f f k k=s s k k+d d k k 通过通过通过通过MM点滑点滑点滑点滑动平均系统的响应，动平均系统的响应，动平均系统的响应，动平均系统的响应，其中其中其中其中s sk k=(2=(2k k)0.9)0.9k k是原始信号，是原始信号，是原始信号，是原始信号，d d k k 是噪声。是噪声。是噪声。是噪声。噪声干扰信号噪声干扰信号fk=sk+dk通过通过M点滑动平均系统的响应点滑动平均系统的响应例例4 求系统求系统yk+3yk-1+2yk-1=10f k的单的单位脉冲响应。位脉冲响应。%program program 3_4 3_4 离散系统的单位脉冲响应离散系统的单位脉冲响应离散系统的单位脉冲响应离散系统的单位脉冲响应 k=0:10;k=0:10;a=1 3 2;a=1 3 2;b=1;b=1;h=h=impzimpz(b,a,k(b,a,k););stem(k,hstem(k,h)根据卷积公式：根据卷积公式：因此编程的过程为：因此编程的过程为：（1）写出）写出h(t)e-0.1t的的MATLAB表达式；表达式；（2）写出）写出u(t)的的MATLAB表达式；表达式；（3）利用）利用MATLAB的卷积语句的卷积语句y=conv(u，h)求解求解（4）画曲线）画曲线plot(t,y)。例例例例5 5 计算计算计算计算x x k k*y y k k 并画出卷积结果，已知并画出卷积结果，已知并画出卷积结果，已知并画出卷积结果，已知x x k k=1,2,3,4;=1,2,3,4;k k=0,1,2,3=0,1,2,3，y y k k=1,1,1,1,1;=1,1,1,1,1;k k=0,1,2,3,4=0,1,2,3,4。%program 3_5%program 3_5 x=1,2,3,4;x=1,2,3,4;y=1,1,1,1,1;y=1,1,1,1,1;z=z=convconv(x,y(x,y););N=N=lengthlength(z(z););stem(0:Nstem(0:N-1,z);1,z);例6 方波分解为多次正弦波之和n图示的周期性方波，其傅里叶级数为n分别计算 n直到9次谐波，并做图。输入周期性方波-x利用利用MATLAB进行信号的频域分析进行信号的频域分析 周期信号周期信号频谱的的MATLAB实现 用数值积分分析非周期信号频谱用数值积分分析非周期信号频谱一、周期信号频谱的一、周期信号频谱的MATLAB实现实现 频谱频谱频谱频谱FnFn一般为复数，可分别利用一般为复数，可分别利用一般为复数，可分别利用一般为复数，可分别利用absabs和和和和angleangle函函函函数获得其幅频特性和相频特性。数获得其幅频特性和相频特性。数获得其幅频特性和相频特性。数获得其幅频特性和相频特性。其调用格式分别为其调用格式分别为其调用格式分别为其调用格式分别为x=x=abs(Fnabs(Fn)y=y=angle(Fnangle(Fn)周期信号的频谱周期信号的频谱周期信号的频谱周期信号的频谱Fn Fn 为离散信号，可以用为离散信号，可以用为离散信号，可以用为离散信号，可以用stemstem例例例例1 1 试用试用试用试用MATLABMATLAB画出图示周期三角波信号的频谱。画出图示周期三角波信号的频谱。画出图示周期三角波信号的频谱。画出图示周期三角波信号的频谱。解：周期信号的频谱为解：周期信号的频谱为画三角波信号频谱的画三角波信号频谱的MATLAB程序程序N=8;N=8;n1=n1=-N:N:-1;1;%计算计算计算计算n=n=-N N到到到到-1 1的的的的FourierFourier系数系数系数系数 c1=c1=-4*j*sin(n1*pi/2)/pi2./n1.2;4*j*sin(n1*pi/2)/pi2./n1.2;c0=0;c0=0;%计算计算计算计算n=0n=0时的时的时的时的FourierFourier系数系数系数系数 n2=1:N;n2=1:N;%计算计算计算计算n=1n=1到到到到N N的的的的FourierFourier系数系数系数系数 c2=c2=-4*j*sin(n2*pi/2)/pi2./n2.2;4*j*sin(n2*pi/2)/pi2./n2.2;cncn=c1 c0 c2;=c1 c0 c2;n=n=-N:N;N:N;subplot(2,1,1);subplot(2,1,1);stem(n,stem(n,absabs(cn);ylabel(cn);ylabel(CnCn的幅度的幅度的幅度的幅度););subplot(2,1,2);subplot(2,1,2);stem(n,stem(n,angleangle(cn(cn););ylabel(ylabel(CnCn的相的相的相的相位位位位);xlabel();xlabel(omega/omega0omega/omega0););程序运行结果程序运行结果例例例例2 2 求周期矩形脉冲的求周期矩形脉冲的求周期矩形脉冲的求周期矩形脉冲的FourierFourier级数表示式。并用级数表示式。并用级数表示式。并用级数表示式。并用MATLABMATLAB求出由前求出由前求出由前求出由前N N项项项项FourierFourier级数系数得出的信号级数系数得出的信号级数系数得出的信号级数系数得出的信号近似波形。近似波形。近似波形。近似波形。取取取取A A=1,=1,T T=2,=2,t t t t=1,=1,0 0=p p p p%Gibbs phenomenont=t=-2:0.001:2;2:0.001:2;N=N=input(input(NumberNumber of harmonics=of harmonics=););c0=0.5;c0=0.5;cNcN=c0*ones(1,length(t);=c0*ones(1,length(t);%dc component%dc component for for n=0:2:Nn=0:2:N-1 1%even harmonics are zero%even harmonics are zero cNcN=cN+cos(picN+cos(pi*n*t)*sinc(n/2);*n*t)*sinc(n/2);end end plot(t,cNplot(t,cN););%Gibbs phenomenonN=5N=15N=50N=500二、用数值积分分析非周期信号频谱二、用数值积分分析非周期信号频谱数值函数积分数值函数积分数值函数积分数值函数积分quad8quad8可用来计算非周期信号频谱可用来计算非周期信号频谱可用来计算非周期信号频谱可用来计算非周期信号频谱F F 是一个字符串，它表示被积函数的文件名。是一个字符串，它表示被积函数的文件名。是一个字符串，它表示被积函数的文件名。是一个字符串，它表示被积函数的文件名。a,ba,b 分别表示定积分的下限和上限分别表示定积分的下限和上限分别表示定积分的下限和上限分别表示定积分的下限和上限quad8quad8的返回是用自适应的返回是用自适应的返回是用自适应的返回是用自适应SimpsonSimpson算法得出的积分值算法得出的积分值算法得出的积分值算法得出的积分值y=quad8(F,a,b)y=quad8(F,a,b)例例例例3 3 试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱F F(j(j )=Sa)=Sa2 2(/2)/2)解解解解:图示三角波可表示为图示三角波可表示为图示三角波可表示为图示三角波可表示为三角波信号频谱的理论值为三角波信号频谱的理论值为三角波信号频谱的理论值为三角波信号频谱的理论值为例例例例3 3 试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱functionfunction y=y=sf1sf1(t,w);(t,w);y=(t=-1&t=-1&t=1).*(1-abs(t).*exp(-j*w*t);w=linspace(-6*pi,6*pi,512);w=linspace(-6*pi,6*pi,512);N=N=length(w);Flength(w);F=zeros(1,N);=zeros(1,N);forfor k=1:N k=1:N F(kF(k)=)=quad8quad8(sf1sf1,-1,1,w(k);,-1,1,w(k);endend figure(1);figure(1);plot(w,real(F);titleplot(w,real(F);title()xlabel(xlabel(omegaomega);ylabel();ylabel(F(jomegaF(jomega););figure(2);figure(2);plot(w,real(F)-sinc(w/2/pi).2);plot(w,real(F)-sinc(w/2/pi).2);xlabel(xlabel(omegaomega);title);title(计算误差计算误差计算误差计算误差););例例例例3 3 试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱试用数值方法近似计算三角波信号的频谱运行结果运行结果利用利用MATLAB进行系统频域分析进行系统频域分析连续系统频响特性的计算连续系统频响特性的计算 周期信号通过系统的响应周期信号通过系统的响应 离散系统频响特性的计算离散系统频响特性的计算一、连续系统频响特性的计算一、连续系统频响特性的计算b b 分子多项式系数分子多项式系数分子多项式系数分子多项式系数 a a 分母多项式系数分母多项式系数分母多项式系数分母多项式系数 w w 需计算的需计算的需计算的需计算的HH(j(j )的抽样点的抽样点的抽样点的抽样点 (数组数组数组数组w w中少需包含两个中少需包含两个中少需包含两个中少需包含两个 的抽样点的抽样点的抽样点的抽样点)。计算频响的计算频响的计算频响的计算频响的MATLABMATLAB函数函数函数函数H=H=freqs(b,a,wfreqs(b,a,w)一、连续系统频响特性的计算一、连续系统频响特性的计算例例例例1 1 三阶归一化的三阶归一化的Butterworth低通滤波器的系统函数为低通滤波器的系统函数为w=linspace(0,5,200);w=linspace(0,5,200);b=1;a=1 2 2 1;b=1;a=1 2 2 1;h=h=freqsfreqs(b,a,w(b,a,w););subplot(2,1,1);subplot(2,1,1);plot(w,abs(hplot(w,abs(h););subplot(2,1,2);subplot(2,1,2);plot(w,angle(hplot(w,angle(h););试画出试画出试画出试画出|HH(j(j )|)|和和和和 ()。一、连续系统频响特性的计算一、连续系统频响特性的计算三阶三阶三阶三阶ButterworthButterworth低通滤波器的幅度响应和相位响应低通滤波器的幅度响应和相位响应低通滤波器的幅度响应和相位响应低通滤波器的幅度响应和相位响应二、周期信号通过系统的响应二、周期信号通过系统的响应例例例例2 2 周期方波通过周期方波通过周期方波通过周期方波通过RCRC系统的响应。系统的响应。系统的响应。系统的响应。二、周期信号通过系统的响应二、周期信号通过系统的响应例例例例2 2 周期方波通过周期方波通过周期方波通过周期方波通过RCRC系统的响应。系统的响应。系统的响应。系统的响应。%p5_2 Periodic signal pass LTI system%p5_2 Periodic signal pass LTI system T=4;w0=2*pi/T;RC=0.1;T=4;w0=2*pi/T;RC=0.1;t=t=-6:0.01:6;N=51;6:0.01:6;N=51;c0=0.5;c0=0.5;xN=c0*ones(1,length(t);xN=c0*ones(1,length(t);%dc%dc forfor n=1:2:N n=1:2:N%even harmonics are zero%even harmonics are zero H=abs(1/(1+j*RC*w0*n);H=abs(1/(1+j*RC*w0*n);phi=angle(1/(1+j*RC*w0*n);phi=angle(1/(1+j*RC*w0*n);xNxN=xN+HxN+H*cos(w0*n*cos(w0*n*t+phit+phi)*)*sinc(nsinc(n*0.5);*0.5);end end plot(t,xNplot(t,xN););xlabel(xlabel(timetime RC=RC=,num2str(RC);grid;,num2str(RC);grid;set(gca,set(gca,xtickxtick,-5 5 -3 3 -1 0 1 3 5);1 0 1 3 5);二、周期信号通过系统的响应二、周期信号通过系统的响应例例例例2 2 周期方波通过周期方波通过周期方波通过周期方波通过RCRC系统的响应。系统的响应。系统的响应。系统的响应。三、离散系统频率响应的计算三、离散系统频率响应的计算计算频率响应的计算频率响应的计算频率响应的计算频率响应的MATLABMATLAB函数函数函数函数b b 分子的系数分子的系数分子的系数分子的系数 a a 分母系数分母系数分母系数分母系数w w 抽样的频率点抽样的频率点抽样的频率点抽样的频率点(至少至少至少至少2 2点点点点)，w w在在在在0202p p p p之间之间之间之间幅频特性幅频特性幅频特性幅频特性:absabs,相频特性相频特性相频特性相频特性:angleangle h=h=freqzfreqz(b,a,w(b,a,w)b=1;b=1;w=linspace(0,2*pi,512);w=linspace(0,2*pi,512);h2=h2=freqzfreqz(b,a2,w);(b,a2,w);plot(w/pi,plot(w/pi,absabs(h1),w/pi,(h1),w/pi,absabs(h2),(h2),:););a1=1 a1=1-0.9;a2=1 0.9;0.9;a2=1 0.9;h1=h1=freqzfreqz(b,a1,w);(b,a1,w);legend(legend(alpha=0.9alpha=0.9,alpha=alpha=-0.90.9););三、离散系统频响特性的计算三、离散系统频响特性的计算利用利用MATLAB进行连续系进行连续系统的复频域分析统的复频域分析部分分式展开的部分分式展开的MATLAB实现实现 H(s)的零极点与系统特性的的零极点与系统特性的MATLAB计算计算一、部分分式展开的一、部分分式展开的MATLAB实现实现 r,p,kr,p,k=residueresidue(num,den(num,den)num,dennum,den分别为分别为分别为分别为F F(s s)分子多项式和分母多项分子多项式和分母多项分子多项式和分母多项分子多项式和分母多项式的系数向量。式的系数向量。式的系数向量。式的系数向量。r r为部分分式的系数，为部分分式的系数，为部分分式的系数，为部分分式的系数，p p为极点，为极点，为极点，为极点，k k为多项式为多项式为多项式为多项式二、二、H(s)的零极点与系统特性的的零极点与系统特性的MATLAB计算计算计算多项式根计算多项式根计算多项式根计算多项式根rootsroots的函数可用于计算的函数可用于计算的函数可用于计算的函数可用于计算HH(s s)的零极点。的零极点。的零极点。的零极点。r=r=rootsroots(N(N)%)%计算多项式计算多项式计算多项式计算多项式N N的根的根的根的根 HH(s s)零极点分布图可用零极点分布图可用零极点分布图可用零极点分布图可用pzmappzmap函数画出，函数画出，函数画出，函数画出，调用形式为调用形式为调用形式为调用形式为pzmappzmap(sys(sys)表示画出表示画出表示画出表示画出syssys所描述系统的零极点图。所描述系统的零极点图。所描述系统的零极点图。所描述系统的零极点图。例例例例1 1 用部分分式展开法求用部分分式展开法求用部分分式展开法求用部分分式展开法求F F(s s)的反变换。的反变换。的反变换。的反变换。%program7_1%program7_1 format rat format rat%将结果数据以分数的形式输出将结果数据以分数的形式输出将结果数据以分数的形式输出将结果数据以分数的形式输出 num=1 2;den=1 4 3 0;num=1 2;den=1 4 3 0;r,pr,p=residueresidue(num,den(num,den)运行结果为运行结果为 r=-1/6 ,-1/2 ,2/3 p=-3 ,-1 ,0故故F F(s s)可展开为可展开为例例例例2 2 用部分分式展开法求用部分分式展开法求用部分分式展开法求用部分分式展开法求F F(s s)的反变换。的反变换。的反变换。的反变换。%program7_2%program7_2 num=2 3 0 5;num=2 3 0 5;den=den=convconv(1 1,1 1 2);(1 1,1 1 2);%将将将将因子相乘的形式转换成多项式的形式因子相乘的形式转换成多项式的形式因子相乘的形式转换成多项式的形式因子相乘的形式转换成多项式的形式 r,p,kr,p,k=residue(num,denresidue(num,den)magrmagr=absabs(r(r)%)%求求求求r r的模的模的模的模 angrangr=angleangle(r(r)%)%求求求求r r的相角的相角的相角的相角例例例例2 2 用部分分式展开法求用部分分式展开法求用部分分式展开法求用部分分式展开法求F F(s s)的反变换。的反变换。的反变换。的反变换。运行结果为运行结果为l lr=r=-2.0000+1.1339i,2.0000+1.1339i,-2.0000 2.0000-1.1339i,3.0000 1.1339i,3.0000 l lp=p=-0.5000+1.3229i,0.5000+1.3229i,-0.5000 0.5000-1.3229i,1.3229i,-1.0000 1.0000 l lk=2 k=2 l lmagrmagr=2.299,2.2991,3.0000 =2.299,2.2991,3.0000 l langrangr=2.6258,=2.6258,-2.6258,0 2.6258,0 故故F F(s s)可展开为可展开为例例例例3 3 试画出系统试画出系统 的零极点分布图，求的零极点分布图，求其单位冲激响应其单位冲激响应h(t)和频率响应和频率响应H(j)，并判断系统是否稳定。并判断系统是否稳定。num=1;den=1 2 2 1;num=1;den=1 2 2 1;sys=sys=tf tf(num,den(num,den););poles=poles=rootsroots(den(den)figure(1);figure(1);pzmappzmap(sys);(sys);t=0:0.02:10;t=0:0.02:10;h=h=impulseimpulse(num,den,t(num,den,t););figure(2);plot(t,h)figure(2);plot(t,h)title(title(ImpulseImpulse ResponeRespone)H,wH,w=freqsfreqs(num,den(num,den););figure(3);plot(w,abs(H)figure(3);plot(w,abs(H)xlabelxlabel(omegaomega)title(title(MagnituMagnitu 利用利用MATLAB进行离散系统的进行离散系统的Z域分析域分析部分分式展开的部分分式展开的MATLAB实现实现 H(z)的零极点与系统特性的的零极点与系统特性的MATLAB计算计算一、部分分式展开的一、部分分式展开的MATLAB实现实现 r,p,kr,p,k=residuezresiduez(num,den(num,den)num,dennum,den分别为分别为分别为分别为F F(z z)分子多项式和分母多项分子多项式和分母多项分子多项式和分母多项分子多项式和分母多项式的系数向量。式的系数向量。式的系数向量。式的系数向量。r r为部分分式的系数，为部分分式的系数，为部分分式的系数，为部分分式的系数，p p为极点，为极点，为极点，为极点，k k为多项式为多项式为多项式为多项式的系数。若的系数。若的系数。若的系数。若为真分式，则为真分式，则为真分式，则为真分式，则k k为零。为零。为零。为零。二、二、H(z)的零极点与系统特性的的零极点与系统特性的MATLAB计算计算利用利用利用利用tf2zptf2zp函数计算函数计算函数计算函数计算HH(z z)的零极点，调用形式为的零极点，调用形式为的零极点，调用形式为的零极点，调用形式为 z,p,kz,p,k=tf2zptf2zp(b,a)(b,a)HH(z z)零极点分布图可用零极点分布图可用零极点分布图可用零极点分布图可用zplanezplane函数画出，函数画出，函数画出，函数画出，调用形式为调用形式为调用形式为调用形式为zplanezplane(b,a(b,a)b b和和和和a a分别为分别为分别为分别为HH(z z)分子多项式和分母多项式的系数向量。分子多项式和分母多项式的系数向量。分子多项式和分母多项式的系数向量。分子多项式和分母多项式的系数向量。返回值返回值返回值返回值z z为零点、为零点、为零点、为零点、p p为极点、为极点、为极点、为极点、k k为增益常数。为增益常数。为增益常数。为增益常数。例例例例1 1 将将将将F F(z z)用部分分式展开。用部分分式展开。用部分分式展开。用部分分式展开。%program8_1%program8_1 num=18;den=18 3 num=18;den=18 3-4 4-1;1;r,p,kr,p,k=residuezresiduez(num,den(num,den)运行结果为运行结果为运行结果为运行结果为r=0.3600,0.2400,0.4000 r=0.3600,0.2400,0.4000 p=0.5000,p=0.5000,-0.3333,0.3333,-0.3333 0.3333 k=k=故故故故F F(z z)可展开为可展开为可展开为可展开为例例例例2 2 试画出系统试画出系统 的零极点的零极点分布图，求其单位冲激响应分布图，求其单位冲激响应hk和频率响应和频率响应H(ejW W)。%program 8_2%program 8_2 b=1 2 1;a=1 b=1 2 1;a=1-0.5 0.5-0.005 0.3;0.005 0.3;figure(1);figure(1);zplanezplane(b,a);(b,a);num=0 1 2 1;num=0 1 2 1;den=1 den=1-0.5 0.5-0.005 0.3;0.005 0.3;h=h=impzimpz(num,den(num,den););figure(2);stem(h)figure(2);stem(h)xlabel(xlabel(kk)title(title(ImpulseImpulse ResponeRespone)H,wH,w=freqzfreqz(num,den(num,den););figure(3);plot(w/pi,abs(H)figure(3);plot(w/pi,abs(H)xlabel(xlabel(FrequencyFrequency omega omega)title(title(MagnitudeMagnitude ResponeRespone)运行结果运行结果 利用利用MATLAB进行系统状态变量分析进行系统状态变量分析 微分方程到状微分方程到状态方程的方程的转换 状态方程系统函数矩阵状态方程系统函数矩阵H(s)的计算的计算 MATLAB求解连续系统的状态方程求解连续系统的状态方程 MATLAB求解离散系统的状态方程求解离散系统的状态方程一、微分方程到状态方程的转换一、微分方程到状态方程的转换 num,dennum,den 分别表示系统函数分别表示系统函数分别表示系统函数分别表示系统函数HH(s s)的分子的分子的分子的分子和分母多项式；和分母多项式；和分母多项式；和分母多项式；A,B,C,D A,B,C,D 分别为状态方程的矩阵。分别为状态方程的矩阵。分别为状态方程的矩阵。分别为状态方程的矩阵。A,B,C,D=A,B,C,D=tf2sstf2ss(num,den)(num,den)二、状态方程系统函数矩阵二、状态方程系统函数矩阵H(s)的计算的计算A,B,C,DA,B,C,D 分别表示状态方程的矩阵。分别表示状态方程的矩阵。分别表示状态方程的矩阵。分别表示状态方程的矩阵。KK 表示函数表示函数表示函数表示函数ss2tfss2tf计算的与第计算的与第计算的与第计算的与第k k个输入相关的系统个输入相关的系统个输入相关的系统个输入相关的系统函数，即函数，即函数，即函数，即HH(s s)的第的第的第的第k k列。列。列。列。numnum 表示表示表示表示H(s s)第第第第k k列的列的列的列的mm个元素的分子多项式个元素的分子多项式个元素的分子多项式个元素的分子多项式 denden 表示表示表示表示HH(s s)公共的分母多项式。公共的分母多项式。公共的分母多项式。公共的分母多项式。num,dennum,den=ss2tf ss2tf(A,B,C,D,kA,B,C,D,k)三、三、MATLAB求解连续求解连续系统的状态方程系统的状态方程sys 由函数由函数ss构造的状态方程模型构造的状态方程模型 t需计算的输出样本点需计算的输出样本点,t=0:dt:Tfinal f(:,k)系统第系统第k个输入在个输入在t上的抽样值上的抽样值 x0系统的初始状态系统的初始状态(可缺省可缺省)y(:,k)系统的第系统的第k个输出个输出 to实际计算时所用的样本点；实际计算时所用的样本点；x 系统的状态系统的状态 获得连续系统状态方程的计算机表示模型获得连续系统状态方程的计算机表示模型获得连续系统状态方程的计算机表示模型获得连续系统状态方程的计算机表示模型sys=ss(A,B,C,D)求解状态方程求解状态方程求解状态方程求解状态方程y,to,x=lsim(sys,f,t,x0)sys 由函数由函数ss构造的状态方程模型构造的状态方程模型 f(:,k)系统第系统第k个输入序列个输入序列 x0系统的初始状态系统的初始状态(可缺省可缺省)y(:,k)系统的第系统的第k个输出个输出 n序列的下标序列的下标;x系统的状态系统的状态获得离散状态方程的计算机表示模型获得离散状态方程的计算机表示模型获得离散状态方程的计算机表示模型获得离散状态方程的计算机表示模型sys=ss(A,B,C,D,)求解求解求解求解状态方程状态方程状态方程状态方程y,n,x=lsim(sys,f,x0)例例例例1 1 写出系统写出系统写出系统写出系统 的状态方程。的状态方程。的状态方程。的状态方程。由由由由A,B,C,D=A,B,C,D=tf2sstf2ss(1,1 5 10)(1,1 5 10)可得可得可得可得所以系统的状态方程为所以系统的状态方程为所以系统的状态方程为所以系统的状态方程为C C=0 =0 D D=0=0例例例例2 2 已知某连续系统的状态方程和输出方程为已知某连续系统的状态方程和输出方程为已知某连续系统的状态方程和输出方程为已知某连续系统的状态方程和输出方程为其初始状态和输入分别为其初始状态和输入分别为其初始状态和输入分别为其初始状态和输入分别为求该系统的系统函数矩阵求该系统的系统函数矩阵求该系统的系统函数矩阵求该系统的系统函数矩阵HH(s s)和输出。和输出。和输出。和输出。计算系统函数矩阵计算系统函数矩阵计算系统函数矩阵计算系统函数矩阵HH(s s)A=2 3;0 A=2 3;0-1;B=0 1;1 0;1;B=0 1;1 0;C=1 1;0 C=1 1;0-1;D=1 0;1 0;1;D=1 0;1 0;num1,den1=num1,den1=ss2tfss2tf(A,B,C,D,1)(A,B,C,D,1)num2,den2=num2,den2=ss2tfss2tf(A,B,C,D,2)(A,B,C,D,2)运行结果运行结果 num1=1 0 num1=1 0 -1 1 1 1 -2 0 2 0 den1=1 den1=1 -1 1 -2 2 num2=0 1 1 num2=0 1 1 0 0 0 0 0 0 den2=1 den2=1 -1 1 -2 2 所以系统函数矩阵所以系统函数矩阵HH(s s)为为计算输出计算输出计算输出计算输出%Program 9_1%Program 9_1 A=2 3;0 A=2 3;0-1;B=0 1;1 0;1;B=0 1;1 0;C=1 1;0 C=1 1;0-1;D=1 0;1 0;1;D=1 0;1 0;x0=2 x0=2-1;1;dtdt=0.01;t=0:dt:2;=0.01;t=0:dt:2;f(:,1)=ones(length(t),1);f(:,1)=ones(length(t),1);f(:,2)=exp(f(:,2)=exp(-3*t);3*t);sys=sys=ssss(A,B,C,D(A,B,C,D););y=y=lsimlsim(sys,f,t,x0);(sys,f,t,x0);subplot(2,1,1);plot(t,y(:,1),subplot(2,1,1);plot(t,y(:,1),rr););ylabel(ylabel(y1(t)y1(t);xlabel();xlabel(tt););subplot(2,1,2);plot(t,y(:,2);subplot(2,1,2);plot(t,y(:,2);ylabel(ylabel(y2(t)y2(t);xlabel();xlabel(tt););运行结果运行结果运行结果运行结果例例例例3 3 已知某离散系统的状态方程和输出方程为已知某离散系统的状态方程和输出方程为已知某离散系统的状态方程和输出方程为已知某离散系统的状态方程和输出方程为其初始状态和输入分别为其初始状态和输入分别为其初始状态和输入分别为其初始状态和输入分别为求该系统的输出。求该系统的输出。求该系统的输出。求该系统的输出。计算输出计算输出计算输出计算输出%Program 9_2%Program 9_2 A=0 1;A=0 1;-2 3;B=0;1;2 3;B=0;1;C=1 1;2 C=1 1;2 -1;D=zeros(2,1);1;D=zeros(2,1);x0=1;x0=1;-1;1;N=10;f=ones(1,N);N=10;f=ones(1,N);sys=sys=ssss(A,B,C,D(A,B,C,D,);,);y=y=lsimlsim(sys,f,x0);(sys,f,x0);subplot(2,1,1);y1=y(:,1);subplot(2,1,1);y1=y(:,1);stem(0:Nstem(0:N-1),y1);1),y1);xlabel(k);ylabel(y1);xlabel(k);ylabel(y1);subplot(2,1,2);y2=y(:,2);subplot(2,1,2);y2=y(:,2);stem(0:Nstem(0:N-1),y2);1),y2);xlabel(k);ylabel(y2);xlabel(k);ylabel(y2);运行结果运行结果运行结果运行结果