MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例

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

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

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

9.4頻譜及其幾何意義9.5場(chǎng)的計(jì)算和偏微分方程的數(shù)值解9.1連續(xù)信號(hào)和系統(tǒng)

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

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

%(1)單位脈沖信號(hào)

%在t1(t0≤t1≤tf)處有一持續(xù)時(shí)間為dt，面積為1的脈沖信號(hào)，其余時(shí)間均為零

t=［t0∶dt∶tf］;st=length(t);n1=floor((t1-t0)/dt);%求t1對(duì)應(yīng)的樣本序號(hào)

x1=zeros(1，st);%把全部信號(hào)先初始化為零

x1(n1)=1/dt;%給出t1處幅度為1/dt的脈沖信號(hào)

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

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

%(2)單位階躍信號(hào)

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

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

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

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

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

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))%繪圖◆程序運(yùn)行結(jié)果

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

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

x2、

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

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

解:

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

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

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

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

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

b=input(′多項(xiàng)式分子系數(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í)間數(shù)組

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

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

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

end

plot(t，h)，grid

else，

end◆程序運(yùn)行結(jié)果例如，給出系統(tǒng)傳遞函數(shù)為求沖擊響應(yīng)。根據(jù)程序提問(wèn)依次輸入: 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，而題中給出的是極點(diǎn)向量p=［0，-1，-2，-5］，因此這里用poly函數(shù)來(lái)作轉(zhuǎn)換。

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

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

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

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

V*C=Y0

其解為C=V-1Y0

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

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í)間函數(shù)并疊加

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

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

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

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

解：

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

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

clear，clf，N=input(′輸入放大器級(jí)數(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)極點(diǎn)分散布置在±5%區(qū)間

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

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

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

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

%下面這幾條語(yǔ)句用來(lái)繪制波特圖

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)在各頻點(diǎn)的值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◆程序運(yùn)行結(jié)果運(yùn)行此程序，設(shè)N=4，可得1級(jí)到4級(jí)放大器的過(guò)渡過(guò)程如圖9-4所示，從中看出輸出信號(hào)達(dá)到0.6處所需的時(shí)間約為單級(jí)時(shí)常數(shù)乘以級(jí)數(shù)。此程序在N>4時(shí)又會(huì)出現(xiàn)很大誤差，讀者可自己編寫(xiě)更好的程序。圖9-4多級(jí)放大器的階躍過(guò)渡過(guò)程為了畫(huà)波特圖，程序可按以下步驟進(jìn)行：

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

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

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

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

9-5多級(jí)放大器的頻率特性

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

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

clear，tf=10;N=256;

t=linspace(0，tf，N); %給出時(shí)間分割

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)］;%給出信號(hào)(此處是方波)

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;%對(duì)中段頻譜求傅里葉逆變換

subplot(1，2，2)，plot(t，f，t，f1，′linewidth′，1.5)，grid◆程序運(yùn)行結(jié)果執(zhí)行這個(gè)程序的結(jié)果如圖9-6所示，因?yàn)榉讲ê泻茇S富的高頻分量，要充分恢復(fù)其原來(lái)波形需要很寬的頻帶，實(shí)踐中不可能完全做到。圖9-6方波信號(hào)的頻譜和取|ω|＜10部分頻譜的逆變換波形(a)方波信號(hào)的頻譜;(b)取|ω|<10部分頻譜的逆變換波形9.2離散信號(hào)和系統(tǒng)信號(hào)可以分為模擬信號(hào)和數(shù)字信號(hào)。模擬信號(hào)用x(t)表示，其中變量t代表時(shí)間。離散信號(hào)用x(n)表示，其中變量n為整數(shù)并代表時(shí)間的離散時(shí)刻，因此它也稱(chēng)為離散時(shí)間信號(hào)。離散信號(hào)是一個(gè)數(shù)字的序列，并可以表述為其中，向上的箭頭表示在n=0處的取樣。在MATLAB中，可以用一個(gè)向量x來(lái)表示一個(gè)有限長(zhǎng)度的序列。然而這樣一個(gè)向量并沒(méi)有包含基準(zhǔn)采樣位置的信息。因此，完全地表示x(n)要用x和n兩個(gè)向量。例如序列x(n)={2，1，-1，，1，4，3，7}(下面的箭頭為第0個(gè)采樣點(diǎn))，在MATLAB中表示為n=[-3,-2,-1,0,1,2,3,4]，

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

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

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

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

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

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

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

