《MATLAB及其在理工課程中的應(yīng)用指南》課件第9章

第9章MATLAB在信號和系統(tǒng)中的應(yīng)用舉例

9.1連續(xù)信號和系統(tǒng)

9.2離散信號和系統(tǒng)

9.3系統(tǒng)函數(shù)

9.4頻譜及其幾何意義

9.1連續(xù)信號和系統(tǒng)

【例9-1-1】列出單位脈沖、單位階躍、復(fù)指數(shù)函數(shù)等連續(xù)信號的MATLAB表達(dá)式。解：◆建模嚴(yán)格說來，MATLAB是不能表示連續(xù)信號的，因為它給出的是各個樣本點的數(shù)據(jù)，只有當(dāng)樣本點取得很密時才可看成連續(xù)信號。所謂“密”，是相對于信號變化的快慢而言，形象地說，在相鄰樣本點之間的數(shù)據(jù)變化必須非常小才能看成“密”，其嚴(yán)格的數(shù)學(xué)定義此處不予討論。在編程中，先設(shè)定共同的時間坐標(biāo)，然后分別列出生成三種信號的程序。◆MATLAB程序

clear，t0=0;tf=5;dt=0.05;t1=1;t=［t0∶dt∶tf］;

%(1)單位脈沖信號

%在t1(t0≤t1≤tf)處有一持續(xù)時間為dt，面積為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處幅度為1/dt的脈沖信號

subplot(2，2，1)，stairs(t，x1) %繪圖，注意為何用stairs而不用plot命令

axis(［0，5，0，1.1/dt］) %為了使脈沖頂部避開圖框，改變圖框坐標(biāo)

%(2)單位階躍信號

%信號從t0到tf，在t1(t0≤t1≤tf)前為0，到t1處有一躍變，以后為1

x2=［zeros(1，n1-1)，ones(1，st-n1+1)］;%產(chǎn)生階躍信號

subplot(2，2，3)，stairs(t，x2)%繪圖

axis(［0，5，0，1.1］) %為了使方波頂部避開圖框，改變圖框坐標(biāo)

%(3)復(fù)數(shù)指數(shù)信號

u=-0.5;w=10;x3=exp((u+j*w)*t);subplot(2，2，2)，plot(t，real(x3))%繪圖

subplot(2，2，4)，plot(t，imag(x3))%繪圖◆程序運行結(jié)果

x1、x2、x3的波形如圖9-1所示。注意：若要顯示連續(xù)信號波形中的不連續(xù)點，用stairs命令；若要使波形光滑些，則用plot命令較好。復(fù)數(shù)指數(shù)信號x3可以分解為余弦和正弦信號，它們分別是復(fù)數(shù)信號的實部和虛部，右圖中的兩個衰減振蕩信號就代表了這兩個相位差90°的分量。圖

9-1例9-1-1中x1、

x2、

x3對應(yīng)的四種波形

【例9-1-2】編寫求任意高階連續(xù)常系數(shù)線性系統(tǒng)沖擊響應(yīng)的程序。

解:

這個問題在第4章4.3.5節(jié)介紹多項式函數(shù)庫時已經(jīng)打下基礎(chǔ)，在第7章機(jī)械振動的例7-3-1又討論過二階常系數(shù)線性微分方程的解法，讀者可以先看懂那些例題的解法，再看本題。任意階次的連續(xù)線性系統(tǒng)可用下列線性常微分方程表述:寫成傳遞函數(shù)形式為

因此，其特性可以用系統(tǒng)傳遞函數(shù)的分子分母系數(shù)向量b和a來表示。向量b和a的長度分別為m+1和n+1，對于一切物理上可實現(xiàn)的系統(tǒng)，必有n≥m。系統(tǒng)的沖擊響應(yīng)等于傳遞函數(shù)的拉普拉斯反變換，問題歸結(jié)為如何求出這個反變換。如果分母系數(shù)多項式?jīng)]有重根，則可以將兩個多項式之比分解成n個一階部分分式之和，即

其中p1，p2，…，pn是分母多項式的n個根，而r1，r2，…，rn則是對應(yīng)于這n個根的留數(shù)。一階分式的反變換可以查表得到，所以可得到?jīng)_擊響應(yīng)的公式:可見只要求出根p和留數(shù)r，線性方程的解就得到了。求根是代數(shù)問題，當(dāng)階次高時，沒有解析解?？上驳氖荕ATLAB提供了用數(shù)值方法求根和留數(shù)的residue函數(shù)，它的調(diào)用方法如下:

