1、1MATLAB語言2本章目標本章目標 熟悉熟悉MATLABMATLAB軟件平臺、工具箱、高效的數值軟件平臺、工具箱、高效的數值計算及符號計算功能；計算及符號計算功能； 熟悉熟悉MATLABMATLAB軟件的信號處理編程方法和結果軟件的信號處理編程方法和結果的可視化；的可視化； 了解數字信號處理的計算機仿真方法；了解數字信號處理的計算機仿真方法； 進一步加深對信號與系統的基本原理、方法進一步加深對信號與系統的基本原理、方法及應用的理解。及應用的理解。MATLAB語言3主要內容主要內容 7.1 7.1 基本信號的表示及可視化基本信號的表示及可視化 7.2 7.2 連續信號的時域運算與時域變換連續信

2、號的時域運算與時域變換 7.3 7.3 線性系統的時域分析線性系統的時域分析 7.4 7.4 連續時間信號的頻域分析連續時間信號的頻域分析 7.5 7.5 連續系統的復頻域分析連續系統的復頻域分析 7.6 7.6 信號采樣與重構信號采樣與重構MATLAB語言47.17.1基本信號的表示及可視化基本信號的表示及可視化 利用利用MatlabMatlab軟件信號處理工具箱軟件信號處理工具箱(Signal Processing (Signal Processing Toolbox)Toolbox)中的專用函數產生信號并繪出波形。中的專用函數產生信號并繪出波形。產生正弦波產生正弦波t=0:0.01:3t

3、=0:0.01:3* *pi;pi;y=sin(2y=sin(2* *t);t);plot(t,y)plot(t,y)產生矩形脈沖信號產生矩形脈沖信號t=-3:0.01:3;t=-3:0.01:3;y=rectpuls(t-1,2);%y=rectpuls(t-1,2);%以以t=1t=1為中心，寬度為為中心，寬度為2 2plot(t,y)plot(t,y)MATLAB語言5產生周期方波產生周期方波t=0:0.01:1;t=0:0.01:1;y=square(4y=square(4* *pipi* *t);t);plot(t,yplot(t,y) )產生周期鋸齒波產生周期鋸齒波t=0:0.00

4、1:2.5;t=0:0.001:2.5;y=sawtooth(2y=sawtooth(2* *pipi* *3030* *t);t);plot(t,yplot(t,y) )axis(0 0.2 -1 1)axis(0 0.2 -1 1)產生指數函數波形產生指數函數波形x=linspace(0,1);x=linspace(0,1);%x=0:0.01:1;%x=0:0.01:1;y=exp(-x);y=exp(-x);plot(x,yplot(x,y) )MATLAB語言6產生抽樣信號產生抽樣信號t=-10:1.5:10; t=-10:1.5:10; % %向量向量t t時間范圍時間范圍t=t1