解:

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

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

n3=n0∶nf;x3=(0.9).^n3;%實(shí)數(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的輸入變?cè)?/p>

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

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

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

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

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

MATLAB只能處理有限長(zhǎng)度的序列，因此，適合于計(jì)算離散傅里葉變換及其逆變換。◆MATLAB程序設(shè)有限信號(hào)序列xn(n)的長(zhǎng)度為Nx，則按定義，求其N(xiāo)點(diǎn)傅里葉變換Xk(k)的程序?yàn)?

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

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)生一個(gè)含nk值的N×N維矩陣

WNnk=WN.^nk; %換算矩陣

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

plot(abs(Xk))，grid %繪幅頻特性圖在N很大時(shí)，這個(gè)程序的運(yùn)算速度比較低。程序中的tic和toc語(yǔ)句用來(lái)測(cè)試它們之間的程序運(yùn)行時(shí)間。實(shí)際上MATLAB已提供了快速離散傅里葉變換的函數(shù)fft，可直接調(diào)用。其調(diào)用格式為X=fft(x，N)其中，x是輸入的時(shí)間序列，N是傅里葉變換取的點(diǎn)數(shù)。若省略N，則它自動(dòng)把x的長(zhǎng)度作為N。當(dāng)N取2的冪時(shí)，變換速度最快，所以要提高fft函數(shù)的運(yùn)行速度，程序應(yīng)編寫(xiě)如下：

xn=input(′x=′);Nx=length(xn)

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

N=pow2(nextpow2(Nx);

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

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

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

【例9-3-1】簡(jiǎn)單信號(hào)流圖模型的矩陣解法[11]。◆建模信號(hào)流圖是用來(lái)表示和分析復(fù)雜系統(tǒng)內(nèi)的信號(hào)變換關(guān)系的工具。其基本概念如下:

(1)系統(tǒng)中每個(gè)信號(hào)用圖上的一個(gè)節(jié)點(diǎn)表示。如圖中的u、x1、x2。

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

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

9-9帶反饋的簡(jiǎn)單信號(hào)流圖

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

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

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

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

x3=G2x2

x4=G3x3將其列為矩陣方程，得到公式W=x/u=inv(I-Q)·P同樣是正確的，不過(guò)這里的Q和P分別為4×4和4×1矩陣，用手工求逆是很麻煩的。圖9-10帶雙重反饋的信號(hào)流圖◆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))◆程序運(yùn)行結(jié)果程序運(yùn)行的結(jié)果為我們關(guān)心的輸出通常是x4，也就是最后那個(gè)傳遞函數(shù)W(4)＝x4)/u，其結(jié)果為當(dāng)系統(tǒng)內(nèi)各個(gè)環(huán)節(jié)都是線(xiàn)性集總參數(shù)，因而它們的傳遞函數(shù)Gi都可以表示為s的多項(xiàng)式有理分式時(shí)，不管其階次有多高，傳遞函數(shù)W都可以很容易地由計(jì)算機(jī)直接自動(dòng)算出。這個(gè)方法還可以推廣到離散系統(tǒng)，用來(lái)計(jì)算任意復(fù)雜的數(shù)字濾波器系統(tǒng)函數(shù)。9.4頻譜及其幾何意義頻譜分析是信號(hào)與系統(tǒng)課程中最重要的內(nèi)容之一,很多讀者在學(xué)習(xí)中感到這部分內(nèi)容很抽象,往往只能從數(shù)學(xué)上了解時(shí)域信號(hào)與其頻譜間的變換關(guān)系,而沒(méi)有理解它的物理意義,而用MATLAB可以幫助讀者建立形象的幾何概念,真正掌握它。首先來(lái)看歐拉公式,它以最簡(jiǎn)明的方式建立了信號(hào)頻域與時(shí)域的關(guān)系:它說(shuō)明一個(gè)最簡(jiǎn)單的實(shí)余弦信號(hào)可以由正、負(fù)兩個(gè)Ω0頻率分量合成。在復(fù)平面上,正的Ω0對(duì)應(yīng)反時(shí)針旋轉(zhuǎn)的向量,負(fù)的Ω0對(duì)應(yīng)順時(shí)針旋轉(zhuǎn)的向量,當(dāng)這兩個(gè)向量的幅度相同而相角符號(hào)相反時(shí),就合成為一個(gè)在實(shí)軸上的向量。它的相角為零,大小按正弦變化,形成了實(shí)信號(hào)cosΩ0t(如圖9-11所示)。推而廣之,任何實(shí)周期信號(hào)必然具有正、負(fù)兩組頻率的頻譜成分,正、負(fù)頻率頻譜的幅度對(duì)稱(chēng)而相位反對(duì)稱(chēng),或者說(shuō)正、負(fù)頻率的頻譜是共軛的。圖9-11實(shí)序列由對(duì)稱(chēng)的正負(fù)頻率合成【例9-4-1】設(shè)計(jì)一個(gè)演示程序,要求它能把用戶(hù)任意給定的四個(gè)集總頻譜合成,并生成對(duì)應(yīng)的時(shí)域信號(hào)。解:◆建模按上述多節(jié)桿合成模型來(lái)考慮這個(gè)問(wèn)題。程序設(shè)計(jì)主要包括三部分:(1)各頻譜分量的輸入,包括其幅度和頻率(有正負(fù)號(hào));(2)將各分量當(dāng)做轉(zhuǎn)動(dòng)的桿件首尾相接;(3)記錄多節(jié)桿系末端的軌跡,并畫(huà)出圖形?！舫绦蜻\(yùn)行結(jié)果運(yùn)行程序并按提示輸入數(shù)據(jù),如果只取兩個(gè)幅度相等而頻率符號(hào)相反的集總頻譜,那么將得到與圖9-11相仿的結(jié)果?，F(xiàn)取四個(gè)集總頻譜,輸入