［r，p］=residue(b，a)只要給出系數(shù)向量b和a，就可以得出根p和留數(shù)r，并得到系統(tǒng)的沖擊響應(yīng)。因此，本題的程序就可方便地編寫如下?！鬗ATLAB程序

a=input(′多項式分母系數(shù)向量a=′)

b=input(′多項式分子系數(shù)向量b=′)

［r，p］=residue(b，a)，%求留數(shù)

disp(′沖擊響應(yīng)的解析式為h(t)=Σr(i)*exp(p(i)*t)′)

k=input(′是否要求波形？是，鍵入1;否，鍵入0′);

ifk==1

dt=input(′dt=′);tf=input(′tf=′);%設(shè)定時間數(shù)組

t=0∶dt∶tf;h=zeros(1，length(t));%h的初始化

fori=1∶length(a)-1%根的數(shù)目等于a的長度減1

h=h+r(i)*exp(p(i)*t);%疊加各根分量

end

plot(t，h)，grid

else，

end◆程序運行結(jié)果例如，給出系統(tǒng)傳遞函數(shù)為求沖擊響應(yīng)。根據(jù)程序提問依次輸入: a=poly(［0，-1，-2，-5］) b=［1，7，1］

dt=0.05 tf=5得出的h(t)如圖9-2所示。圖

9-2高階系統(tǒng)的沖擊響應(yīng)

程序中要的是系數(shù)向量a，而題中給出的是極點向量p=［0，-1，-2，-5］，因此這里用poly函數(shù)來作轉(zhuǎn)換。

【例9-1-3】線性時不變系統(tǒng)的特性可用常系數(shù)線性微分方程表示:求輸入u為0時，由初始狀態(tài)決定的輸出，即其零輸入響應(yīng)。

解：◆建模在零輸入條件下，系統(tǒng)的響應(yīng)取決于微分方程左端特征方程的根，與右端無關(guān)，其通解為其中，p1、p2、…、pn是特征方程a1λn+a2λn-1+…+anλ+an+1=0的根。每個分量的系數(shù)C1、C2、…、Cn

應(yīng)由y及其各階導(dǎo)數(shù)的初始條件來確定。

(Dy0表示y的導(dǎo)數(shù)的初始值)初始條件數(shù)應(yīng)該和待定系數(shù)的數(shù)目相等，構(gòu)成一個確定C1、…、Cn的線性代數(shù)方程組，寫成

V*C=Y0

其解為C=V-1Y0

其中C=［C1，C2，…，Cn］′;Y0=［y0，Dy0，…，Dny0］′這種形式的矩陣稱為范德蒙特矩陣。在MATLAB的特殊矩陣庫中提供了專門的函數(shù)vander，給出p向量，就可由V=vander(p)生成范德蒙特矩陣。從helpvander知道，需要將它旋轉(zhuǎn)90°，才與本題的形式相符?！鬗ATLAB程序

a=input(′輸入分母系數(shù)向量a=［a1，a2，...］=′);

n=length(a)-1;

Y0=input(′輸入初始條件向量Y0=[y0,Dy0,D2y0，...]=′);

p=roots(a);V=rot90(vander(p));%求根，構(gòu)成范德蒙特矩陣

C=V＼Y0′;%求出C

dt=input(′dt=′);tf=input(′tf=′)

t=0∶dt∶tf;y=zeros(1，length(t));%給出自變量數(shù)據(jù)組

fork=1∶ny=y+c(k)*exp(p(k)*t);end%求各分量的時間函數(shù)并疊加

plot(t，y)，grid下面利用這一程序來解一個三階系統(tǒng)?！舫绦蜻\行結(jié)果運行此程序并輸入：

a=［1，2，9，3］;

dt=0.1;tf=5;而初始值Y0分別?。?，0，0］;［0，1，0］;［0，0，1］用holdon語句使三次運行生成的圖形畫在一幅圖上，得到圖9-3。圖9-3三階系統(tǒng)零輸入分量的解

【例9-1-4】n級放大器，每級的轉(zhuǎn)移函數(shù)均為，求階躍輸入下的過渡過程，畫出n不同時的波形及頻率特性進(jìn)行比較。

解：

