信號(hào)信息系統(tǒng)教案l28_matlab

1、信號(hào)與系統(tǒng),Signals and Systems,國家精品課程主教材、北京市精品教材《信號(hào)與系統(tǒng)》(第2版)陳后金，胡健，薛健清華大學(xué)出版社，2005年,利用MATLAB進(jìn)行信號(hào)與系統(tǒng)分析,,MATLAB簡介信號(hào)的MATLAB表示利用MATLAB進(jìn)行系統(tǒng)的時(shí)域分析利用MATLAB進(jìn)行信號(hào)的頻域分析利用MATLAB分析系統(tǒng)的頻率特性利用MATLAB進(jìn)行連續(xù)系統(tǒng)的S域分析利用MATLAB進(jìn)行離散系統(tǒng)的Z域分析利用MAT

2、LAB進(jìn)行系統(tǒng)的狀態(tài)變量分析,MATLAB簡介,,( Matrix Laboratory ),MATLAB的工作方式如何獲取幫助表達(dá)式——變量、數(shù)值、算數(shù)運(yùn)算符、關(guān)系運(yùn)算符、邏輯運(yùn)算符、冒號(hào)運(yùn)算符數(shù)組及其運(yùn)算函數(shù)文件循環(huán)（FOR、 WHILE 循環(huán)）基本繪圖語句,一、MATLAB的工作方式,（1）窗口命令方式（2）運(yùn)行以 .M 為擴(kuò)展名磁盤文件,工作方式舉例,%用plot函數(shù)畫一個(gè)方波t=[-1 0 0 1 1 3];

3、f=[0 0 1 1 0 0];plot(t,f);xlabel('t');ylabel('f(t)');axis([-1 3 0 2]);,直接在命令窗口輸入以上命令,建一個(gè)名字為my_file.M的文件，然后在命令窗口輸入文件名及回車。,二、獲取幫助,命令窗口輸入： help+函數(shù)名例如 help plot,三、表達(dá)式,不需要變量的類型說明變量名的第一個(gè)字符必須是字母變量名長度：不超過

4、31個(gè)字符大寫和小寫的字母視為不同的字符 例如：num_students = 25特殊變量： pi 表示圓周率，inf 表示無窮大，NaN(Not a Number)表示不定量，如0/0。,變量,三、表達(dá)式,數(shù)值,MATLAB用常規(guī)的十進(jìn)制表示數(shù)值 用i或j作為后綴來表示復(fù)數(shù)的虛部 例 1.235e5表示1.235?105，x=2+3jabs(x) 求復(fù)數(shù)x的模angle(x

5、) 求復(fù)數(shù)x的相角(弧度)real(x) 求復(fù)數(shù)x的實(shí)部imag(x) 求復(fù)數(shù)x的虛部conj(x) 求復(fù)數(shù)x的共軛,三、表達(dá)式,運(yùn)算符號(hào),算數(shù)運(yùn)算符,+ 加 - 減 * 乘 / 除 ^ 乘方 ‘ 矩陣的復(fù)共軛轉(zhuǎn)置,三、表達(dá)式,運(yùn)算符號(hào),邏輯運(yùn)算符,A & B邏輯與(and) A | B邏輯或(or)

6、 ~A 邏輯非(not),值為0時(shí)表示邏輯假(F)，其它任何非零值表示邏輯真。,三、表達(dá)式,運(yùn)算符號(hào),關(guān)系運(yùn)算符,A B 大于A = B 大于等于A == B 等于A ~= B不等于,三、表達(dá)式,運(yùn)算符號(hào),冒號(hào)運(yùn)算符,表達(dá)式 1:10 表示產(chǎn)生一個(gè)行向量，它的值為 1 2 3 4 5 6 7 8 9

7、 10表達(dá)式 10:-2:1 表示產(chǎn)生一個(gè)遞減的行向量，它的值為 10 8 6 4 2,四、數(shù)組,1. 數(shù)組的構(gòu)造,用冒號(hào):產(chǎn)生數(shù)組,例 x=2:5 產(chǎn)生一個(gè)數(shù)組，它的值為x(1)=2, x(2)=3, x(3)=4, x(4)=5,x=linspace(0,2,11)將區(qū)間[0，2]均勻抽樣11點(diǎn)作為數(shù)組x,給2維數(shù)組賦值時(shí)，用分號(hào)表示一行的結(jié)束，如

