信號與系統(tǒng)陳后金MATLABPPT學習教案

1、會計學1信號與系統(tǒng)陳后金信號與系統(tǒng)陳后金MATLAB第1頁/共100頁第2頁/共100頁（1）窗口命令方式（2）運行以 .M 為擴展名磁盤文件第3頁/共100頁%用plot函數畫一個方波t=-1 0 0 1 1 3;x=0 0 1 1 0 0;plot(t,x);xlabel(t);ylabel(x(t);axis(-1 3 0 2);直接在命令窗口輸入以上命令建一個名字為的文件，然后在命令窗口輸入文件名及回車。第4頁/共100頁命令窗口輸入： help+函數名例如 help plot第5頁/共100頁不需要變量的類型說明變量名的第一個字符必須是字母變量名長度：不超過31個字符大寫和小寫的字母

2、視為不同的字符 例如：num_students = 25特殊變量： pi 表示圓周率，inf 表示無窮大，NaN(Not a Number)表示不定量，如0/0。第6頁/共100頁 MATLAB用常規(guī)的十進制表示數值 用i或j作為后綴來表示復數的虛部 例 表示105，x=2+3jabs(x) 求復數x的模angle(x) 求復數x的相角(弧度)real(x) 求復數x的實部imag(x) 求復數x的虛部conj(x) 求復數x的共軛第7頁/共100頁u + 加u - 減u * 乘u 除u 乘方u 矩陣的復共軛轉置第8頁/共100頁u A & B邏輯與(and)u A | B邏輯或(or

3、)u A 邏輯非(not) 值為0時表示邏輯假(F)，其它任何非零值表示邏輯真。第9頁/共100頁u A B 大于uA = B 大于等于uA = B 等于uA = B 不等于第10頁/共100頁u表達式 1:10 表示產生一個行向量，它的值為 1 2 3 4 5 6 7 8 9 10u表達式 10:-2:1 表示產生一個遞減的行向量，它的值為 10 8 6 4 2 第11頁/共100頁用冒號:產生數組例 x=2:5 產生一個數組，它的值為x(1)=2, x(2)=3, x(3)=4, x(4)=5x=linspace(0,2,11)將區(qū)間0，2均勻抽樣11點作為數組x 給2維數組賦值時，用分號

4、表示一行的結束，如：z=1 2; 3 4。用linspace產生數組第12頁/共100頁MATLAB 提供了一些產生基本矩陣的函數 zeros 產生矩陣元素全為0的矩陣 ones 產生矩陣元素全為1的矩陣 rand 產生(0,1)均勻分布隨機數矩陣 randn 產生正態(tài)分布隨機數矩陣第13頁/共100頁size(A) 返回值數組A的行數和列數（二維） length(B) 確定數組B的元素個數（一維）第14頁/共100頁第15頁/共100頁例：編一個繪制圖示波形的函數。function y=tri(t)y= abs(t) eps s=s+1/(n*n); n=n+1; end+222312111

5、s例例: 計算 的值，且誤差小于10-6fprintf(s=%.5fn,s)第18頁/共100頁t=linspace(0,4*pi,512);plot(t,sin(t),t,cos(t),-.);title(my figure);xlabel(t);legend(sin(t),cos(t);第19頁/共100頁第20頁/共100頁k=0:39;stem(k,cos(0.9*pi*k);title(cos(0.9pik);第21頁/共100頁0510152025303540-1-0.8-0.6-0.4-0.60.81cos(0.9k)k)波形波形第22頁/共100頁指數信號A

6、eat 、指數序列ak 、抽樣函數Sa(t)、 正弦型信號、矩形脈沖信號、三角脈沖信號尺度變換、翻轉、時移、尺度變換、翻轉、時移、 相加、相乘、相加、相乘、 差分與求和、微分與積分差分與求和、微分與積分第23頁/共100頁指數信號Aeat y = A*exp(a*t);指數序列ak冪運算a.k實現正弦型信號內部函數cos( ) 和sin( )抽樣函數Sa(t) sinc(t) 矩形脈沖信號 y = rectpuls(t,width)三角波脈沖信號y = tripuls(t, width,skew) 第24頁/共100頁%decaying exponentialt=0:001:10;A=1;a=

7、-0.4;ft=A*exp(a*t);plot(t,ft)01234567891001t=0:0.1:10;A=1;a=-0.4;ft=A*exp(a*t);stem(t,ft)01234567891001第25頁/共100頁% rectpulst=0:0.001:4;T=1;ft=rectpuls(t-2*T,T);plot(t,ft)axis(0,4,-0.5,1.5)00.511.522.533.54-0.500.511.5第26頁/共100頁-3-2-1012300.10.20.

8、1-3-2-1012301第27頁/共100頁-50-40-30-20-100102030405001第28頁/共100頁-50-40-30-20-100102030405001第29頁/共100頁第30頁/共100頁-3-2-1012301第31頁/共100頁相加相加用算術運算符用算術運算符“+”實現實現相乘相乘用數組運算符用數組運算符“.*

9、”實現實現 例:畫信號Aeatcos(w0t+f)的波形 t=0:0.001:8; A=1; a=-0.4; w0=2*pi;phi=0; ft1=A*exp(a*t).*sin(w0*t+phi); plot(t,ft1)012345678-0.8-0.6-0.4-0.60.81第32頁/共100頁第33頁/共100頁:已知三角波x(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

10、 x=1:length(t) y2(x)=quad(f2_2, -3,t(x);endplot(t,y2)第34頁/共100頁-3-2-1012301第35頁/共100頁第36頁/共100頁t 表示計算系統(tǒng)響應的抽樣點向量)()( )( )( )()( )( )( 01230123txbtxbtxbtxbtyatyatyatya+a=a3, a2, a1, a0; b=b3, b2, b1, b0; sys=tf(b,a)y=lsim(sys,x,t)sys=tf(b,a)b和a分別為微分方程右端和左端各項的系數向量x 是系統(tǒng)輸入信號向量，

11、sys 是LTI系統(tǒng)模型，借助tf函數獲得第37頁/共100頁 連續(xù)時間系統(tǒng)沖激響應可用impulse函數直接求出，其調用形式為y=impulse(sys, t) 連續(xù)時間系統(tǒng)階躍響應可用step函數直接求出，其調用形式為y=step(sys, t)t 表示計算系統(tǒng)響應的抽樣點向量sys 是LTI系統(tǒng)模型第38頁/共100頁b , a 分別是差分方程左、右端的系數向量00jkxbikyajmjini-b=b0,b1,b2,bM;a=a0,a1,a2, ,aN;可用MATLAB表示為y=filter(b,a,x)x 表示輸入序列, y 表示輸出序列第39頁/共100頁b, a 分別是差分方程左、

12、右端的系數向量k 表示輸出序列的取值范圍h 就是單位脈沖響應h=impz(b,a,k)第40頁/共100頁例：(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是卷積結果。conv函數也可用于計算兩個多項式的積第41頁/共100頁%program3_1微分方程求解ts=0;te=5;dt=0.01;sys=tf(10,1 2 100);t=ts:dt:te;x=sin(2*pi*t);y=lsim(sys,x,t);plot(t,y);xlabel(Ti

13、me(sec)ylabel(y(t)012345-0.2-Time(sec)y(t)第42頁/共100頁%program3_2連續(xù)時間系統(tǒ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)012345-1-0.500.51Time(sec)h(t)第43頁/共100頁R =51 ; d = rand(1,R) - 0.5;k=0:R-1;s=2*k.*(0.9.k); x=s+d;figure(1); plot(

14、k,d,r-.,k,s,b-,k,x,g-);M =5; b = ones(M,1)/M; a = 1;y = filter(b,a,x);figure(2); plot(k,s,b-,k,y,r-);110nkxMkyMn-第44頁/共100頁例例3 3 分析噪聲干擾的信號xk=sk+dk通過M點滑動平均系統(tǒng)的響應，其中sk=(2kk是原始信號，dk是噪聲。110nkxMkyMn- 0 10 20 30 40 50 -2 0 2 4 6 8 Time index k dk sk xk 01020304050-202468Time index kskyk噪聲干擾信號xk=sk+dk通過M點滑動

15、平均系統(tǒng)的響應第45頁/共100頁% program 3_4 離散系統(tǒng)的單位脈沖響應k=0:10;a=1 3 2;b=10; h=impz(b,a,k);stem(k,h)012345678910-2000-10000100020003000單位脈沖響應第46頁/共100頁例例5 計算xk* yk并畫出卷積結果，已知xk=1,2,3,4; k=0,1,2,3，yk=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,y);N=length(z);stem(0:N-1,z);012345670246810第4

16、7頁/共100頁第48頁/共100頁 頻譜Cn一般為復數，可分別利用abs和angle函數獲得其幅度頻譜和相位頻譜。其調用格式分別為x=abs(Cn) y=angle(Cn) 周期信號的頻譜Cn 為離散信號，可以用stem畫出其頻譜圖。第49頁/共100頁例例1 試用MATLAB畫出圖示周期三角波信號的頻譜。解：解：周期信號的頻譜為-000 )2sin(j422nnnnCn第50頁/共100頁N=8;n1= -N:-1; %計算n=-N到-1的Fourier系數c1= -4*j*sin(n1*pi/2)/pi2./n1.2;c0=0; %計算n=0時的Fourier系數n2=1:N; %計算n

17、=1到N的Fourier系數c2= -4*j*sin(n2*pi/2)/pi2./n2.2;cn=c1 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);第51頁/共100頁-8-6-4-202468-2-1012Cn的相位w /w 0-8-6-4-20246800.40.5Cn的幅度第52頁/共100頁例例2 求周期矩形脈沖的Fourier級數表示式。并用MATLAB求出由前N

18、項Fourier級數系數得出的信號近似波形。)2/(Sa0wnTACntnnTATACtxnntnnT001jcos)2(Sa2e)(0www-+取A=1, T=2, =1, w0=)cos()2(Sa5 . 0)(1tnntxnT+第53頁/共100頁t=-2:0.001:2;N=input(Number of harmonics= );c0=0.5;xN=c0*ones(1,length(t); %dc componentfor n=0:2:N % even harmonics are zero xN=xN+cos(pi*n*t)*sinc(n/2);endplot(t,xN);第54頁/

19、共100頁-2-1.5-1-0.500.511.52-0.60.811.2-2-1.5-1-0.500.511.52-0.60.811.2-2-1.5-1-0.500.511.52-0.60.811.2-2-1.5-1-0.500.511.52-0.60.811.2N=5N=15N=50N=500第55頁/共100頁數值函數積分quad8可用來計算非周期信號頻譜F 是一個字符串，它表示被積函數的文件名。a,b 分別表示定積分的下限和上限quad8的返回是用自適應Simpson算法得出的積分值y = quad8(F,

20、a,b) 第56頁/共100頁)()1 ()(21tpttx-解解: 圖示三角波可表示為三角波信號頻譜的理論值為第57頁/共100頁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);X=zeros(1,N);for k=1:N X(k)=quad8(sf1,-1,1,w(k);endfigure(1);plot(w,real(X);title()xlabel(omega);ylabel(X(jomega);figure(2);plot(w,r

21、eal(X)-sinc(w/2/pi).2);xlabel(omega);title(計算誤差);第58頁/共100頁運行結果運行結果-20-15-10-505101520-20246x 10-9w計算誤差第59頁/共100頁函數fft可用來計算離散周期信號頻譜x 是離散周期信號0N-1 一個周期的序列值Xm = fft(x) Xm 是離散周期信號頻譜在0N-1 的值函數fft還可用來計算離散非周期信號頻譜、連續(xù)周期信號和連續(xù)非周期信號的頻譜。第60頁/共100頁xkk0M-M%Program 4_4 計算離散周期矩形序列的頻譜N=32; M=4; %定義周期矩形序列的參數x=ones(1,M

22、+1) zeros(1,N-2*M-1) ones(1,M); %產生周期矩形序列X=fft(x); %計算DFS系數m=0:N-1; stem(m,real(X); %畫出頻譜X的實部title(Xm的實部);xlabel(m)figurestem(m,imag(X); %畫出頻譜X的虛部title(Xm的虛部);xlabel(m)第61頁/共100頁第62頁/共100頁第63頁/共100頁) 1()j)(2()j)(1 () 1()j)(2()j)(1 ()()()j (11+-NaaaMbbbABHNNMMwwwwwwwb 分子多項式系數 a 分母多項式系數 w 需計算的H(jw)的抽樣

23、點 (數組w中少需包含兩個w的抽樣點)。計算頻響的MATLAB函數H=freqs(b,a,w)第64頁/共100頁例例1 三階歸一化的Butterworth低通濾波器的系統(tǒng)函數為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);1)j (2)j (2)j (1)j (23+wwwwH 試畫出|H(jw)| 和(w)。第65頁/共100頁三階Butterworth低通濾波器的幅度響應和相位響應第66頁/共100頁)(cosj11)

24、2(Sa21001www+tnRCnnnCRCHwwwj/1j/1)j (+RCwj11+tnnnnHCty0j0)e(j)(ww-2Sa0wnTACn2Sa21n第67頁/共100頁%p5_2 Periodic signal pass LTI systemT=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: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*c

25、os(w0*n*t+phi)*sinc(n*0.5);endplot(t,xN);xlabel(time RC=,num2str(RC);grid;set(gca,xtick,-5 -3 -1 0 1 3 5);第68頁/共100頁-5-3-113500.51RC=1t(秒 )y(t)-5-3-113500.51RC=0.1y(t)t(秒 )第69頁/共100頁NNMMaaabbbABHjj10jj10jjjeeee)e ()e ()e (-+計算頻率響應的MATLAB函數b 分子的系數 a 分母系數w 抽樣的頻率點(至少2點)， w在02之間幅度響應: abs, 相位響應: angleh =

26、 freqz(b,a,w)第70頁/共100頁b=1;a1=1 -0.9; 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);的幅度響應曲線畫例：jje11)e (-H第71頁/共100頁的幅度響應曲線畫例：jje11)e (-H第72頁/共100頁第73頁/共100頁r,p,k=residue(num,den) num,den分別為X(s)分子多項式和分母多項式的系數向量。 r為部分分式

27、的系數，p為極點，k為多項式的系數。若為真分式，則k為零。第74頁/共100頁計算多項式根roots的函數可用于計算H(s)的零極點。r=roots(D) %計算多項式D(s)的根 H(s)零極點分布圖可用pzmap函數畫出，調用形式為pzmap(sys)表示畫出sys所描述系統(tǒng)的零極點圖。第75頁/共100頁用部分分式展開法求X(s)的反變換。%program6_1format rat %將結果數據以分數的形式輸出num=1 2; den=1 4 3 0; r,p=residue(num,den) sssssX342)(23+r = -1/6 ,-1/2 ,2/3 p = -3 ,-1 ,0

28、36/115 . 03/2)(+-+-+ssssX)(e61)(e21)(32)()(31tututusXLtxtt-第76頁/共100頁用部分分式展開法求X(s)的反變換。% program6_2num=2 3 0 5;den=conv(1 1,1 1 2); %將因子相乘的形式轉換成多項式的形式r,p,k=residue(num,den)magr=abs(r) %求r的模angr=angle(r) %求r的相角)2(1532)(223+ssssssX）（第77頁/共100頁用部分分式展開法求X(s)的反變換。)2(1532)(223+ssssssX）（lr =-2.0000 + 1.133

29、9i, -2.0000 -lp =-0.5000 + 1.3229i, -0.5000 - 1.3229i, -lk =2lmagr =2.299, 2.2991, 3.0000 langr =2.6258, -2.6258, 0 3229. 1 j5 . 0e2991. 23229. 1 j5 . 0e2991. 2132)(6258. 2j6258. 2j-+-ssssX)()6258. 23229. 1cos(e1495. 1)(e3)(2)(5 . 0tuttuttxtt+-第78頁/共100頁例3 試畫出系統(tǒng) 的零極點分布圖，求其單位沖激響應h(t)和頻率響應H(jw)，并判斷系統(tǒng)是

30、否穩(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(w,abs(H)xlabel(omega)title(Magnitude Respone)1221)(23+ssssH第79頁/共100頁-1.4-1.2-1-0.8-0.6-0.4-0.20-1-0.8-0.6-0.4-0.200

31、.1Real AxisImag AxisPole-zero map012345678910-0.0500.050.40.45Impulse Responetime (s)01234567891001FrequencywMagnitude Respone第80頁/共100頁第81頁/共100頁r,p,k=residuez(num,den) num,den分別為X(z)分子多項式和分母多項式的系數向量。 r為部分分式的系數，p為極點，k為多項式的系數。若為真分式，則k為零。第82頁/

32、共100頁利用tf2zp函數計算H(z)的零極點，調用形式為z,p,k=tf2zp(b,a) H(z)零極點分布圖可用zplane函數畫出，調用形式為zplane(b,a)b和a分別為H(z)分子多項式和分母多項式的系數向量。返回值z為零點、p為極點、 k為增益常數。第83頁/共100頁例例1 1 將X(z)用部分分式展開。%program7_1num = 18; den = 18 3 -4 -1;r,p,k = residuez(num,den)運行結果為r =0.3600 , 0.2400 , 0.4000 p =0.5000 , -0.3333 , -0.3333 k =故X(z)可展開

33、為321431818)(-+zzzzX2111)3333. 01 (4 . 03333. 0124. 05 . 0136. 0)(-+-zzzzX第84頁/共100頁例2 試畫出系統(tǒng) 的零極點分布圖，求其單位沖激響應hk和頻率響應H(ej) 。% program 7_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 Respone)H,w=freqz(num,d

34、en);figure(3);plot(w/pi,abs(H)xlabel(Frequency omega)title(Magnitude Respone)3213213 . 0005. 05 . 012)(-+-+zzzzzzzH第85頁/共100頁2Real part-1-0.500.51Imaginary part-1-0.500.5105101520253035-1-0.500.511.522.5kImpulse Respone01012345678FrequencyMagnitude Respone第86頁/共100頁第87頁/共1

35、00頁 num,den 分別表示系統(tǒng)函數H(s)的分子和分母多項式； A,B,C,D 分別為狀態(tài)方程的矩陣。 A,B,C,D= tf2ss(num,den)第88頁/共100頁A, B, C, D 分別表示狀態(tài)方程的矩陣。K 表示函數ss2tf計算的與第k個輸入相關的系統(tǒng)函數，即H(s)的第k列。num 表示H(s)第k列的m個元素的分子多項式den 表示H(s)公共的分母多項式。num,den= ss2tf (A,B,C,D,k)第89頁/共100頁sys 由函數ss構造的狀態(tài)方程模型t 需計算的輸出樣本點, t=0:dt:Tfinalx(:,k) 系統(tǒng)第k個輸入在t上的抽樣值q0 系統(tǒng)的初

36、始狀態(tài)(可缺省)y(:,k) 系統(tǒng)的第k個輸出 to 實際計算時所用的樣本點； q 系統(tǒng)的狀態(tài) 獲得連續(xù)系統(tǒng)狀態(tài)方程的計算機表示模型sys = ss(A,B,C,D)求解狀態(tài)方程y,to,q=lsim(sys,x,t,q0)第90頁/共100頁sys 由函數ss構造的狀態(tài)方程模型x(:,k)系統(tǒng)第k個輸入序列q0 系統(tǒng)的初始狀態(tài)(可缺省)y(:,k)系統(tǒng)的第k個輸出n 序列的下標; q 系統(tǒng)的狀態(tài)獲得離散狀態(tài)方程的計算機表示模型sys = ss(A,B,C,D, )求解狀態(tài)方程y,n,q=lsim(sys,x, ,q0)或直接利用y, q=dlsim(A,B,C,D,x,q0)第91頁/共100頁由A,B,C,D=tf2ss(1,1 5 10)可得)()(10)( 5)( txtytyty+-01105A01B所以系統(tǒng)的狀態(tài)方程為)(01)()(01105)()(2121