◆建模系統(tǒng)的轉(zhuǎn)移函數(shù)為，階躍輸入的拉普拉斯變換為，因此輸出為兩者的乘積，即求Y(s)的拉普拉斯反變換，即可得到輸出過渡過程y(t)。這里我們遇到了一個有多重極點-ωn的H(s)求拉普拉斯反變換的問題。從原理上說，同樣可以用4.3節(jié)中的留數(shù)極點分解法來求。只是在有n重根時，分解出的部分分式的分母將不再是一次極點，而有(s+ωn)、…、(s+ωn)n-1、(s+ωn)n等項，而(s+ωn)-q的反變換可用如下解析式表示:按照這個思路，應(yīng)該先求出Y(s)的極點留數(shù)。注意分母中除了有n個重極點外，還有一個零極點（即1/s），故共有n+1個極點。先用poly函數(shù)求出Y(s)分母多項式的系數(shù)向量，并求出極點及相應(yīng)的留數(shù)：

by=wn^n;ay=［poly(-ones(1，n)*wn)，0］［r，p］=residue(by，ay);再把各個分量疊加在一起。然而，實際上，這樣編程不僅非常麻煩，而且難以得到正確的結(jié)果。其原因是在重極點處，MATLAB的residue算法遇到了病態(tài)問題，數(shù)據(jù)中小小的舍入誤差會使結(jié)果產(chǎn)生很大誤差，即使n取2都得不出正確結(jié)果。為了避開重極點問題，可以有意把極點拉開一些。例如，設(shè)n個極點散布在0.98ωn～1.02ωn之間，那樣就可以全部作為非重極點來列程序。這種處理在工程上是完全沒有問題的，一般電阻的標(biāo)稱誤差為±5%，電容則更大，在工程實踐中，使各個放大器時常數(shù)完全相同是不可能的，即便要把其誤差控制到±2%以內(nèi)也非易事，所以不必自找麻煩，干脆用非重極點的下列程序來求解?！鬗ATLAB程序

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∶N

p0=-linspace(.95，1.05，n)*wn; %將H(s)極點分散布置在±5%區(qū)間

ay=poly(［p0，0］); %由Y(s)的極點（比H(s)多一個零極點）求分母系數(shù)

by=prod(abs(p0));%求Y(s)的分子系數(shù)［r，p］=residue(by，ay);%求Y(s)的留數(shù)極點

fork=1∶n+1%把各部分分式對應(yīng)的時域分量相加