8、：z=[1 2; 3 4]。,用linspace產(chǎn)生數(shù)組,四、數(shù)組,1. 數(shù)組的構(gòu)造,MATLAB 提供了一些產(chǎn)生基本矩陣的函數(shù) zeros 產(chǎn)生矩陣元素全為0的矩陣 ones 產(chǎn)生矩陣元素全為1的矩陣 rand 產(chǎn)生(0,1)均勻分布隨機(jī)數(shù)矩陣 randn 產(chǎn)生正態(tài)分布隨機(jī)數(shù)矩陣,四、數(shù)組,2. 數(shù)組的運(yùn)算,數(shù)組和一個(gè)標(biāo)量相加或相

9、乘 例 y=x-1 z=3*x2個(gè)數(shù)組的對(duì)應(yīng)元素相乘除 .* ./ 例 z=x.*y確定數(shù)組大小的函數(shù) size(A) 返回值數(shù)組A的行數(shù)和列數(shù)（二維） length(B) 確定數(shù)組B的元素個(gè)數(shù)（一維）,五、函數(shù)文件,M文件的第一行包含function功能: 建立一個(gè)函數(shù)，可以同MATLAB的庫函數(shù)一樣使用。,五、函數(shù)文件,

10、例：編一個(gè)繪制圖示波形的函數(shù)。,function y=tri(t)y=[ abs(t)<=1].*(1-abs(t));,調(diào)用函數(shù)tri，并畫出它的波形,t=-2:0.05:2;plot(t,tri(t));,六、For 循環(huán),例: 編寫計(jì)算s=1+2+3+…+100的MATLAB程序s=0;,for n=1:100 s=s+n;end,七、While 循環(huán),s=0;n=1;eps=1e-6;,while

11、 1/(n*n) > eps s=s+1/(n*n); n=n+1; end,例: 計(jì)算 的值，且誤差小于10-6,fprintf('s=%.5f\n',s),八、plot函數(shù)——繪圖函數(shù)(continuous),t=linspace(0,4*pi,512);plot(t,sin(t),t,cos(t),'-.&

12、#39;);title('my figure');xlabel('t');legend('sin(t)','cos(t)');,,八、plot函數(shù)——繪圖函數(shù)(continuous),九、stem函數(shù)——繪圖函數(shù)(discrete),k=0:39;stem(k,cos(0.9*pi*k));title('cos(0.9\pik)');,,九、ste

13、m函數(shù)——繪圖函數(shù)(discrete),cos(0.9?k)波形,信號(hào)的MATLAB表示,,基本信號(hào)的MATLAB表示 指數(shù)信號(hào)Aeat 、指數(shù)序列ak 、抽樣函數(shù)Sa(t)、 正弦型信號(hào)、矩形脈沖信號(hào)、三角脈沖信號(hào)信號(hào)基本運(yùn)算的MATLAB實(shí)現(xiàn) 尺度變換、翻轉(zhuǎn)、時(shí)移、 相加、相乘、 差分與求和、微分與積分,一、基本信號(hào)的MATLAB表示,指數(shù)信號(hào)Aeat y = A*exp(a*

14、t);指數(shù)序列ak冪運(yùn)算a.^k實(shí)現(xiàn)正弦型信號(hào)內(nèi)部函數(shù)cos( ) 和sin( )抽樣函數(shù)Sa(t) sinc(t) 矩形脈沖信號(hào) y = rectpuls(t,width)三角波脈沖信號(hào)y = tripuls(t, width,skew),一、基本信號(hào)的MATLAB表示,,%decaying exponential,t=0:001:10;A=1;a=-0.4;ft=A*exp(a*t);pl

15、ot(t,ft),t=0:0.1:10;A=1;a=-0.4;ft=A*exp(a*t);stem(t,ft),一、基本信號(hào)的MATLAB表示,,% rectpuls,t=0:0.001:4;T=1;ft=rectpuls(t-2*T,T);plot(t,ft)axis([0,4,-0.5,1.5]),一、基本信號(hào)的MATLAB表示,,% tripuls,t=-3:0.001:3;ft=tripuls(t,4,0.5)

16、;plot(t,ft),ft=tripuls(t,4,1);,一、基本信號(hào)的MATLAB表示,,% unit impuls sequence,k=-50:50;delta=[zeros(1,50),1,zeros(1,50)];stem(k,delta),function [f,k]=impseq(k0,k1,k2)%產(chǎn)生 f[k]=delta(k-k0)；k1<=k<=k2k=[k1:k2];f=[(k-k0)=

17、=0];,k0=0;k1=-50;k2=50;[f,k]=impseq(k0,k1,k2);stem(k,f),一、基本信號(hào)的MATLAB表示,,% unit step sequence,k=-50:50;uk=[zeros(1,50), ones(1,51)];stem(k,uk),function [f,k]=stepseq(k0,k1,k2) %產(chǎn)生 f[k]=u(k-k0);k1=0];,k0=0;k1=-50;k2=