5、:p:t2,pt=t1:p:t2,p為時間間隔為時間間隔f=f=sinc(tsinc(t); ); % %sinc(tsinc(t)=)=sin(pisin(pi* *t)/pit)/pi* *t tplot(t,fplot(t,f); %); %顯示該信號的時域波形顯示該信號的時域波形title(f(t)=Sa(t);title(f(t)=Sa(t);xlabel(txlabel(t)axis(-10,10,-0.4,1.1)axis(-10,10,-0.4,1.1)產生單位階躍信號產生單位階躍信號 定義階躍函數定義階躍函數 function f=function f=heaviside(t

6、heaviside(t) )f=(t0);f=(t0); 調用階躍函數調用階躍函數t=-1:0.01:3;t=-1:0.01:3;f=f=heaviside(theaviside(t););plot(t,fplot(t,f););axis(-1,3,-0.2,1.2);axis(-1,3,-0.2,1.2); 1000ttt MATLAB語言7產生單位沖激信號產生單位沖激信號 定義沖激函數定義沖激函數function chongji(t1,t2,t0)function chongji(t1,t2,t0)dt=0.01;t=t1:dt:t2;dt=0.01;t=t1:dt:t2;n=length

7、(t);n=length(t);x=zeros(1,n);x=zeros(1,n);x(1,(-t0-t1)/dt+1)=1/dt;%x(1,(-t0-t1)/dt+1)=1/dt;%表示在表示在t=t+t0t=t+t0時刻存在沖激時刻存在沖激stairs(t,xstairs(t,x););axis(t1,t2,0,1.2/dt);title(axis(t1,t2,0,1.2/dt);title(單位沖擊信號單位沖擊信號(t(t)調用調用chongji(-1,5,0)chongji(-1,5,0)；可以試著給別的可以試著給別的t1,t2,t0.t1,t2,t0. 10,0t dtttMATLA

8、B語言87.2 7.2 連續信號的時域運算與時域變換連續信號的時域運算與時域變換 1.1.加加( (減減) )、乘運算、乘運算要求二個信號序列長度相同要求二個信號序列長度相同. .t=0:0.01:2;t=0:0.01:2;f1=exp(-3f1=exp(-3* *t);t);f2=0.2f2=0.2* *sin(4sin(4* *pipi* *t);t);f3=f1+f2;f3=f1+f2;f4=f1.f4=f1.* *f2;f2;subplot(2,2,1);plot(t,f1);title(f1(t);subplot(2,2,1);plot(t,f1);title(f1(t);subpl

9、ot(2,2,2);plot(t,f2);title(f2(t);subplot(2,2,2);plot(t,f2);title(f2(t);subplot(2,2,3);plot(t,f3);title(f1+f2);subplot(2,2,3);plot(t,f3);title(f1+f2);subplot(2,2,4);plot(t,f4);title(f1subplot(2,2,4);plot(t,f4);title(f1* *f2);f2);MATLAB語言92.2.用用matlabmatlab的符號函數實現信號的反折、移位、的符號函數實現信號的反折、移位、尺度變換尺度變換. .由由

10、f(tf(t) )到到f(-at+b)(af(-at+b)(a0)0)步驟步驟: : f(t)f(tb)f(atb)f( atb) 移位尺度反折MATLAB語言10例例: :已知已知f(tf(t)=)=sin(t)/tsin(t)/t, ,試通過反折、移位、尺度變換由試通過反折、移位、尺度變換由f(tf(t) )得到得到f(-2t+3) f(-2t+3) 的波形的波形. .symssyms t; t;f=f=sin(t)/tsin(t)/t; %; %定義符號函數定義符號函數f(tf(t)=)=sin(t)/tsin(t)/tf1=subs(f,t,t+3);f1=subs(f,t,t+3);

11、 % %對對f f進行移位進行移位f2=subs(f1,t,2f2=subs(f1,t,2* *t);t); % %對對f1f1進行尺度變換進行尺度變換f3=subs(f2,t,-t);f3=subs(f2,t,-t); % %對對f2f2進行反褶進行反褶subplot(2,2,1);ezplot(f,-8,8);grid on;subplot(2,2,1);ezplot(f,-8,8);grid on;% % ezplotezplot是符號函數繪圖命令是符號函數繪圖命令subplot(2,2,2);ezplot(f1,-8,8);grid on;subplot(2,2,2);ezplot(f

12、1,-8,8);grid on;subplot(2,2,3);ezplot(f2,-8,8);grid on;subplot(2,2,3);ezplot(f2,-8,8);grid on;subplot(2,2,4);ezplot(f3,-8,8);grid on;subplot(2,2,4);ezplot(f3,-8,8);grid on;注注: :也可用一條指令也可用一條指令:subs(f,t,-2:subs(f,t,-2* *t+3)t+3)實現實現f(tf(t) )到到f(-2t+3)f(-2t+3)的變換的變換MATLAB語言11例：已知信號例：已知信號用用matlabmatlab求

13、求f(t+2),f(t-2),f(-t),f(2t),-f(t),f(t+2),f(t-2),f(-t),f(2t),-f(t),并繪出時域波并繪出時域波形。形。symssyms t tf=sym(t/2+1)f=sym(t/2+1)* *(heaviside(t+2)-heaviside(t-2);(heaviside(t+2)-heaviside(t-2);subplot(2,3,1);ezplot(f,-3,3);subplot(2,3,1);ezplot(f,-3,3);y1=subs(f,t,t+2);y1=subs(f,t,t+2); subplot(2,3,2),ezplot(y

14、1,-5,1); subplot(2,3,2),ezplot(y1,-5,1);y2=subs(f,t,t-2);y2=subs(f,t,t-2); subplot(2,3,3),ezplot(y2,-1,5); subplot(2,3,3),ezplot(y2,-1,5);y3=y3=subs(f,t,-tsubs(f,t,-t);); subplot(2,3,4),ezplot(y3,-3,3); subplot(2,3,4),ezplot(y3,-3,3);y4=subs(f,t,2y4=subs(f,t,2* *t);t); subplot(2,3,5),ezplot(y4,-2,2)

15、; subplot(2,3,5),ezplot(y4,-2,2);y5=-f;y5=-f; subplot(2,3,6),ezplot(y5,-3,3) subplot(2,3,6),ezplot(y5,-3,3) 1222tf tttMATLAB語言123. 3. 卷積運算卷積運算Y=Y=conv(x,hconv(x,h) )實現實現x,hx,h二個序列的卷積二個序列的卷積. .例例: :求二個信號的卷積求二個信號的卷積. .t11=0; t12=1;t11=0; t12=1;t21=0; t22=2;t21=0; t22=2;t1=t11:0.001:t12;ft1=2.t1=t11:0.

16、001:t12;ft1=2.* *rectpuls(rectpuls(t1-0.5,1);t1-0.5,1);subplot(2,2,1);plot(t1,ft1);axis(0 3 0 4)subplot(2,2,1);plot(t1,ft1);axis(0 3 0 4)t2=t21:0.001:t22;ft2=t2;t2=t21:0.001:t22;ft2=t2;subplot(2,2,2);plot(t2,ft2);axis(0 3 0 4)subplot(2,2,2);plot(t2,ft2);axis(0 3 0 4)t3=t11+t21:0.001:t12+t22;t3=t11+t

17、21:0.001:t12+t22;ft3=conv(ft1,ft2)ft3=conv(ft1,ft2)ft3=ft3ft3=ft3* *0.0010.001subplot(2,2,3);plot(t3,ft3);axis(0 3 0 4)subplot(2,2,3);plot(t3,ft3);axis(0 3 0 4) 12,01( )0,tf telse 2,020,ttftelse MATLAB語言137.37.3線性系統的時域分析線性系統的時域分析1. 1. 脈沖響應脈沖響應impulseimpulse(num,den,T)(num,den,T)y=y=impulseimpulse(nu

18、m,den,T)(num,den,T)T:T:為等間隔的時間向量為等間隔的時間向量, ,指明要計算響應的時間點指明要計算響應的時間點. .2. 2. 階躍響應階躍響應stepstep(num,den,T)(num,den,T)y=step(y=step(num,den,T)num,den,T)T:T:為等間隔的時間向量為等間隔的時間向量, ,指明要計算響應的時間點指明要計算響應的時間點. .3)3)對任意輸入的響應對任意輸入的響應lsimlsim(num,den,U,T)(num,den,U,T)y=y=lsimlsim(num,den,U,T)(num,den,U,T)U:U:任意輸入信號任

19、意輸入信號. . T:T:為等間隔的時間向量為等間隔的時間向量, ,指明要計算響應的時間點指明要計算響應的時間點. .MATLAB語言14例：已知描述某連續系統的微分方程為例：已知描述某連續系統的微分方程為：試用試用MatMatlablab繪出該系統沖激響應和階躍響應；輸入信號繪出該系統沖激響應和階躍響應；輸入信號為為 , ,該系統零狀態響應該系統零狀態響應y y（t t）。b=1 2;b=1 2;a=1 2 1;a=1 2 1;subplot(1,3,1);subplot(1,3,1);impulseimpulse(b,a); %(b,a); %沖激響應沖激響應subplot(1,3,2);

20、subplot(1,3,2);stepstep(b,a) %(b,a) %階躍響應階躍響應p=0.5; %p=0.5; %定義取樣時間間隔定義取樣時間間隔t=0:p:5; % t=0:p:5; % 定義時間范圍定義時間范圍x=exp(-2x=exp(-2* *t);t); % %定義輸入信號定義輸入信號subplot(1,3,3);subplot(1,3,3);lsim(b,a,x,tlsim(b,a,x,t);); % %對系統的輸出信號進行仿真對系統的輸出信號進行仿真 22yty ty tftf t 2tf tetMATLAB語言15例例: :對系統對系統 分別求脈沖響應、階躍響應及對輸入

21、分別求脈沖響應、階躍響應及對輸入u(tu(t)=)=sin(tsin(t) )的響應的響應. .num=1,1;num=1,1;den=1,1.3,0.8;den=1,1.3,0.8;T=0:0.1:3;T=0:0.1:3;y1=y1=impulseimpulse(num,den,T);(num,den,T);y2=y2=stepstep(num,den,T);(num,den,T);U=sin(T);U=sin(T);y3=y3=lsimlsim(num,den,U,T);(num,den,U,T);subplot(1,3,1);plot(T,y1);title(subplot(1,3,1)

22、;plot(T,y1);title(脈沖響應脈沖響應)subplot(1,3,2);plot(T,y2);title(subplot(1,3,2);plot(T,y2);title(階躍響應階躍響應)subplot(1,3,3);plot(T,y3);title(subplot(1,3,3);plot(T,y3);title(輸入為輸入為u=sintu=sint的響應的響應) 21( )1.30.8sHsssMATLAB語言167.4 7.4 連續時間信號的頻域分析連續時間信號的頻域分析1 1連續信號的傅立葉變換連續信號的傅立葉變換信號的傅立葉變換定義為：信號的傅立葉變換定義為： ()( )j

23、 tF jf t edt ( )()j tf tF jedMATLAB語言17（1 1）符號解法）符號解法調用函數：調用函數：fourierfourier、ifourierifourier例：求函數例：求函數 的傅立葉變換。的傅立葉變換。symssyms t tF=fourier(exp(-2F=fourier(exp(-2* *abs(tabs(t)F=4/(4+w2)F=4/(4+w2)例：求例：求傅立葉逆變換傅立葉逆變換syms t wsyms t wf =ifourier(1/( 1+ w2),t)f =ifourier(1/( 1+ w2),t)f =f = 1/2 1/2* *ex

24、p(-exp(-t)t)* *Heaviside(t)+1/2Heaviside(t)+1/2* *exp(t)exp(t)* *Heaviside(-t)Heaviside(-t) 2( )tf te 21()1F jMATLAB語言18例：設例：設 ，畫出，畫出f(tf(t) )及其幅頻譜。及其幅頻譜。syms t;syms t;f=1/2f=1/2* *exp(-2exp(-2* *t)t)* *sym(Heaviside(t);sym(Heaviside(t);F=F=fourierfourier(f(f););subplot(211);ezplot(f);subplot(211);e

25、zplot(f);subplot(212);ezplot(abs(F)subplot(212);ezplot(abs(F)符號解法的局限符號解法的局限： 運算結果是符號表達式，必須使用運算結果是符號表達式，必須使用ezplotezplot來作圖。如果來作圖。如果結果中含有沖激函數等項，那么結果中含有沖激函數等項，那么ezplotezplot無法作圖。無法作圖。 返回結果有可能包含一些不能直接表達的式子，甚至會返回結果有可能包含一些不能直接表達的式子，甚至會出現出現“未被定義函數或變量未被定義函數或變量”項，更加無法作圖。項，更加無法作圖。 很多場合，盡管連續，但卻不可能表示成符號表達式。很多場

26、合，盡管連續，但卻不可能表示成符號表達式。 21( )( )2tf tet( )f tMATLAB語言19（2 2）數值解法）數值解法從定義出發：從定義出發： 用用MATLABMATLAB作數值計算，它不能計算無限區間。可以根作數值計算，它不能計算無限區間。可以根據波形情況，確定積分上下限，將據波形情況，確定積分上下限，將t t分分N N等份，用求和等份，用求和代替積分代替積分。 0()( )lim()j tj nnF jf t edtf n e 12112()( ) ( ),( ),. (),.iNNj tiij tj tj tNF jf t ef tf tf teeeMATLAB語言20

27、同樣，在頻域內也是取一系列的樣本點，然后計算出不同同樣，在頻域內也是取一系列的樣本點，然后計算出不同w w處處的的F F值，然后就得到了傅立葉變換的數值解。而不同值，然后就得到了傅立葉變換的數值解。而不同w w處的處的F F值值都用上面同一個公式計算，因此可以使用都用上面同一個公式計算，因此可以使用MATLABMATLAB的矩陣計算功的矩陣計算功能。將能。將w w設為一個行向量，則設為一個行向量，則: : 其中，其中，F F是與是與w w等長的行向量，等長的行向量，tt是列向量，是列向量，w w是行向量，是行向量，tt* *w w是一個矩陣，其行數與是一個矩陣，其行數與t t相同，列數與相同，

28、列數與w w相同。相同。 *exp(* * )*Ffj twdtMATLAB語言21例：求門信號例：求門信號的傅立葉變換。的傅立葉變換。dtdt=0.02; =0.02; taotao=2;=2;t=-t=-tao:dt:taotao:dt:tao; ;f=Heaviside(t+1)-Heaviside(t-1);f=Heaviside(t+1)-Heaviside(t-1);w1=10w1=10* *pi;N=500;k=0:N;pi;N=500;k=0:N;w=kw=k* *w1/N;w1/N;F=fF=f* *exp(-jexp(-j* *tt* *w)w)* *dt;dt;F=F=r

29、eal(Freal(F););w=-fliplr(w),w(2:501); %w=-fliplr(w),w(2:501); %雙邊譜雙邊譜F=fliplr(F),F(2:501);F=fliplr(F),F(2:501);subplot(211),plot(t,f)subplot(211),plot(t,f)subplot(212),plot(w,F)subplot(212),plot(w,F) 211( )( )01tf tg tt MATLAB語言222.2.頻響特性頻響特性系統在正弦激勵下穩態響應隨信號頻率變化的特性系統在正弦激勵下穩態響應隨信號頻率變化的特性. .| |H(jH(j )

30、|:)|:幅頻響應特性幅頻響應特性. . ( ( ):):相頻響應特性相頻響應特性( (或相移特性或相移特性).).MatlabMatlab求系統頻響特性函數求系統頻響特性函數freqsfreqs的調用格式的調用格式: :h=h=freqs(num,denfreqs(num,den, , ) ) 為等間隔的角頻率向量為等間隔的角頻率向量, ,指明要計算響應的頻率點指明要計算響應的頻率點. . )()()()(jjsejHsHjHMATLAB語言23例例: :求系統求系統的頻響特性的頻響特性. . num=1,1; num=1,1;den=1,1.3,0.8;den=1,1.3,0.8;W=0:

31、0.1:100;W=0:0.1:100;h=freqs(num,den,W);h=freqs(num,den,W);subplot(2,1,1);plot(W,abs(h);title(subplot(2,1,1);plot(W,abs(h);title(幅頻特性幅頻特性)axis(0,20,0,1.5);axis(0,20,0,1.5);set(gca,xtick,0,10,20);set(gca,xtick,0,10,20);set(gca,ytick,0,1/sqrt(2),1.25);grid on;set(gca,ytick,0,1/sqrt(2),1.25);grid on;sub

32、plot(2,1,2);plot(W,angle(h);title(subplot(2,1,2);plot(W,angle(h);title(相頻特性相頻特性)axis(0,20,-pi/2,0.2);axis(0,20,-pi/2,0.2);set(gca,xtick,0,10,20);set(gca,xtick,0,10,20);set(gca,ytick,-pi/2,-pi/4,0);grid on;set(gca,ytick,-pi/2,-pi/4,0);grid on; 21( )1.30.8sH sssMATLAB語言247.5 7.5 連續系統的復頻域分析連續系統的復頻域分析1.

33、 1. 系統函數系統函數H(sH(s) )定義定義系統零狀態響應的拉氏變換與激勵的拉氏變換之比系統零狀態響應的拉氏變換與激勵的拉氏變換之比. .H(sH(s)=)=R(s)/E(sR(s)/E(s) )在在matlabmatlab中中, , 傳遞函數描述法是通過傳遞函數分子和分母關傳遞函數描述法是通過傳遞函數分子和分母關于于s s降冪排列的多項式系數來表示的降冪排列的多項式系數來表示的. .例如例如, ,某系統傳遞函數：某系統傳遞函數：則可用如下二個向量則可用如下二個向量numnum和和denden來表示來表示: :num=1,1num=1,1den=1,1.3,0.8den=1,1.3,0.

34、8 21( )1.30.8sH sssMATLAB語言252. 2. 系統零、極點分布與系統穩定性關系系統零、極點分布與系統穩定性關系 系統函數系統函數H(sH(s) )集中表現了系統的性能集中表現了系統的性能, ,研究研究H(sH(s) )在在S S平面中極平面中極點分布的位置點分布的位置, ,可很方面地判斷系統穩定性可很方面地判斷系統穩定性. .u穩定系統穩定系統: : H(sH(s) )全部極點落于全部極點落于S S左半平面左半平面( (不包括虛軸不包括虛軸),),即若即若滿足，則系統是穩定的滿足，則系統是穩定的. .u不穩定系統不穩定系統: : H(sH(s) )極點落于極點落于S S

35、右半平面右半平面, ,或在虛軸上具有二階或在虛軸上具有二階以上極點以上極點, ,則在足夠長時間后則在足夠長時間后, ,h(th(t) )仍繼續增長仍繼續增長, , 系統是不穩定系統是不穩定u臨界穩定系統臨界穩定系統: : H(sH(s) )極點落于極點落于S S平面虛軸上平面虛軸上, ,且只有一階且只有一階, ,則在則在足夠長時間后足夠長時間后, ,h(th(t) )趨于一個非零數值或形成一個等幅振蕩趨于一個非零數值或形成一個等幅振蕩. . H(sH(s) )零、極點零、極點可用可用matlabmatlab的多項式求根函數的多項式求根函數rootsroots求得：求得：極點極點: :p=p=r

36、oots(denroots(den) )零點零點: :z=z=roots(numroots(num) )根據根據p p和和z z，用，用plotplot命令命令即可畫出系統零極點分布圖即可畫出系統零極點分布圖, ,進而分析判進而分析判斷系統穩定性。斷系統穩定性。MATLAB語言26例例: : 系統函數系統函數畫出系統零、極點分布圖畫出系統零、極點分布圖, , 判斷該系統穩定性。判斷該系統穩定性。num=1,0,-4;num=1,0,-4;den=1,2,-3,2,1;den=1,2,-3,2,1;p=p=roots(denroots(den) )z=z=roots(numroots(num)

37、)plot(real(p),imag(pplot(real(p),imag(p),),* *);hold on;);hold on;plot(real(z),imag(z),o);gridplot(real(z),imag(z),o);grid on on 由系統零極點分布圖可知由系統零極點分布圖可知, ,該系統有一對極點位于該系統有一對極點位于s s右半右半平面平面, ,故系統是不穩定的故系統是不穩定的. . 24324( )2321sH sssssMATLAB語言277.6 7.6 信號采樣與重構信號采樣與重構Nyquist 采樣定理： 對于帶限信號 f t ,如果其帶寬為mw , 那么以

38、大于 2mw 的頻率對其進行采樣，得到采樣信號 sft，則該信號 f t 就可由采樣信號 sft 完全確定，將采樣信號通過一個低通濾波器，可以恢復出原始信號 f t。 那么，若采樣間隔為sT ，應有smTw，則 ssssnnftf ttnTf nTtnT MATLAB語言28信號重構： 下面我們看如何由 sft ，來重構 f t 。 sft 經過一低通濾波器，那么 *sf tfth t 對于低通濾波器，有： 1cscwh tFHjwTSa w t 則有：( )( )()()()cscsscscsnwT wf tf tTSa w tf nT Sa w tnT 因為sin()sin ( )tc t

39、Satt，所以t( )sincSa t（ ） 這樣，上面的式子可以寫成：( )()sin ()ccsssnwwf tTf nTctnT MATLAB語言29例：信號的采樣與重構。例：信號的采樣與重構。wm=1 ; %wm=1 ; %信號帶寬信號帶寬wcwc=wm; %=wm; %濾波器截至頻率濾波器截至頻率tsts=pi/wm %=pi/wm %采樣間隔，臨界采樣采樣間隔，臨界采樣 wsws=2=2* *pi/pi/tsts; %; %采樣角頻率采樣角頻率n=-100:100; %n=-100:100; %時域采樣點數目時域采樣點數目ntsnts=n=n* *tsts; %; %時域采樣點時域采樣點f=f=sinc(ntssinc(nts/pi); %/pi); %采樣信號采樣信號dtdt=0.005;t=-15:dt:15; =