y(n,∶)=y(n,∶)+r(k)*exp(p(t);endfigure(1)，plot(t，y(n,∶));grid，holdon%繪制過渡過程曲線

%下面這幾條語句用來繪制波特圖

figure(2)，bode(prod(abs(p0))，poly(p0));holdonbh=by;ah=poly(p0);%求H(s)的分子分母系數(shù)

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

H=polyval(bh，j*w)./polyval(ah，j*w); %求H(s)在各頻點的值H(jw)

aH=unwrap(angle(H))*180/pi;%求出以度為單位的連續(xù)相角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，holdoff◆程序運行結(jié)果運行此程序，設(shè)N=4，可得1級到4級放大器的過渡過程如圖9-4所示，從中看出輸出信號達(dá)到0.6處所需的時間約為單級時常數(shù)乘以級數(shù)。此程序在N>4時又會出現(xiàn)很大誤差，讀者可自己編寫更好的程序。圖9-4多級放大器的階躍過渡過程為了畫波特圖，程序可按以下步驟進(jìn)行：

(1)求頻率特性。用多項式求值函數(shù)polyval，并且用了元素群運算，輸入頻率數(shù)組作為自變量，一次就求出全部的頻率特性。注意頻率特性是復(fù)數(shù)，通常關(guān)心的是它們的振幅和相角。

(2)振幅和相位特性橫坐標(biāo)都要用對數(shù)坐標(biāo)并且上下對齊。

(3)振幅的縱坐標(biāo)單位應(yīng)化為分貝,這里用了20*log10(abs(H))。

(4)相位的縱坐標(biāo)單位應(yīng)為度，并且應(yīng)連續(xù)變化，不取主角，所以這里加了unwrap命令。在控制系統(tǒng)工具箱中有一個bode命令可以直接完成這些功能，但本書遵循的原則是不用工具箱，以便讓讀者知道工具箱是怎么編程的。當(dāng)然讀者也應(yīng)知道，工具箱中的bode函數(shù)，遠(yuǎn)不止這幾條語句，作為商用軟件庫中的一個正式函數(shù)，它必須考慮自動確定頻率區(qū)間，自動檢查輸入有無錯誤等，程序要復(fù)雜得多。圖9-5繪出了多級放大器的頻率特性。其幅頻特性顯示了低通的特點，隨級數(shù)的增加，通帶減??；相頻特性說明，隨級數(shù)的增加，負(fù)相移成比例地增加。圖

9-5多級放大器的頻率特性

【例9-1-5】設(shè)方波信號的寬度為5s，信號持續(xù)期為10s，試求其在0～20(1/s)頻段間的頻譜特性。如只取0～10(1/s)的頻譜分量（相當(dāng)于通過了一個低通濾波器），求其輸出波形。解：

◆建模設(shè)信號的時域波形為f(t)，在0~10s的區(qū)間外信號為0，則其傅里葉變換為按MATLAB作數(shù)值計算的要求，必須把t分成N份，用相加來代替積分，對于任一給定的ω，可寫成這說明求和的問題可以用f(t)行向量乘以ejωt列向量來實現(xiàn)。此處的Δt是t的增量，在程序中，將用dt來代替。由于要求出一系列不同的ω處的F值都用同一公式，因此可以利用MATLAB中的元素群運算功能，把ω設(shè)成一個行數(shù)組，分別代入本公式左右端的ω中去，寫成(程序中把ω寫成w)F=f*exp(-j*t′*w)·Δt其中，F(xiàn)是與w等長的行向量，exp中的t′是列向量，w是行向量，t′*w是一個矩陣，其行數(shù)與t相同，列數(shù)與w相同。這個矩陣乘式就完成了傅里葉變換。類似地可以得到傅里葉逆變換表示式。由此得到下面的傅里葉變換程序?！鬗ATLAB程序

clear，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)］;%補(bǔ)上負(fù)頻率F=［fliplr(F1)，F(xiàn)1(2∶N)］; %補(bǔ)上負(fù)頻率區(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◆程序運行結(jié)果執(zhí)行這個程序的結(jié)果如圖9-6所示，因為方波含有很豐富的高頻分量，要充分恢復(fù)其原來波形需要很寬的頻帶，實踐中不可能完全做到。圖9-6方波信號的頻譜和取|ω|＜10部分頻譜的逆變換波形(a)方波信號的頻譜;(b)取|ω|<10部分頻譜的逆變換波形9.2離散信號和系統(tǒng)信號可以分為模擬信號和數(shù)字信號。模擬信號用x(t)表示，其中變量t代表時間。離散信號用x(n)表示，其中變量n為整數(shù)并代表時間的離散時刻，因此它也稱為離散時間信號。離散信號是一個數(shù)字的序列，并可以表述為其中，向上的箭頭表示在n=0處的取樣。在MATLAB中，可以用一個向量x來表示一個有限長度的序列。然而這樣一個向量并沒有包含基準(zhǔn)采樣位置的信息。因此，完全地表示x(n)要用x和n兩個向量。例如序列x(n)={2，1，-1，，1，4，3，7}(下面的箭頭為第0個采樣點)，在MATLAB中表示為n=[-3,-2,-1,0,1,2,3,4]，

x=[2,1,-1,5,1,4,3,7]

當(dāng)不需要采樣位置信息時，可以只用x向量來表示。由于內(nèi)存有限，因此MATLAB無法表示無限序列。

【例9-2-1】編寫MATLAB程序來產(chǎn)生下列基本脈沖序列。

(1)單位脈沖序列：起點n0，終點nf，在ns處有一單位脈沖(n0≤ns≤nf)。

(2)單位階躍序列：起點n0，終點nf，在ns前為0，在ns處及以后為1(n0≤ns≤nf)。

(3)實數(shù)指數(shù)序列:x3=(0.9)n。

(4)復(fù)數(shù)指數(shù)序列:x4=e(-0.2+0.3j)n。

解:

◆建模這些基本序列的表達(dá)式比較簡明，編寫程序也不難。對單位脈沖序列，我們提供了直接賦值和邏輯關(guān)系兩種方法,其中用邏輯關(guān)系的編法比較簡潔，讀者從中可看到MATLAB編程的靈活性和技巧性。通常用stem語句來繪制離散序列?！鬗ATLAB程序

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)］;