18、50;[f,k]=stepseq(k0,k1,k2);stem(k,f),二、信號(hào)基本運(yùn)算的MATLAB實(shí)現(xiàn),,t=-3:0.001:3;ft1=tripuls(2*t,4,0.5);subplot(2,1,1)plot(t,ft1)title('f(2t)')ft2=tripuls((2-2*t),4,0.5);subplot(2,1,2)plot(t,ft2)title('f(2-2t)&

19、#39;),1. 信號(hào)的尺度變換、翻轉(zhuǎn)、時(shí)移（平移）,已知三角波f(t)，用MATLAB畫出的f(2t)和f(2-2t) 波形,二、信號(hào)基本運(yùn)算的MATLAB實(shí)現(xiàn),,2. 信號(hào)的相加與相乘,相加用算術(shù)運(yùn)算符“+”實(shí)現(xiàn)相乘用數(shù)組運(yùn)算符“.*”實(shí)現(xiàn) 例:畫信號(hào)Aeatcos(w0t+f)的波形 t=0:0.001:8; A=1; a=-0.4; w0=2*pi;phi=0; ft1=A*exp(a*t

20、).*sin(w0*t+phi); plot(t,ft1),二、信號(hào)基本運(yùn)算的MATLAB實(shí)現(xiàn),3. 離散序列的差分與求和 連續(xù)信號(hào)的微分與積分,差分y=diff(f);求和y=sum(f(k1:k2));,微分 y=diff(f)/h; h為數(shù)值計(jì)算所取時(shí)間間隔,定積分 quad(‘function_name’,a,b);,function_name為被積函數(shù)名，a和b指定積分區(qū)間。,,二、

