matlab在信號與系統(tǒng)中的應(yīng)用

第9章在信號與系統(tǒng)中的應(yīng)用9.1連續(xù)信號和系統(tǒng)嚴格說來，MATLAB（就基本部分而言）是不能表示連續(xù)信號的，因為它給出的是各個樣本點的數(shù)據(jù)。只是當(dāng)樣本點取得很密時可看成連續(xù)信號，什么叫密，要相對于信號變化的快慢而言，形象地說，在所有相鄰樣本點之間的數(shù)據(jù)變化必須非常小才能看成‘密’，其嚴格的數(shù)學(xué)定義此處不予討論。以下均假定相對于采樣點密度而言,信號變化足夠慢。例9-1-1

連續(xù)信號的MATLAB描述（144）列出單位脈沖、單位階躍、復(fù)指數(shù)函數(shù)等連續(xù)信號的MATLAB表達式。程序exn911(1)clear,t0=0;tf=5;dt=0.05;t1=1;t=[t0:dt:tf];%(1)單位脈沖信號,%在t1(t0≤t1≤tf)處有面積為1的脈沖信號。

t=[t0:dt:tf];st=length(t); n1=floor((t1-t0)/dt); %求t1對應(yīng)的樣本序號

x1=zeros(1,st); %把信號先初始化為零

x1(n1)=1/dt; %給出t1處的脈沖信號

subplot(2,2,1),stairs(t,x1)%繪圖，用stairs命令

axis([0,5,0,22]) %為了使脈沖頂部避開圖框程序exn911(2)%(2)單位階躍信號，%信號從t0到tf，在t1前為0，到t1處躍變?yōu)?.% 程序前幾句即求t,st,n1的語句與上同%產(chǎn)生階躍信號x2=[zeros(1,n1-1),ones(1,st-n1+1)]; subplot(2,2,3),stairs(t,x2) %繪圖axis([0,5,0,1.1]) %改變圖框坐標%(3)復(fù)數(shù)指數(shù)信號u=-0.5;w=10;x4=exp((u+j*w)*t);subplot(2,2,2),plot(t,real(x4)) %繪圖,subplot(2,2,4),plot(t,imag(x4)) %繪圖,程序exn911運行結(jié)果x1,x2,x3,x4的波形見右圖.注意若要顯示連續(xù)信號波形中的不連續(xù)點，用stairs命令；而要使波形光滑些，則用plot命令較好。復(fù)數(shù)指數(shù)信號可以分解為余弦和正弦信號，它們分別是復(fù)數(shù)信號的實部和虛部。右圖中的兩個衰減振蕩信號就代表了這兩個相位差90度的分量。例9-1-2線性系統(tǒng)沖擊響應(yīng)（145）編寫求任意高階連續(xù)常系數(shù)線性系統(tǒng)沖擊響應(yīng)的程序。解:◆建模這個問題在第四章介紹多項式函數(shù)庫時已打下了基礎(chǔ)，在第七章例7-3-1中又給出了解法，讀者可先看懂那些內(nèi)容，然后再看本題。任意階次的連續(xù)線性系統(tǒng)可用下列微分方程表述：寫成傳遞函數(shù)形式為其特性可用系統(tǒng)傳遞函數(shù)的分子分母系數(shù)向量b和a來表示。

傳遞函數(shù)反變換的求法

如果分母系數(shù)多項式?jīng)]有重根，則可以把兩個多項式之比分解成n個一階部分分式之和。即： 其中p1,p2,…,pn是分母多項式的n個根，而r1,r2,…,rn是則對應(yīng)于這n個根的留數(shù)。一階分式的反變換可以查表得到，容易寫出沖擊響應(yīng)的公式如下： 可見只要求出根向量p和留數(shù)向量r，線性方程的解就得到了。求根是代數(shù)問題，當(dāng)階次很高時，代數(shù)方程沒有解析解?？上驳氖荕ATLAB提供了用數(shù)值方法求根和留數(shù)的函數(shù)residue.m，它的調(diào)用格式為：