%也有類似的用邏輯比較語句產(chǎn)生單位階躍序列的方法，留給讀者思考

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));

%注意subplot的輸入變元

subplot(4，2，8)，stem(n4，imag(x4));

line(［0，10］，［0，0］)，%畫橫坐標(biāo)◆程序運行結(jié)果程序運行結(jié)果如圖9-7所示。圖9-7基本脈沖序列的波形

【例9-2-2】離散傅里葉變換的計算。解：◆建模一個時間序列x(n)的離散時間傅里葉變換的定義為

如果序列的長度是有限的，可以把它看做是周期性無限序列中的一個周期，其長度為N。對這個周期性序列可以用離散傅里葉變換（注意少了“時間”兩字）進(jìn)行研究，它的定義為

其中用例9-1-4中的方法，引入矩陣乘法來實現(xiàn)求和運算。用元素群算法來求不同k時的X，把n和k都設(shè)成1×N的行數(shù)組，令nk=n′*k，它就成為N×N的方陣，因而也是N×N方陣。由此得出離散傅里葉變換的算式

MATLAB只能處理有限長度的序列，因此，適合于計算離散傅里葉變換及其逆變換?！鬗ATLAB程序設(shè)有限信號序列xn(n)的長度為Nx，則按定義，求其N點傅里葉變換Xk(k)的程序為:

xn=input(′x=′);Nx=length(xn);N=Nx%取N為x的長度

tic，n=［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; %換算矩陣

Xk=xn*WNnk;toc %DFT系數(shù)向量，即離散傅里葉變換的結(jié)果

plot(abs(Xk))，grid %繪幅頻特性圖在N很大時，這個程序的運算速度比較低。程序中的tic和toc語句用來測試它們之間的程序運行時間。實際上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=′);Nx=length(xn)

%取N為大于Nx而最接近于Nx的2的冪

N=pow2(nextpow2(Nx);

tic，X=fft(xn，N);toc

Nx＜N，x長度不足N的部分，程序會自動補(bǔ)0。要注意X是一個長度為N的復(fù)數(shù)數(shù)組，可以分解出它的振幅和相位，分別繪圖?！舫绦蜻\行結(jié)果按程序提示輸入

x=sin(0.1*［1∶700］)+randn(1，700);所得的fft的幅度特性是一樣的，如圖9-8所示，其中有效信號與噪聲可明顯區(qū)分。在作者的計算機(jī)上，前一程序的時間測試結(jié)果約為8s，后一程序為0.05s。圖9-8信號的fft的振幅頻率特性9.3系統(tǒng)函數(shù)

【例9-3-1】簡單信號流圖模型的矩陣解法[11]?！艚Ｐ盘柫鲌D是用來表示和分析復(fù)雜系統(tǒng)內(nèi)的信號變換關(guān)系的工具。其基本概念如下:

(1)系統(tǒng)中每個信號用圖上的一個節(jié)點表示。如圖中的u、x1、x2。

(2)系統(tǒng)部件對信號實施的變換關(guān)系用有向線段表示，箭尾為輸入信號，箭頭為輸出信號，箭身標(biāo)注對此信號進(jìn)行變換的乘子。如圖上的G1、G2。如果乘子為1，可以不必標(biāo)注。

(3)每個節(jié)點信號的值等于所有指向此節(jié)點的箭頭信號之和，每個節(jié)點信號可以向外輸出給多個部件，其值相同。根據(jù)這幾個概念，可以列出圖9-9的方程如下：x1=u-G2x2，x2=G1x1圖

9-9帶反饋的簡單信號流圖

寫成矩陣方程或x=Qx+Pu移項整理，可以得到求所有未知信號向量x的公式定義系統(tǒng)的傳遞函數(shù)W為輸出信號與輸入信號之比x/u，則W可按下式求得:W=x/u=inv(I-Q)*P因為求這個二階矩陣的逆可以直接用下面的公式：若則所以即u到x1的傳遞函數(shù)為，u到x2的傳遞函數(shù)為。對于階次高的情況，求逆就必須用軟件工具了。如果信號流圖中有G1那樣的符號變量,那么它的求解要用符號運算工具箱?！鬗ATLAB程序

symsG1G2Q=［0，-G2;G1，0］，P=［1;0］

W=inv(eye(2)-Q)*P◆程序運行結(jié)果程序運行的結(jié)果是與前面的結(jié)果相同。這是一個簡單的問題，用一些其他的數(shù)學(xué)方法也能得到同樣的結(jié)果，很出名的“梅森公式”就是用圖形拓樸的方法得到信號流圖的公式，但這個公式非常繁瑣，而且無法機(jī)械化。到現(xiàn)在為止，我們還沒有見過任何一本書籍用矩陣方法來推導(dǎo)這個公式。用矩陣代數(shù)方法的最大好處是可以向任意高的階次、任意復(fù)雜的信號流圖推廣，實現(xiàn)復(fù)雜系統(tǒng)傳遞函數(shù)推導(dǎo)的自動化[11-12]。

【例9-3-2】較復(fù)雜信號流圖模型的矩陣解法?！艚D9-10是一個較復(fù)雜的信號流圖。照上述方法列出它的方程如下: x1=-G4x3+u x2=G1x1-G5x4

x3=G2x2

x4=G3x3將其列為矩陣方程，得到公式W=x/u=inv(I-Q)·P同樣是正確的，不過這里的Q和P分別為4×4和4×1矩陣，用手工求逆是很麻煩的。圖9-10帶雙重反饋的信號流圖◆MATLAB程序

symsG1G2G3G4G5

Q=［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)*P