21、信號(hào)基本運(yùn)算的MATLAB實(shí)現(xiàn),3. 離散序列的差分與求和 連續(xù)信號(hào)的微分與積分,例:已知三角波f(t)，畫出其微分與積分的波形,%differentiation h=0.001;t= -3:h:3;y1=diff(f2_2(t))*1/h;plot(t(1:length(t)-1),y1),%integrationt= -3:0.1:3;for x=1:length(t) y2(x)=quad('f2_

22、2', -3,t(x));endplot(t,y2),三角波f(t)微分與積分的波形,利用MATLAB進(jìn)行系統(tǒng)的時(shí)域分析,,連續(xù)時(shí)間系統(tǒng)零狀態(tài)響應(yīng)的求解連續(xù)時(shí)間系統(tǒng)沖激響應(yīng)和階躍響應(yīng)的求解離散時(shí)間系統(tǒng)零狀態(tài)響應(yīng)的求解離散時(shí)間系統(tǒng)單位脈沖響應(yīng)的求解離散卷積的計(jì)算,,一、連續(xù)時(shí)間系統(tǒng)零狀態(tài)響應(yīng)的求解,t 表示計(jì)算系統(tǒng)響應(yīng)的抽樣點(diǎn)向量,a=[a3, a2, a1, a0]; b=[b3, b2, b1, b0];

23、sys=tf(b,a),y=lsim(sys,f,t),sys=tf(b,a),b和a分別為微分方程右端和左端各項(xiàng)的系數(shù)向量,f 是系統(tǒng)輸入信號(hào)向量，,sys 是LTI系統(tǒng)模型，借助tf函數(shù)獲得,二、連續(xù)時(shí)間系統(tǒng)沖激響應(yīng)和階躍響應(yīng)的求解,連續(xù)時(shí)間系統(tǒng)沖激響應(yīng)可用impulse函數(shù)直接求出，其調(diào)用形式為,y=impulse(sys,t),連續(xù)時(shí)間系統(tǒng)階躍響應(yīng)可用step函數(shù)直接求出，其調(diào)用形式為,y=step(sys,t),t 表示

24、計(jì)算系統(tǒng)響應(yīng)的抽樣點(diǎn)向量sys 是LTI系統(tǒng)模型,三、離散時(shí)間系統(tǒng)零狀態(tài)響應(yīng)的求解,b , a 分別是差分方程左、右端的系數(shù)向量,b=[b0,b1,b2,??,bM];a=[a0,a1,a2, ??,aN];,可用MATLAB表示為,y=filter(b,a,f),f 表示輸入序列, y 表示輸出序列,四、離散時(shí)間系統(tǒng)單位脈沖響應(yīng)的求解,b, a 分別是差分方程左、右端的系數(shù)向量k 表示輸出序列的取值范圍h 就是單位脈沖

25、響應(yīng),h=impz(b,a,k),五、離散卷積的計(jì)算,例：(s3+2s+3)(s2+3s+2)可用下面MATLAB語句求出 a =[1,0,2,3]; b =[1,3,2]; c=conv(a,b),c=conv(a,b),式中a,b為待卷積兩序列的向量表示，c是卷積結(jié)果。,conv函數(shù)也可用于計(jì)算兩個(gè)多項(xiàng)式的積,例1 求系統(tǒng) y"(t)+2y'(t)+100y(t)=10f(t) 的零狀

26、態(tài)響應(yīng)，已知f(t)=sin(2pt) u(t)。,%program3_1微分方程求解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 求系統(tǒng) y" (t)+2y 

27、9; (t)+100y(t)=10f(t) 的零狀態(tài)響應(yīng)，已知f(t) =d (t) 。,%program3_2連續(xù)時(shí)間系統(tǒng)的沖激響應(yīng)ts=0;te=5;dt=0.01;sys=tf([10],[1 2 100]);t=ts:dt:te;y=impulse(sys,t);plot(t,y);xlabel('Time(sec)')ylabel('h(t)'),例3 分析噪聲干擾的信號(hào)f[k]

28、=s[k]+d[k]通過M點(diǎn)滑動(dòng)平均系統(tǒng)的響應(yīng)，其中s[k]=(2k)0.9k是原始信號(hào)，d[k]是噪聲。,R =51 ; d = rand(1,R) - 0.5;k=0:R-1;s=2*k.*(0.9.^k); f=s+d;figure(1); plot(k,d,'r-.',k,s,'b--',k,f,'g-');M =5; b = ones(M,1)/M; a = 1;y

29、 = filter(b,a,f);figure(2); plot(k,s,'b--',k,y,'r-');,例3 分析噪聲干擾的信號(hào)f[k]=s[k]+d[k]通過M點(diǎn)滑動(dòng)平均系統(tǒng)的響應(yīng)，其中s[k]=(2k)0.9k是原始信號(hào)，d[k]是噪聲。,噪聲干擾信號(hào)f[k]=s[k]+d[k]通過M點(diǎn)滑動(dòng)平均系統(tǒng)的響應(yīng),例4 求系統(tǒng)y[k]+3y[k-1]+2y[k-1]=10f [k]的單位脈沖響

30、應(yīng)。,% program 3_4 離散系統(tǒng)的單位脈沖響應(yīng)k=0:10;a=[1 3 2];b=[1]; h=impz(b,a,k);stem(k,h),例5 計(jì)算x[k]* y[k]并畫出卷積結(jié)果，已知x[k]={1,2,3,4; k=0,1,2,3}，y[k]={1,1,1,1,1; k=0,1,2,3,4} 。,% program 3_5x=[1,2,3,4]; y=[1,1,1,1,1]; z=conv(x,

31、y);N=length(z);stem(0:N-1,z);,利用MATLAB進(jìn)行信號(hào)的頻域分析,,周期信號(hào)頻譜的MATLAB實(shí)現(xiàn)用數(shù)值積分分析非周期信號(hào)頻譜,一、周期信號(hào)頻譜的MATLAB實(shí)現(xiàn),頻譜Fn一般為復(fù)數(shù)，可分別利用abs和angle函數(shù)獲得其幅頻特性和相頻特性。,其調(diào)用格式分別為,x=abs(Fn) y=angle(Fn),周期信號(hào)的頻譜Fn 為離散信號(hào)，可以用stem畫出其頻譜圖。,例

32、1 試用MATLAB畫出圖示周期三角波信號(hào)的頻譜。,解：周期信號(hào)的頻譜為,畫三角波信號(hào)頻譜的MATLAB程序,N=8;n1= -N:-1; %計(jì)算n=-N到-1的Fourier系數(shù)c1= -4*j*sin(n1*pi/2)/pi^2./n1.^2;c0=0; %計(jì)算n=0時(shí)的Fourier系數(shù)n2=1:N; %計(jì)算n=1到N的Fourier系數(shù)c2= -4*j*sin(n2*pi/2)/pi^2./n2.^2;cn=[c1

33、 c0 c2];n= -N:N;subplot(2,1,1);stem(n,abs(cn));ylabel('Cn的幅度');subplot(2,1,2);stem(n,angle(cn));ylabel('Cn的相位');xlabel('\omega/\omega0');,程序運(yùn)行結(jié)果,例2 求周期矩形脈沖的Fourier級(jí)數(shù)表示式。并用MATLAB求出由前N項(xiàng)Four

34、ier級(jí)數(shù)系數(shù)得出的信號(hào)近似波形。,取A=1, T=2, t=1, w0=p,% Gibbs phenomenon,t=-2:0.001:2;N=input('Number of harmonics= ');c0=0.5;cN=c0*ones(1,length(t)); %dc componentfor n=0:2:N-1 % even harmonics are zero cN=cN+

35、cos(pi*n*t)*sinc(n/2);endplot(t,xN);,% Gibbs phenomenon,N=5,N=15,N=50,N=500,二、用數(shù)值積分分析非周期信號(hào)頻譜,數(shù)值函數(shù)積分quad8可用來計(jì)算非周期信號(hào)頻譜,F 是一個(gè)字符串，它表示被積函數(shù)的文件名。,a,b 分別表示定積分的下限和上限,quad8的返回是用自適應(yīng)Simpson算法得出的積分值,y = quad8('F',a,b),例3

36、 試用數(shù)值方法近似計(jì)算三角波信號(hào)的頻譜,F(jw)= Sa2(w / 2),解: 圖示三角波可表示為,三角波信號(hào)頻譜的理論值為,例3 試用數(shù)值方法近似計(jì)算三角波信號(hào)的頻譜,function y=sf1(t,w);y=(t>=-1 & t<=1).*(1-abs(t)).*exp(-j*w*t);,w=linspace(-6*pi,6*pi,512);N=length(w);F=zeros(1,N);

37、for k=1:N F(k)=quad8('sf1',-1,1,[],[],w(k));endfigure(1);plot(w,real(F));title('')xlabel('\omega');ylabel('F(j\omega)');figure(2);plot(w,real(F)-sinc(w/2/pi).^2);xlabel('\om

38、ega');title('計(jì)算誤差');,例3 試用數(shù)值方法近似計(jì)算三角波信號(hào)的頻譜,運(yùn)行結(jié)果,利用MATLAB進(jìn)行系統(tǒng)頻域分析,,連續(xù)系統(tǒng)頻響特性的計(jì)算 周期信號(hào)通過系統(tǒng)的響應(yīng) 離散系統(tǒng)頻響特性的計(jì)算,一、連續(xù)系統(tǒng)頻響特性的計(jì)算,b 分子多項(xiàng)式系數(shù) a 分母多項(xiàng)式系數(shù) w 需計(jì)算的H(jw)的抽樣點(diǎn) (數(shù)組w中少需包含兩個(gè)w的抽樣點(diǎn))。,計(jì)算頻響的MATLAB函數(shù),H

39、=freqs(b,a,w),,一、連續(xù)系統(tǒng)頻響特性的計(jì)算,例1 三階歸一化的Butterworth低通濾波器的系統(tǒng)函數(shù)為,w=linspace(0,5,200);b=[1];a=[1 2 2 1];h=freqs(b,a,w);subplot(2,1,1);plot(w,abs(h));subplot(2,1,2);plot(w,angle(h));,試畫出|H(jw)| 和?(w)。,,一、連續(xù)系統(tǒng)頻響特性的計(jì)算,三階Bu

40、tterworth低通濾波器的幅度響應(yīng)和相位響應(yīng),,二、周期信號(hào)通過系統(tǒng)的響應(yīng),例2 周期方波通過RC系統(tǒng)的響應(yīng)。,,二、周期信號(hào)通過系統(tǒng)的響應(yīng),例2 周期方波通過RC系統(tǒng)的響應(yīng)。,%p5_2 Periodic signal pass LTI system,T=4;w0=2*pi/T;RC=0.1;t= -6:0.01:6;N=51;c0=0.5;xN=c0*ones(1,length(t)); %dc for n=1

41、:2:N % even harmonics are zero H=abs(1/(1+j*RC*w0*n)); phi=angle(1/(1+j*RC*w0*n)); xN=xN+H*cos(w0*n*t+phi)*sinc(n*0.5);endplot(t,xN);xlabel(['time RC=',num2str(RC)]);grid;set(gca,'xtick',[-

42、5 -3 -1 0 1 3 5]);,,二、周期信號(hào)通過系統(tǒng)的響應(yīng),例2 周期方波通過RC系統(tǒng)的響應(yīng)。,三、離散系統(tǒng)頻率響應(yīng)的計(jì)算,計(jì)算頻率響應(yīng)的MATLAB函數(shù),b 分子的系數(shù) a 分母系數(shù),w 抽樣的頻率點(diǎn)(至少2點(diǎn))， w在0~2p之間,幅頻特性: abs, 相頻特性: angle,h = freqz(b,a,w),三、離散系統(tǒng)頻響特性的計(jì)算,,b=[1];,a1=[1 -0.9]

43、; a2=[1 0.9];,w=linspace(0,2*pi,512);,h1=freqz(b,a1,w);,h2=freqz(b,a2,w);,plot(w/pi,abs(h1),w/pi,abs(h2),':');,legend('\alpha=0.9','\alpha=-0.9');,三、離散系統(tǒng)頻響特性的計(jì)算,,利用MATLAB進(jìn)行連續(xù)系統(tǒng)的復(fù)頻域分析,,部分分式展開的MATLA

44、B實(shí)現(xiàn)H(s)的零極點(diǎn)與系統(tǒng)特性的MATLAB計(jì)算,一、部分分式展開的MATLAB實(shí)現(xiàn),[r,p,k]=residue(num,den),num,den分別為F(s)分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)向量。,r為部分分式的系數(shù)，p為極點(diǎn)，k為多項(xiàng)式的系數(shù)。若為真分式，則k為零。,二、H(s)的零極點(diǎn)與系統(tǒng)特性的MATLAB計(jì)算,計(jì)算多項(xiàng)式根roots的函數(shù)可用于計(jì)算H(s)的零極點(diǎn)。,r=roots(N) %計(jì)算多項(xiàng)式N的根,H(s)零極點(diǎn)

45、分布圖可用pzmap函數(shù)畫出，調(diào)用形式為,pzmap(sys),表示畫出sys所描述系統(tǒng)的零極點(diǎn)圖。,例1 用部分分式展開法求F(s)的反變換。,%program7_1format rat %將結(jié)果數(shù)據(jù)以分?jǐn)?shù)的形式輸出num=[1 2]; den=[1 4 3 0]; [r,p]=residue(num,den),運(yùn)行結(jié)果為 r = -1/6 ,-1/2 ,2/3 p = -3 ,-1 ,0,故F(s

46、)可展開為,例2 用部分分式展開法求F(s)的反變換。,% program7_2num=[2 3 0 5];den=conv([1 1],[1 1 2]); %將因子相乘的形式轉(zhuǎn)換成多項(xiàng)式的形式[r,p,k]=residue(num,den)magr=abs(r) %求r的模angr=angle(r) %求r的相角,例2 用部分分式展開法求F(s)的反變換。,運(yùn)行結(jié)果為,r =-2.0000 + 1.1339i, -