a(1)=1,w(1)=-1;

a(2)=1,w(2)=-1;

a(3)=0.5,w(3)=3;

a(4)=0.5,w(1)=-4;圖9-12四個(gè)頻譜向量組成的多節(jié)桿及其端點(diǎn)軌跡圖9-13四個(gè)集總頻譜生成的復(fù)信號(hào)及其實(shí)信號(hào)分量(a)合成的復(fù)平面信號(hào);(b)合成的虛部信號(hào);(c)合成的實(shí)部信號(hào)輸入頻譜的幅度可以是負(fù)數(shù),也可以是虛數(shù),甚至可以是復(fù)數(shù),它不僅反映了頻譜的大小,還反映了該向量的起始相位;頻譜的頻率則只能是可正可負(fù)的實(shí)數(shù),正頻率和負(fù)頻率以及在該頻率上頻譜的意義在此不言自明。讀者可以做各種各樣的試驗(yàn),例如當(dāng)兩組頻率具有倍頻關(guān)系時(shí),得到的是周期信號(hào);如果頻率比是任意小數(shù),那將得出非周期的信號(hào)。另外,這樣的演示適用于集總頻譜,對(duì)于分布的頻譜密度,就要把它想象為若干小的集總頻譜的疊合。總之,有了這樣的形象演示,便可以大大擴(kuò)展時(shí)域信號(hào)與頻域譜之間關(guān)系的思維空間。可見(jiàn),只有在復(fù)信號(hào)平面上,才能看到頻率的正負(fù)。許多人往往看不到負(fù)頻率頻譜的意義,這是因?yàn)樗麄兛偸峭Ａ粼趯?shí)信號(hào)的范疇來(lái)思考問(wèn)題。要想知道象有沒(méi)有鼻子(負(fù)頻率),必須去摸象頭(復(fù)信號(hào)),只摸象腿(實(shí)信號(hào))是得不到真切的認(rèn)識(shí)(信號(hào)理論普遍而科學(xué)的規(guī)則)的。9.5場(chǎng)的計(jì)算和偏微分方程的數(shù)值解現(xiàn)代科學(xué)和工程中大量遇到偏微分方程,這類(lèi)方程通常用來(lái)表示分布在三維空間的