Pretty(W(4))◆程序運行結(jié)果程序運行的結(jié)果為我們關(guān)心的輸出通常是x4，也就是最后那個傳遞函數(shù)W(4)＝x4)/u，其結(jié)果為當(dāng)系統(tǒng)內(nèi)各個環(huán)節(jié)都是線性集總參數(shù)，因而它們的傳遞函數(shù)Gi都可以表示為s的多項式有理分式時，不管其階次有多高，傳遞函數(shù)W都可以很容易地由計算機(jī)直接自動算出。這個方法還可以推廣到離散系統(tǒng)，用來計算任意復(fù)雜的數(shù)字濾波器系統(tǒng)函數(shù)。9.4頻譜及其幾何意義頻譜分析是信號與系統(tǒng)課程中最重要的內(nèi)容之一，很多讀者在學(xué)習(xí)中感到這部分內(nèi)容很抽象，往往只能從數(shù)學(xué)上了解時域信號與其頻譜間的變換關(guān)系，而沒有理解它的物理意義，而用MATLAB可以幫助讀者建立形象的幾何概念，真正掌握它。首先來看歐拉公式，它是以最簡明的方式建立了信號頻域與時域的關(guān)系:它說明一個最簡單的實余弦信號可以由正、負(fù)兩個Ω0頻率分量合成。在復(fù)平面上，正的Ω0對應(yīng)反時針旋轉(zhuǎn)的向量，負(fù)的Ω0對應(yīng)順時針旋轉(zhuǎn)的向量，當(dāng)這兩個向量的幅度相同而相角符號相反時，就合成為一個在實軸上的向量。它的相角為零，大小按正弦變化，形成了實信號cosΩ0t(如圖9-11所示)。推而廣之，任何實周期信號必然具有正、負(fù)兩組頻率的頻譜成分，正、負(fù)頻率頻譜的幅度對稱而相位反對稱，或者說正、負(fù)頻率的頻譜是共軛的。圖

9-11實序列由對稱的正負(fù)頻率合成

如果頻譜不止這兩項，而是有四項或更多，它們的合成仍然可以用幾何動畫來表示?？梢园衙總€頻譜看做一根長度等于頻譜幅度、按頻率Ω旋轉(zhuǎn)的桿件，頻譜的相加等價于多節(jié)桿件首尾相接，桿件末端的軌跡就描述了生成的時域波形。因為這個端點是在平面上運動，所以它將產(chǎn)生復(fù)信號，在實軸和虛軸上的投影分別為實信號和虛信號。

【例9-4-1】設(shè)計一個演示程序，要求它能把用戶任意給定的四個集總頻譜合成，并生成對應(yīng)的時域信號。

解:◆建模按上述多節(jié)桿合成模型來考慮這個問題。程序設(shè)計主要包括三部分:(1)各頻譜分量的輸入，包括其幅度和頻率(有正負(fù)號);(2)將各分量當(dāng)做轉(zhuǎn)動的桿件首尾相接;(3)記錄多節(jié)桿系末端的軌跡，并畫出圖形?！鬗ATLAB程序

%(1)給頻譜向量賦值

N=input(′N(輸入向量個數(shù)，限定N不大于4