47、2.0000 - 1.1339i, 3.0000p =-0.5000 + 1.3229i, -0.5000 - 1.3229i, -1.0000k =2magr =2.299, 2.2991, 3.0000 angr =2.6258, -2.6258, 0,故F(s)可展開為,例3 試畫出系統(tǒng) 的零極點(diǎn)分布圖，求其單位沖激響應(yīng)h(t)和頻率響應(yīng)H(

48、j?)，并判斷系統(tǒng)是否穩(wěn)定。,num=[1];den=[1 2 2 1];sys=tf(num,den); poles=roots(den)figure(1);pzmap(sys);t=0:0.02:10;h=impulse(num,den,t);figure(2);plot(t,h)title('Impulse Respone')[H,w]=freqs(num,den);figure(3);plot(

49、w,abs(H))xlabel('\omega')title('Magnitude Respone'),運(yùn)行結(jié)果,利用MATLAB進(jìn)行離散系統(tǒng)的Z域分析,,部分分式展開的MATLAB實(shí)現(xiàn)H(z)的零極點(diǎn)與系統(tǒng)特性的MATLAB計(jì)算,一、部分分式展開的MATLAB實(shí)現(xiàn),[r,p,k]=residuez(num,den),num,den分別為F(z)分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)向量。,r為部分分式的系數(shù)