[r,p]=residue(b,a)程序exn912a=input('分母系數(shù)向量a=(書上取poly([0,-1,-2,-5])）')b=input('分子系數(shù)向量b=(書上取[1,7,1]）')[r,p]=residue(b,a), k=input('是否要求波形？是,鍵入1;否,鍵入0');ifk==1dt=input('dt=(書上取0.05)');tf=input('tf=(書上取5)');

t=0:dt:tf;h=zeros(1,length(t));

fori=1:length(a)-1 h=h+r(i)*exp(p(i)*t);

end,plot(t,h),gridelse,end程序exn912運行結(jié)果給出系統(tǒng)的傳遞函數(shù)為運行程序exn912，依次輸入：

（注意用poly函數(shù)把極點向量p=[0,-1,-2,-5]轉(zhuǎn)換成系數(shù)向量a）a=poly([0,-1,-2,-5])B=[1,7,1]，dt=0.05，tf=5得出的h(t)如右圖所示。9-1-3線性系統(tǒng)零輸入響應(yīng)的計算（146）線性時不變連續(xù)系統(tǒng)的特性可用常微分方程表示為：求其零輸入響應(yīng)。解：在零輸入條件u=0時，等式右端為零。系統(tǒng)響應(yīng)的通解為：

其中，p是特征方程的n個根組成的向量[p1,p2,…,pn]，其每個分量的系數(shù)Cn則由y及其各階導(dǎo)數(shù)的初始條件y0,Dy0,…,D(n-1)y0來確定。代入初始條件得到的矩陣方程初始條件數(shù)應(yīng)該和常數(shù)數(shù)相等，由此構(gòu)成一個確定C1,…,Cn的線性代數(shù)方程組，寫成:矩陣V由特征根向量p確定，這種矩陣稱為范德蒙特矩陣。在MATLAB中，有生成它的函數(shù)vander(p)。