50、，p為極點(diǎn)，k為多項(xiàng)式的系數(shù)。若為真分式，則k為零。,二、H(z)的零極點(diǎn)與系統(tǒng)特性的MATLAB計(jì)算,利用tf2zp函數(shù)計(jì)算H(z)的零極點(diǎn)，調(diào)用形式為,[z,p,k]=tf2zp(b,a),H(z)零極點(diǎn)分布圖可用zplane函數(shù)畫出，調(diào)用形式為,zplane(b,a),b和a分別為H(z)分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)向量。,返回值z為零點(diǎn)、p為極點(diǎn)、 k為增益常數(shù)。,例1 將F(z)用部分分式展開。,%program8_1nu

51、m = [18]; den = [18 3 -4 -1];[r,p,k] = residuez(num,den),運(yùn)行結(jié)果為r =0.3600 , 0.2400 , 0.4000 p =0.5000 , -0.3333 , -0.3333 k =[],故F(z)可展開為,例2 試畫出系統(tǒng) 的零極點(diǎn)分布

52、圖，求其單位沖激響應(yīng)h[k]和頻率響應(yīng)H(ejW) 。,% program 8_2b =[1 2 1];a =[1 -0.5 -0.005 0.3];figure(1);zplane(b,a);num=[0 1 2 1];den=[1 -0.5 -0.005 0.3];h=impz(num,den);figure(2);stem(h)xlabel('k')title('Impulse Respon

53、e')[H,w]=freqz(num,den);figure(3);plot(w/pi,abs(H))xlabel('Frequency \omega')title('Magnitude Respone'),運(yùn)行結(jié)果,利用MATLAB進(jìn)行系統(tǒng)狀態(tài)變量分析,,微分方程到狀態(tài)方程的轉(zhuǎn)換 狀態(tài)方程系統(tǒng)函數(shù)矩陣H(s)的計(jì)算 MATLAB求解連續(xù)系統(tǒng)的狀態(tài)方程 MATLAB求解離散系統(tǒng)的狀態(tài)

54、方程,一、微分方程到狀態(tài)方程的轉(zhuǎn)換,num,den 分別表示系統(tǒng)函數(shù)H(s)的分子和分母多項(xiàng)式； A,B,C,D 分別為狀態(tài)方程的矩陣。,[A,B,C,D]= tf2ss(num,den),二、狀態(tài)方程系統(tǒng)函數(shù)矩陣H(s)的計(jì)算,A, B, C, D 分別表示狀態(tài)方程的矩陣。K 表示函數(shù)ss2tf計(jì)算的與第k個(gè)輸入相關(guān)的系統(tǒng)函數(shù)，即H(s)的第k列。num 表示H(s)第k列的m個(gè)元素的分子多項(xiàng)式den

55、表示H(s)公共的分母多項(xiàng)式。,[num,den]= ss2tf (A,B,C,D,k),,三、MATLAB求解連續(xù)系統(tǒng)的狀態(tài)方程,sys 由函數(shù)ss構(gòu)造的狀態(tài)方程模型t 需計(jì)算的輸出樣本點(diǎn), t=0:dt:Tfinalf(:,k) 系統(tǒng)第k個(gè)輸入在t上的抽樣值x0 系統(tǒng)的初始狀態(tài)(可缺省)y(:,k) 系統(tǒng)的第k個(gè)輸出 to 實(shí)際計(jì)算時(shí)所用的樣本點(diǎn)； x 系統(tǒng)的狀態(tài),獲得連續(xù)

56、系統(tǒng)狀態(tài)方程的計(jì)算機(jī)表示模型,sys = ss(A,B,C,D),求解狀態(tài)方程,[y,to,x]=lsim(sys,f,t,x0),,四、MATLAB求解離散系統(tǒng)的狀態(tài)方程,sys 由函數(shù)ss構(gòu)造的狀態(tài)方程模型f(:,k)系統(tǒng)第k個(gè)輸入序列x0 系統(tǒng)的初始狀態(tài)(可缺省)y(:,k)系統(tǒng)的第k個(gè)輸出n 序列的下標(biāo); x 系統(tǒng)的狀態(tài),獲得離散狀態(tài)方程的計(jì)算機(jī)表示模型,sys

57、 = ss(A,B,C,D,[ ]),求解狀態(tài)方程,[y,n,x]=lsim(sys,f,[ ],x0),例1 寫出系統(tǒng) 的狀態(tài)方程。,由[A,B,C,D]=tf2ss([1],[1 5 10])可得,所以系統(tǒng)的狀態(tài)方程為,C=0 D=0,例2 已知某連續(xù)系統(tǒng)的狀態(tài)方程和輸出方程為,其初始狀態(tài)和輸入分別為,求該

58、系統(tǒng)的系統(tǒng)函數(shù)矩陣H(s) 和輸出。,計(jì)算系統(tǒng)函數(shù)矩陣H(s),A=[2 3;0 -1];B=[0 1; 1 0];C=[1 1; 0 -1];D=[1 0; 1 0];[num1,den1]=ss2tf(A,B,C,D,1)[num2,den2]=ss2tf(A,B,C,D,2)運(yùn)行結(jié)果num1 =1 0 -1 1 -2 0den1 =1 -1 -2num

59、2 =0 1 1 0 0 0den2 =1 -1 -2所以系統(tǒng)函數(shù)矩陣H(s)為,計(jì)算輸出,%Program 9_1A=[2 3;0 -1];B=[0 1; 1 0];C=[1 1; 0 -1];D=[1 0; 1 0];x0=[2 -1];dt=0.01;t=0:dt:2; f(:,1)=ones(length(t),1);f(:,2)=exp(-3

60、*t)';sys=ss(A,B,C,D);y=lsim(sys,f,t,x0);subplot(2,1,1);plot(t,y(:,1),'r');ylabel('y1(t)');xlabel('t');subplot(2,1,2);plot(t,y(:,2));ylabel('y2(t)');xlabel('t');,運(yùn)行結(jié)果,例3

61、 已知某離散系統(tǒng)的狀態(tài)方程和輸出方程為,其初始狀態(tài)和輸入分別為,求該系統(tǒng)的輸出。,計(jì)算輸出,%Program 9_2A=[0 1; -2 3];B=[0;1];C=[1 1; 2 -1];D=zeros(2,1);x0=[1; -1];N=10;f=ones(1,N);sys=ss(A,B,C,D,[]);y=lsim(sys,f,[],x0);subplot(2,1,1);y1=y(:,1)';stem((0