求零輸入響應(yīng)程序exn913它產(chǎn)生的矩陣與上述矩陣排列轉(zhuǎn)了90度，故用V=rot90(vander(p))，按此思路編成程序exn913:a=input('輸入分母系數(shù)向量a=[a1,a2,...]=');n=length(a)-1;Y0=input('初始條件

Y0=[y0,Dy0,D2y0,...]=');p=roots(a);V=rot90(vander(p))；%生成系數(shù)矩陣Vc=V\Y0‘; %用左除求出系數(shù)向量C%以下是計算并畫出時間響應(yīng)的程序段dt=input('dt=');tf=input('tf=')%生成t向量和y初始向量t=0:dt:tf;y=zeros(1,length(t));fork=1:ny=y+c(k)*exp(p(k)*t);endplot(t,y),grid %繪圖數(shù)字實例用這一個普遍程序來解一個三階系統(tǒng)，設(shè)其微分方程為：

初始條件為： ，求零輸入響應(yīng)。解：運行exn913，按提示輸入a為[1,2,9,3];Y0為[1,0,0];dt為0.1;tf為5；即可得到t=0~5秒的零輸入響應(yīng)曲線。

再分別取Y0為[0,1,0];[0,0,1]，用holdon語句使三次運行生成的圖形疊畫在一幅圖上，得到右圖。

例9-1-4重極點求反變換（147）命題：n級放大器，每級的轉(zhuǎn)移函數(shù)均為，求階躍輸入下的過渡過程，畫出n不同時的波形及頻率特性。解：◆建模系統(tǒng)的轉(zhuǎn)移函數(shù)為H(s)=，階躍輸入的拉普拉斯變換為1/s，因此輸出為兩者的乘積，即求Y(s)的拉普拉斯反變換，即可得到輸出過渡過程y(t)。這里我們遇到了一個有多重極點-wn的H(s)求拉普拉斯反變換的問題，數(shù)學(xué)上比較麻煩。為了避開重極點問題，可以有意把極點拉開一些，例如設(shè)n個極點散布在-0.98wn到1.02wn之間，由此就可用下面的程序來求：程序exn914(1)clear,clf,N=input('輸入放大器級數(shù)N='); wn=1000;dt=1e-4;tf=0.01;t=0:dt:tf;y=zeros(N,length(t)); %輸出初始化forn=1:Np0=-linspace(.95,1.05,n)*wn; %將H(s)極點分散布置

ay=poly([p0,0]); %由Y(s)的極點求分母系數(shù)

by=prod(abs(p0)); %求Y(s)的分子系數(shù)

[r,p]=residue(by,ay); %求Y(s)的留數(shù)極點

fork=1:n+1 %把各時域分量相加

y(n,:)=y(n,:)+r(k)*exp(p(k)*t); endfigure(1),plot(t,y(n,:));grid,holdon %繪制過渡過程end程序exn914(2)-畫波德圖%下面的語句用來繪制波德圖，如果用bode函數(shù),只要一句% figure(2),bode(prod(abs(p0)),poly(p0));holdon bh=by;ah=poly(p0); %求H(s)的分子分母系數(shù)

w=logspace(2,4); %給出頻率范圍和分度

H=polyval(bh,j*w)./polyval(ah,j*w);%求H(jw)aH=unwrap(angle(H))*180/pi; %求出相角

fH=20*log10(abs(H)); %求出振幅

figure(2),subplot(2,1,1),semilogx(w,fH),gridon,holdon %繪幅頻圖

subplot(2,1,2),semilogx(w,aH),gridon,holdon%相頻圖end程序exn914運行結(jié)果(時域)運行此程序，設(shè)N=4,可得過渡過程如右圖,從中看出輸出信號達到0.6處所需的時間約為單級時常數(shù)乘以級數(shù)。由于極點互相接近，此程序在N>4時又會出現(xiàn)很大誤差。程序exn914運行結(jié)果(頻域)右圖繪出了多級放大器的頻率特性，其幅特性（圖上為分貝數(shù)），顯示了低通特性，隨級數(shù)的增加,通帶減小,從相特性看出，隨級數(shù)的增加,負相移成比例地增加。例9-1-5方波通過濾波器（148）設(shè)方波信號的寬度為5秒，信號持續(xù)期為10秒，試求其在0～20【1/秒】頻段間的頻譜特性。如只取從0～10【1/秒】的頻譜分量作反變換（相當(dāng)于通過了一個低通濾波器），求其輸出波形。解：◆建模設(shè)信號的時域波形f(t),在0到10秒的區(qū)間外信號為零,則其付利葉變換為:按MATLAB作數(shù)值計算的要求,必須把t分成N份,用相加來代替積分,對于任一ω,可寫成:積分轉(zhuǎn)化為求和運算這說明求和的問題可以用f(t)行向量乘以ejωt列向量來實現(xiàn).此處的Δt是t的增量,在程序中,將用dt來代替.由于要求出一系列不同的ω處的F值,而都可用同一公式,這就可以利用MATLAB中的元素群運算能力,把ω設(shè)成一個行數(shù)組,分別代入本公式左右端的ω中去,寫成: F=f*exp(-j*t’*w)*Δt其中,F是與ω同長的行向量,exp中的t’是列向量,w是行向量,t’*w是一個矩陣,其行數(shù)與t相同,列數(shù)與w相同.這個式子就完成了付利葉變換,類比地可以得出付利葉逆變換表示式.由此得到下面的付利葉變換程序。程序exn915clear,tf=10;N=256;t=linspace(0,tf,N); %給出時間分割w1=linspace(eps,20,N);dw=20/(N-1);%dw=1/4/tf;w1=[eps:dw:(N-1)/4/tf]; %給出頻率分割f=[ones(1,N/2),zeros(1,N/2)]; %給出信號（此處是方波）F1=f*exp(-j*t'*w1)*tf/(N-1); %求付利葉變換w=[-fliplr(w1),w1(2:N)]; %補上負頻率F=[fliplr(F1),F1(2:N)]; %補上負頻率區(qū)的頻譜w2=w(N/2:3*N/2); %取出中段頻率F2=F(N/2:3*N/2); %取出中段頻譜subplot(1,2,1),plot(w,abs(F),'linewidth',1.5),gridf1=F2*exp(j*w2'*t)/pi*dw; %對中段頻譜求付利葉逆變換subplot(1,2,2),plot(t,f,t,f1,'linewidth',1.5),grid程序exn915運行結(jié)果◆執(zhí)行這個程序的結(jié)果見下圖,因為方波含有很豐富的高頻分量，要充分恢復(fù)其原來波形需要很寬的頻帶，實踐中不太可能做到。9.2離散信號和系統(tǒng)信號可以粗略地分為模擬信號和數(shù)字信號.模擬信號將用x(t)表示,其中變量t可以表示任何物理量,但我們假定它代表以秒為單位的時間.離散信號用x(n)表示.其中變量n為整數(shù)并代表時間的離散時刻.因此它也稱為離散時間信號.他是一個數(shù)字的序列并可用以下符號之一來表述:x(n)={x(n)}={...,x(-1),x(0),x(1),...} ↑其中,向上的箭頭表示在n=0處的取樣.離散信號的MATLAB表示在MATLAB中,我們可以用一個列向量來表示一個有限長度的序列.然而這樣一個向量并沒有包含采樣位置的信息.因此,完全地表示x(n)要用x和n兩個向量，例如序列x(n)={2,1,-1,5,1,4,3,7}（下面的

↑箭頭為第0個采樣點），在MATLAB中表示為:n=[-3,-2,-1,0,1,2,3,4];x=[2,1,-1,0,1,4,3,7];當(dāng)不需要采樣位置信息或這個信息是多余的時候(例如該序列從n=0開始),我們可以只用x向量來表示.由于有限的內(nèi)存,MATLAB無法表示無限序列.例9-2-1基本脈沖序列程序（149）1.單位脈沖序列起點n0,終點nf,在ns處有一單位脈沖(n0≤ns≤nf),2.單位階躍序列起點n0,終點nf,在ns前為0,在ns后為1(n0≤ns≤nf),3.復(fù)數(shù)指數(shù)序列:4.復(fù)數(shù)指數(shù)序列:

例9-2-1的程序解:◆建模 這些基本序列的表達式比較簡明,編寫程序也不難,對單位脈沖序列,此處提供了兩種方法,其中用邏輯關(guān)系的編法比較簡潔.讀者可從中看到MATLAB編程的靈活性和技巧性,要多用多想才能編出高明簡潔的程序。繪制脈沖序列,通常用stem語句.◆MATLAB程序clear,no=0;nf=10;ns=3;n1=n0:nf;x1=[zeros(1,ns-n0),1,zeros(1,nf-ns)]; %n1=n0:nf;x1=[(n1-ns)==0]; %顯然,用邏輯式是比較高明的方法n2=n0:nf;x2=[zeros(1,ns-n0),ones(1,nf-ns+1)]; %也有類似的用邏輯比較語句的方法,留給讀者思考例9-2-1的程序(續(xù))n3=n0:nf;x3=(0.9).^n3; %實數(shù)指數(shù)序列n4=n0:nf;x4=exp((-0.2+0.3j)*n3); %復(fù)數(shù)指數(shù)序列subplot(2,2,1),stem(n1,x1);subplot(2,2,2),stem(n2,x2);subplot(2,2,3),stem(n3,x3);subplot(4,2,6),stem(n4,real(x4)); %注意sunplot的輸入變元subplot(4,2,8),stem(n4,imag(x4));line([0,10],[0,0]), %畫橫坐標例9-2-1的程序運行結(jié)果例9-2-2離散付利葉變換的計算（150）解：◆建模一個時間序列x(n)的離散時間付利葉變換的定義為：如果序列的長度是有限的，可以把它看作是周期性無限序列中的一個周期，其長度為N,對這個周期性序列可以用離散付利葉變換（注意少了‘時間’兩字）進行研究。它的定義為：其中

離散傅里葉變換化為矩陣乘積用例9-1-5中的方法,引入矩陣乘法來實現(xiàn)求和運算,用元素群算法來求不同k時的X,把n和k都設(shè)成1×N的行數(shù)組，則nk=n’*k(以及WN.^nk)成為一個N×N的方陣則有 X=x*WN.^nk

即MATLAB只能處理有限長度的序列，因此，適合于計算離散付利葉變換及其逆變換。離散傅里葉變換程序exn922設(shè)有限信號序列xn(n)的長度為Nx，則按定義，其N點付利葉變換Xk(k)的的程序為：xn=input(‘x=‘);N=(length(xn)); %取N為Nxn=[0:1:N-1];k=[0:1:N-1];%設(shè)定n和k行向量WN=exp(-j*2*pi/N);%Wn因子nk=n'*k;%產(chǎn)生一個含nk值的N乘N維矩陣WNnk=WN.^nk;%DFT矩陣Xk=xn*WNnk;%DFT系數(shù)的行向量plot(abs(Xk)),grid %繪幅頻特性圖MATLAB中的fft子程序這個程序的運算速度比較低，實際上MATLAB已提供了快速離散付利葉變換的函數(shù)fft可直接調(diào)用。其調(diào)用格式為：

X=fft(x,N)x是輸入的時間序列，N則是付利葉變換取的點數(shù)，若省略N，則它自動把x的長度作為N。當(dāng)N取2的冪時，變換速度最快，所以要提高fft函數(shù)的運行速度，程序應(yīng)編寫如下：xn=input(‘x=‘);N=pow2(nextpow2((length(xn)))); %取N為大于Nx而最接近的2的冪tic,X=fft(xn,N);toc要注意X是一個長度為N的復(fù)數(shù)數(shù)組，可以分解出它的振幅和相位頻率分別繪圖。程序exn922運行結(jié)果令Nx=700;n=1:Nx;x=sin(0.1*n)+randn(1,Nx);然后調(diào)用上述程序之一,輸入此x,所得其fft的幅度特性如右圖所示。在第一個程序的末兩行之前加上tic，在其后加上toc,則屏幕上會運行這兩語句句所需的時間，在作者的計算機上約為8秒；而同樣的fft只需0.05秒。計算逆離散付利葉變換的程序與此相仿，留待讀者自行編寫。信號的fft

的振幅頻率特性9.3系統(tǒng)函數(shù)的計算機求法例9-3-1簡單信號流圖模型的矩陣解法[11]信號流圖是用來表示和分析復(fù)雜系統(tǒng)內(nèi)的信號變換關(guān)系的工具。其基本概念如下：（1）系統(tǒng)中每個信號用圖上的一個節(jié)點表示。如圖中的u,x1,x2。（一般物流圖中是把物流標在箭桿上的。）（2）系統(tǒng)部件對信號實施的變換關(guān)系用有向線段表示，箭尾為輸入信號，箭頭為輸出信號，箭身標注對此信號進行變換的乘子。如圖上的G1，G2。如果乘子為1，可以不必標注。（3）每個節(jié)點信號的值等于所有指向此節(jié)點的箭頭信號之和，每個節(jié)點信號可以向外輸出給多個部件，其值不變。簡單信號流圖的數(shù)學(xué)模型根據(jù)這幾個概念，可以列出右圖的方程如下。寫成矩陣方程或 x=QxPu帶反饋的簡單信號流圖矩陣方程的化簡(1)移項整理，可以得到信號向量x的公式(I–Q)x=Pux=inv(I–Q)*Pu

定義系統(tǒng)的傳遞函數(shù)W為輸出信號與輸入信號之比x/u，則W可按下式求得：W=x/u=inv(I–Q)*P矩陣方程的化簡(2)因為求二階矩陣的逆可以直接用下面的公式：

所以即對x1的傳遞函數(shù)為，對x2的傳遞函數(shù)為。用MATLAB解信號流圖對于階次高的情況，求逆就必須用軟件工具了。如果信號流圖中有G1那樣的符號變量，那么它的求解要用符號運算工具箱，對于本題，其MATLAB程序exn931為：symsG1G2Q=[0,-G2;G1,0],P=[1;0]W=inv(eye(2)Q)*P程序運行的結(jié)果是 與前面的結(jié)果相同。用經(jīng)典的“梅森公式”求信號流圖的解非常繁瑣，而且無法機械化。用矩陣代數(shù)方法的最大好處是可以向任意高的階次、任意復(fù)雜的信號流圖推廣，實現(xiàn)復(fù)雜系統(tǒng)傳遞函數(shù)推導(dǎo)的自動化。到現(xiàn)在為止，我們還沒有見過任何一本書籍用矩陣方法來推導(dǎo)這個公式。例9-3-2

較復(fù)雜的信號流圖圖9-11就是一個較復(fù)雜些的信號流圖。照上述方法列出它的方程如下：x1=G4x3ux2=G1x1G5x4x3=G2x2x4=G3x3列為矩陣方程，得到帶雙重反饋的信號流圖程序exn932公式W=x/u=inv(I–Q)*P同樣是正確的，不過這里的Q和P分別為44和41矩陣，用手工求逆是辦不到了。我們可采用類似的MATLAB程序exn931：symsG1G2G3G4G5Q=[0,0,-G4,0;G1,0,0,-G5;0,G2,0,0;0,0,G3,0],P=[1;0;0;0]W=inv(eye(4)-Q)*PPretty(W(4))程序運行的結(jié)果為：exn932的運行結(jié)果當(dāng)系統(tǒng)內(nèi)各個環(huán)節(jié)都是線性集總參數(shù)，因而它們的傳遞函數(shù)Gi都可以表為s的多項式有理分式時，不管其階次有多高，傳遞函數(shù)W可以很容易由計算機直接自動算出。這個方法還可以推廣到離散系統(tǒng)，用來計算任意復(fù)雜的數(shù)字濾波器系統(tǒng)函數(shù)。9.4頻譜及其幾何意義頻譜分析是信號與系統(tǒng)課程中最重要的內(nèi)容之一，許多讀者在學(xué)習(xí)中感到抽象，往往只能從數(shù)學(xué)上承認時域信號與它的頻譜之間的變換關(guān)系，而沒有理解它的物理意義。用MATLAB可以幫助讀者建立形象的幾何概念，真正掌握它。首先來看歐拉公式，它是以最簡明的方式建立了信號頻域與時域的關(guān)系：它說明一個最簡單的實余弦信號可以由正、負兩個Ω0頻率分量合成。在復(fù)平面上，正的Ω0對應(yīng)于反時針旋轉(zhuǎn)的向量，負的Ω0對應(yīng)于順時針旋轉(zhuǎn)的向量，當(dāng)這兩個向量的幅度相同，而相角符號相反時，就合成為一個在實軸上的向量。它的相角為零，大小按正弦變化，形成了實信號cosΩ0t。實周期信號是頻譜向量的合成推而廣之，任何實周期信號必然具有正、負兩組頻頻率的頻譜成分，正、負頻率頻譜的幅度對稱而相位反對稱，或者說，是共軛的。如果頻譜不止這兩項，而是有四項或更多，它們的合成仍然可以用幾何動畫來表示?？梢园衙總€頻譜看作一根長度等于頻譜幅度、按頻率Ω旋轉(zhuǎn)的桿件，頻譜的相加等價于多節(jié)桿件首尾相接，桿件末端的軌跡就描述了生成的時域波形。因為這個端點是在平面上運動，所以它將產(chǎn)生復(fù)信號，在實軸和虛軸上的投影分別為實信號和虛信號。多個頻率向量的合成程序【例9-4-1】設(shè)計一個演示程序，它能把四個用戶任意給定集總頻譜合成并生成對應(yīng)的時域信號。解：建模按上述多節(jié)桿合成模型程序設(shè)計包括三個主要部分：（1）各頻譜分量的輸入，包括其幅度和頻率（有正負號）；（2）將各分量當(dāng)作轉(zhuǎn)動的桿件首尾相接；（3）記錄多節(jié)桿系末端的軌跡畫出圖形。程序exn941%（1）給頻譜向量賦值N=input('N(輸入向量個數(shù)，限定N不大于4)='); fori=1:N