MATLAB仿真在通信與電子工程中的應(yīng)用 - 第2章 MATLAB仿真

第2章MATLAB仿真2.1MATLAB介紹

2.2MATLAB編程仿真2.3Simulink仿真

2.1MATLAB介紹 2.1.1MATLAB的特點

MATLAB是MathWorks公司開發(fā)的一種跨平臺的，用于矩陣數(shù)值計算的簡單高效的數(shù)學語言，與其它計算機高級語言如C，C++，F(xiàn)ortran，Basic，Pascal等相比，MATLAB語言編程要簡潔得多，編程語句更加接近數(shù)學描述，可讀性好，其強大的圖形功能和可視化數(shù)據(jù)處理能力也是其它高級語言望塵莫及的。對于具有任何一門高級語言基礎(chǔ)的讀者來說，學習MATLAB十分容易。 但是，要用好MATLAB卻不是在短時間就可以達到的。這并不是因為MATLAB語言復(fù)雜難懂，而是實際問題的求解往往更多的是需要使用者具備數(shù)學知識和專業(yè)知識。MATLAB使得人們擺脫了常規(guī)計算機編程的繁瑣，讓人們能夠?qū)⒋蟛糠志ν度氲窖芯繂栴}的數(shù)學建模上?？梢哉f，應(yīng)用MATLAB這一數(shù)學計算和系統(tǒng)仿真的強大工具，可以使科學研究的效率得以成百倍的提高。 目前，MATLAB已經(jīng)廣泛用于理工科大學從高等數(shù)學到幾乎各門專業(yè)課程之中，成為這些課程進行虛擬實驗的有效工具。在科研部門，MATLAB更是極為廣泛地得到應(yīng)用，成為全球科學家和工程師進行學術(shù)交流首選的共同語言。在國內(nèi)外許多著名學術(shù)期刊上登載的論文，大部分的數(shù)值結(jié)果和圖形都是借助MATLAB來完成的。

與其它高級語言相比較，MATLAB具有獨特的優(yōu)勢： (1)MATLAB是一種跨平臺的數(shù)學語言。采用MATLAB編寫的程序可以在目前所有的操作系統(tǒng)上運行(只要這些系統(tǒng)上安裝了MATLAB平臺)。MATLAB程序不依賴于計算機類型和操作系統(tǒng)類型。 (2)MATLAB是一種超高級語言。MATLAB平臺本身是用C語言寫成的，其中匯集了當前最新的數(shù)學算法庫，是許多專業(yè)數(shù)學家和工程學者多年的勞動結(jié)晶。使用MATLAB意味著站在巨人的肩膀上觀察和處理問題，所以在編程效率，程序的可讀性、可靠性和可移植性上遠遠超過了常規(guī)的高級語言。這使得MATLAB成為了進行科學研究和數(shù)值計算的首選語言。

(3)MATLAB語法簡單，編程風格接近數(shù)學語言描述，是數(shù)學算法開發(fā)和驗證的最佳工具。MATLAB以復(fù)數(shù)矩陣運算為基礎(chǔ)，其基本編程單位是矩陣，使得編程簡單，而功能極為強大。對于常規(guī)語言中必須使用許多語句才能實現(xiàn)的功能，如矩陣分解、矩陣求逆、積分、快速傅立葉變換，甚至串口操作、聲音的輸入輸出等，在MATLAB中均用一兩句指令即可實現(xiàn)。而且，MATLAB中的數(shù)值算法是經(jīng)過千錘百煉的，比用戶自己編程實現(xiàn)的算法的可信度和可靠性都大為提高。 (4)MATLAB計算精度很高。MATLAB中數(shù)據(jù)是以雙精度存儲的，一個實數(shù)采用8字節(jié)存儲，而一個復(fù)數(shù)則采用16字節(jié)存儲。通常矩陣運算精度高達1015以上，完全能夠滿足一般工程和科學計算的需要。與其它語言相比，MATLAB對計算機內(nèi)存、硬盤空間的要求也是比較高的。 (5)MATLAB具有強大的繪圖功能。利用MATLAB的繪圖功能，可以輕易地獲得高質(zhì)量的(印刷級)曲線圖。具有多種形式來表達二維、三維圖形，并具有強大的動畫功能，可以非常直觀地表現(xiàn)抽象的數(shù)值結(jié)果。這也是MATLAB廣為流行的重要原因之一。 (6)MATLAB具有串口操作、聲音輸入輸出等硬件操控能力。隨著版本的提高，這種能力還會不斷加強，使得人們利用計算機和實際硬件相連接的半實物仿真的夢想得以輕易實現(xiàn)。 (7)MATLAB程序可以直接映射為DSP芯片可接受的代碼，大大提高了現(xiàn)代電子通信設(shè)備的研發(fā)效率。 (8)MATLAB的程序執(zhí)行效率比其它語言低。MATLAB程序通常是解釋執(zhí)行的，在執(zhí)行效率和速度上低于其它高級語言，當然如果對執(zhí)行效率有特別要求，可以采用C語言編制算法，然后通過MATLAB接口在MATLAB中執(zhí)行。事實上，MATLAB自帶的許多內(nèi)部函數(shù)均是用C語言編寫并編譯的，因此利用MATLAB內(nèi)部函數(shù)的程序部分運行速度并不比其它語言中相應(yīng)函數(shù)低。 本書給出的程序和例子均在MATLABReleaseR13（6.5.1）（完全安裝）版本下驗證通過。

2.1.2MATLAB快速入門 1.MATLAB的啟動和退出 本書以Windows2000下的MATLAB6.5.1為例，在其它Windows版本下操作是類似的。假定已經(jīng)正確安裝了MATLAB的Windows版本，從Windows的“開始|程序”菜單中找到MATLAB圖標，單擊該圖標即可進入MATLAB環(huán)境。進入MATLAB后操作環(huán)境通常顯示的是MATLAB的命令窗口，在該窗口中可以輸入各種MATLAB命令和語句，通過命令方式與MATLAB進行交互。雖然MATLAB也像通常的Windows程序一樣提供了菜單和快捷工具欄，通過它們可以很方便地對MATLAB進行操作，但是建議讀者盡可能使用命令方式去操作MATLAB，雖然剛開始可能覺得不太方便，但是與菜單和快捷工具方式相比，命令方式的功能最為強大，也最能體現(xiàn)MATLAB的精髓，而且命令方式本身也是跨平臺的。

MATLAB命令窗口中，輸入命令的提示符為“>>”。 在MATLAB命令窗口中，鍵入命令“quit”或“exit”并按回車鍵(回車)，即可退出MATLAB環(huán)境。 2.在MATLAB中如何獲取幫助

MATLAB幾乎涉及了所有工程領(lǐng)域的數(shù)學問題，沒有一本書能夠完全覆蓋讀者所需要解決的問題。因為MATLAB幫助文檔是惟一完全覆蓋MATLAB功能函數(shù)的權(quán)威技術(shù)文件，所以，善于利用MATLAB的幫助文檔將是非常重要的。 在MATLAB命令窗口中，使用命令“intro”可以進入MATLAB簡介演示。而使用命令“demo”可以打開MATLAB的演示窗口，其中包含了大量MATLAB程序的演示實例，對于初學者具有很高的參考價值。如果讀者需要打開MATLAB的幫助文檔，使用命令“helpwin”即可。如果要用瀏覽器打開html形式的幫助文檔，可使用命令“doc”。若需要打開MATLAB某條命令或函數(shù)用法的html幫助文檔，則只需輸入命令

>>doc命令或函數(shù)名［回車］即可。 如果要在命令窗口中顯示幫助信息，則只需輸入命令“help”即可。若需要了解MATLAB某條命令或函數(shù)的用法，使用命令 >>help命令或函數(shù)名［回車］即可。例如，要查詢自然對數(shù)函數(shù)命令“l(fā)og”的用法，應(yīng)使用命令 >>helplog［回車］ 另外，MATLAB幫助文檔還以PDF電子文件格式提供，存放在MATLAB的安裝目錄中，讀者可以用AcrobatReader等軟件閱讀。讀者也可以通過介紹MATLAB的若干網(wǎng)站獲取和交流應(yīng)用MATLAB的心得體會。

3.在MATLAB命令窗口中計算數(shù)學表達式

MATLAB語言是一種解釋性語言，它提供了方便的演算紙式的數(shù)學計算方式。在MATLAB命令窗口中輸入數(shù)學表達式，然后回車即可得出計算結(jié)果。MATLAB的數(shù)學表達式與數(shù)學公式表達極為相似，也非常類似于C語言的表達。例如，計算表達式2sin(0.3π)/(1+[KF(]5[KF)])的值，在MATLAB命令窗口中輸入以下語句并回車就能得到結(jié)果。

>>2*sin(0.3*pi)/(1+sqrt(5))［回車］

ans= 0.5000 >> 以上語句中，“pi”是MATLAB已定義的常數(shù)，即圓周率π；“ans”是表達式計算結(jié)果的默認存儲變量。MATLAB定義的常數(shù)中，常用的有以下一些：

ans最新表達式計算結(jié)果的默認存儲變量

inf表示正無窮大+∞

NaN非數(shù)，例如0/0將得到非數(shù)

i或j

虛數(shù)單位,即

eps從1.0到下一個最大浮點數(shù)的距離,常用來作為浮點計算相對誤差使用

pi

圓周率π realminMATLAB所能表示的最小正實數(shù)

realmaxMATLAB所能表示的最大正實數(shù)

MATLAB中常用的算術(shù)運算符有+、-、*、/、\、^、′等，含義如下： +數(shù)量加法,矩陣加法 -數(shù)量減法,矩陣減法 *數(shù)量乘法,矩陣乘法 .*數(shù)組乘法 /數(shù)量除法,矩陣右除 ＼數(shù)量除法(左除),矩陣左除 ./數(shù)組除法 ^矩陣乘方 .^數(shù)組乘方 ′矩陣的共軛轉(zhuǎn)置,對于一個復(fù)數(shù)而言將得到其共軛復(fù)數(shù) .′矩陣轉(zhuǎn)置(不共軛) 關(guān)于這些算術(shù)運算符的詳細用法以及示例，可以用

>>docarithmeticoperators［回車］ 打開html格式的幫助文檔查看。 4.在MATLAB命令窗口中輸入簡單矩陣 例如，輸入矩陣

時可以采用下面的方法：

>>A=［1,2,3;4,5,6;7,8,9］［回車］

A= 123 456 789(2-1)

MATLAB中可以采用逗號或空格來分隔矩陣中的列元素，而采用分號或回車符來分隔矩陣的行，整個矩陣包含在方括號“［］”內(nèi)。采用命令“whos”可以查看用戶在MATLAB工作空間(內(nèi)存)中所存儲的變量情況。使用命令“workspace”就可以打開工作空間瀏覽器窗口，雙擊其中的變量可以對其值進行修改。 命令“clear變量名”可以清除相應(yīng)的變量，而命令“clear”則清除所有用戶自定義的變量。系統(tǒng)的默認變量是不會被“clear”清除的。為了避免前面的程序?qū)罄m(xù)程序的影響，通常在程序的開始使用“clear”語句來復(fù)位MATLAB的內(nèi)存空間。使用命令“clc”可以清除命令窗口中的顯示字符。使用命令“home”可以使得命令窗口中提示符光標回到窗口的左上角。這兩個命令僅僅影響屏幕的顯示，不會清除內(nèi)存中的變量。

5.MATLAB的語句和變量

MATLAB語句可以有兩種形式： （1)表達式； （2)變量名＝表達式。 在第一種形式中，表達式計算的值將存放于默認變量“ans”中；而在第二種形式中，表達式的值將存放于變量名所指定的變量中。MATLAB中多條語句可以在一行內(nèi)書寫，以逗號“,”或分號“;”相互隔開。如果是以分號隔開的，則計算結(jié)果不顯示在屏幕上，否則回車后將顯示計算結(jié)果。例如，計算1+2+3+…+100的值以及100的階乘100!。 >>s=［1:100］;sum(s),p=prod(s)［回車］

ans= 5050 p= 9.3326e+157 MATLAB中的變量是區(qū)分大小寫的，變量、函數(shù)名必須以字母開頭，其后最多可接19個字母、數(shù)字或下劃線。例如：a和A是不同的變量，p121-6，yinyue3是合法的變量名，而3sd則是非法的變量。在變量或函數(shù)的命名時，應(yīng)該養(yǎng)成良好的命名習慣，命名不要和MATLAB中的系統(tǒng)函數(shù)或變量相同。特別地，如果計算中存在復(fù)數(shù)運算，那么就應(yīng)該避免采用“i”和“j”作為循環(huán)變量。

MATLAB中可以方便地進行復(fù)數(shù)運算，例如計算

，其中，a=15+j3，b=5ej2。

可在命令窗口中輸入：

>>a=15+j*3,b=5*exp(j*2),(a.^2+b).^(1/5)［回車］ 6.繪制簡單的函數(shù)曲線

MATLAB提供了極為便利的數(shù)據(jù)可視化手段，可以作出任意函數(shù)的圖像。作為快速入門，在此以一個二維作圖為例，作出函數(shù)y=e-x/10sinx在x∈［-1,10］范圍的圖像。 >>x=-1:0.1:20; %定義x的范圍和步進［回車］ >>y=exp(-x./10).*sin(x);%計算函數(shù)值［回車］ >>plot(x,y);grid;%作出函數(shù)圖像,并在坐標上畫出網(wǎng)格［回車］結(jié)果如圖2-1所示。注意，在程序語句中,以百分號“%”開始的是注釋部分。關(guān)于“plot”的詳細用法和例子可以用“docplot”命令查看幫助文檔。圖2-1函數(shù)y=e-x/10sinx在x∈［-1,10］范圍的圖像 7.編寫簡單的MATLAB程序

MATLAB是解釋性語言，輸入一行語句后回車，就會立即執(zhí)行得出結(jié)果。如果要實現(xiàn)比較復(fù)雜的功能，單靠一條一條地在命令窗口中輸入指令執(zhí)行，效率是很低的。如何解決這個問題呢？為此MATLAB提供了擴展名為“.m”的文本文件，在文件中事先寫入一行行的MATLAB命令，存盤后從MATLAB的命令窗口調(diào)入執(zhí)行（類似于DOS下的批處理），這種文件稱為底稿文件或MATLAB腳本文件。 用前面作圖的三條語句為例，在文本編輯器中輸入這些語句，然后將文件存盤，例如命名文件名為“my1stprg.m”，保存于MATLAB的默認工作路徑中。然后回到MATLAB命令窗口，在提示符“>>”下鍵入文件名(可以省略擴展名)后回車，即可運行程序得到結(jié)果，即>>my1stprg［回車］在任何一種純文本編輯器中均可以書寫MATLAB程序，只要以“.m”為擴展名保存，即可在MATLAB中調(diào)用運行。MATLAB命令窗口中，在提示符“>>”后輸入的命令語句要回車才能執(zhí)行，本書以后的例子中，將省略［回車］字樣，在此提醒讀者注意。 2.1.3MATLAB程序設(shè)計 1.M文件簡介

MATLAB除了如前所述的在命令窗口進行的直接交互的指令操作方式外，另外一種更為重要的工作方式就是M文件的編程工作方式。M文件有兩種形式，一種是腳本文件，另一種是函數(shù)文件。M文件的擴展名為“.m”。M文件可以通過任何純文本編輯器進行編輯，MATLAB中也有自帶的文本編輯器，使用“edit”命令即可開啟。 2.程序控制流語句 任何計算機語言，只要利用順序結(jié)構(gòu)、循環(huán)結(jié)構(gòu)以及分支結(jié)構(gòu)，就可以完成任何程序功能。在MATLAB中也有這三種基本的程序結(jié)構(gòu)。但值得注意的是，由于MATLAB語言矩陣計算功能十分強大，常常僅使用順序結(jié)構(gòu)借以矩陣的邏輯運算就可以完成計算任務(wù)。由于循環(huán)結(jié)構(gòu)和分支結(jié)構(gòu)在MATLAB語言中的運行速度相對較慢，因此在算法優(yōu)化的編程中應(yīng)當盡可能避免使用，而代之以矩陣運算，從而提高程序運行速度，簡化程序代碼，并使得程序代碼更加接近于數(shù)學表達。當然，采用矩陣編程的編程方法需要讀者具有更多的關(guān)于線性代數(shù)和矩陣數(shù)學的知識和思維方式。 順序結(jié)構(gòu)是MATLAB中最常用的程序結(jié)構(gòu)，也是執(zhí)行效率最高的程序結(jié)構(gòu)。順序結(jié)構(gòu)的語句是按照書寫的前后順序來執(zhí)行的。

MATLAB用于循環(huán)結(jié)構(gòu)的語句有兩種：“for…end”循環(huán)和“while…end”循環(huán)。 1)for…end語句 利用“helpfor”或“docfor”可以獲得關(guān)于該語句的使用手冊?！癴or…end”語句適用于循環(huán)次數(shù)確定的情況，將循環(huán)變量的初值、判別和變化放在循環(huán)開頭?！癴or…end”語句的調(diào)用形式是：

forv=表達式 語句1; … 語句n; end

例如，最簡單的for…end循環(huán)：

fork=1:10 x(k)=k.^2; end 事實上，采用矩陣思想也可以獲得相同的結(jié)果，但編程更加簡單明了：

k=1:10;x=k.^2

x= 149162536496481100

絕大部分循環(huán)都可以遵從這樣的方法變成向量化(矩陣化)的算法，避免采用循環(huán)語句，從而大大提高程序的執(zhí)行效率。 2)while…end語句 對于循環(huán)次數(shù)不能預(yù)先確定，而是由某個邏輯條件來控制循環(huán)次數(shù)的情況，MATLAB提供了“while…end”來實現(xiàn)。和“for…end”語句類似，“while…end”語句也允許嵌套。 “while…end”語句的一般形式是：

whilev=表達式 語句1; … 語句n; end 例如，求當整數(shù)n的階乘值是一個50位數(shù)的第一個數(shù)時，n為多少？（程序jiechen50.m）

n=1;%n的初值

whileprod(1:n)<1e50 n=n+1; end n n= 42 3)條件分支結(jié)構(gòu)“if”、“break”、“switch”語句 “if”分支結(jié)構(gòu)的一般形式是：

if表達式 語句段1;

else

語句段2;

end

“break”語句一般出現(xiàn)在循環(huán)體中，它表示跳出循環(huán)。仍然以前例加以說明：求整數(shù)n的階乘值是一個50位數(shù)的第一個數(shù)時，n為多少？現(xiàn)在使用“if”和“break”語句來編程。 (程序jiechen50b.m） n=1;%n的初值

while1%構(gòu)成一個死循環(huán)

n=n+1;

ifprod(1:n)>1e50,break;end%滿足條件,則跳出循環(huán)

end n%顯示結(jié)果

n= 42 “switch”語句來實現(xiàn)多重分支結(jié)構(gòu)。其用法是：

switch開關(guān)表達式

case表達式1 語句,…,語句

case{表達式1,表達式2,表達式3,…} 語句,…,語句 …

otherwise

語句,…,語句

end

3.數(shù)據(jù)和文件的輸入輸出 1)“input”指令 “input”指令提示用戶從鍵盤輸入數(shù)據(jù)、字符串或表達式，并接受該輸入?！癷nput”指令的調(diào)用格式有兩種：

user-entry=input(′prompt′)%輸入數(shù)據(jù)或表達式

user-entry=input(′prompt′,′s′)%輸入字符串 例如： >>a=input(′請輸入矩陣或表達式,賦值到a:′)[DW]%輸入數(shù)據(jù)的例子 請輸入矩陣或表達式,賦值到a:［12;34］ a= 12 34

>>a=input(′請輸入矩陣或表達式,賦值到a:′)[DW]%輸入表達式的例子 請輸入矩陣或表達式,賦值到a:sin(1)+8 a= 8.8415 >>s=input(′請輸入一個字符串:′,′s′)[DW]%輸入字符串的例子 請輸入一個字符串:thisisastring

s=

thisisastring

2)“pause”指令 “pause”指令可以使程序暫停運行，等待用戶按任意鍵繼續(xù)?！皃ause”指令主要用于程序調(diào)試或顯示中間結(jié)果。“pause(n)”指令使得程序暫停執(zhí)行n秒。 3)利用文件輸入輸出數(shù)據(jù) 指令“save”和“l(fā)oad”用于MATLAB與磁盤的數(shù)據(jù)交換?！皊ave”指令將MATLAB工作空間的數(shù)據(jù)存入磁盤；“l(fā)oad”指令則將數(shù)據(jù)從磁盤送入到MATLAB工作空間。舉例說明如下：

clear; A=［1,2;3,4］; b=81; save(′c:＼mymatlabdat.mat′,′A′,′b′);%將變量A和b保存在文件C:＼mymatlabdat.m

at中

clear%清除內(nèi)存變量

whos%查看,應(yīng)該無變量顯示

load(′c:＼mymatlabdat.mat′);%調(diào)入數(shù)據(jù)文件c:＼mymatlabdat.mat

whos%顯示調(diào)入的變量

NameSizeBytesClass

A

2x232doublearray b1x18doublearray Grandtotalis5elementsusing40bytes 4.MATLAB編程特點

MATLAB有兩種工作方式：一種是在命令窗口進行的指令操作方式，前面我們所使用的大部分就是這種方式；另外一種是M文件編程方式，這種方式特別適合于復(fù)雜問題的求解，是MATLAB高級應(yīng)用的一種常用方式。MATLAB編程中，要特別注意程序的書寫風格，一個好的程序，必須思路清晰，注釋詳細，而且是運行速度較快的。M文件編程中開頭的注釋行將作為該程序的幫助信息，可以在命令窗口中用“help”命令顯示出來。例如，編程計算函數(shù)

f(x)=x3+x+lnxsinx+

當x=1,3,5時的值。 編寫腳本文件calcfx.m如下： 程序2-1 %這是開頭的注釋行,可以用help命令顯示 %程序calcfx.m的功能是計算表達式……

clear;%清除MATLAB內(nèi)存空間,這一命令常用于腳本文件的首句 %以避免前面命令在內(nèi)存空間形成的變量的意外影響

k=1;%數(shù)組下標變量,在MATLAB編程中,注釋要盡可能詳盡!

Int-F=inline(′t′,′t′);%用inline函數(shù)建立積分的被積函數(shù)

forx=［1,3,5］

f-x(k)=x^3+x+log(x)*sin(x)+quad8(Int-F,0,x);%計算表達式值

k=k+1;%數(shù)組下標加1

end f-x%顯示計算結(jié)果

文件編輯存盤后,在MATLAB命令窗口執(zhí)行:

>>helpcalcfx

這是開頭的注釋行,可以用help命令顯示 程序calcfx.m的功能是計算表達式…… >>calcfx

f-x= 2.500034.6550140.9567 >> 5.MATLAB函數(shù)編程 如果M文件的第一行是以關(guān)鍵字“function”開頭的，則就是函數(shù)文件。函數(shù)文件是MATLAB程序設(shè)計的主流，MATLAB自身所帶的許多函數(shù)(指令)都是由相應(yīng)的函數(shù)文件來定義的。函數(shù)文件好像一個黑箱，將數(shù)據(jù)送進去，經(jīng)過函數(shù)處理，然后將結(jié)果數(shù)據(jù)輸出。函數(shù)文件和腳本文件在內(nèi)存使用上存在重要的區(qū)別，函數(shù)文件內(nèi)部所定義的變量僅僅在該函數(shù)文件內(nèi)部有效，函數(shù)返回后這些內(nèi)部變量將自動被清除，也就是說，函數(shù)內(nèi)部所定義的變量僅在函數(shù)內(nèi)部起作用，是局部變量。 而腳本文件中所定義的變量，在使用“clear”命令清除之前，始終存在于工作空間當中，是全局變量。另外，函數(shù)文件的文件名必須和函數(shù)名相同，而且注意切記不要和MATLAB已經(jīng)定義的系統(tǒng)函數(shù)和其它自定義的函數(shù)同名。下面是函數(shù)編程的例子，函數(shù)的文件名與函數(shù)名相同，為stat.m。

程序2-2

function［mean,stdev］=stat(x) n=length(x);%求輸入向量x的長度

mean=sum(x)/n;%求平均值

stdev=sqrt(sum((x-mean).^2/n));%求均方根值 6.測定程序執(zhí)行時間和時間分配 利用“tic”和“toc”指令可以對程序段的執(zhí)行時間進行測定，從而估計出程序執(zhí)行效率，并找出改進程序、提高效率的方法?！皌ic”用于計時開始，而“toc”用于計時結(jié)束并顯示計時結(jié)果。MATLAB還提供了對程序執(zhí)行的耗時剖析功能“profile”指令。用戶通過調(diào)用該功能函數(shù)，可以輕松地觀察程序中各條語句的執(zhí)行耗時情況，從而為提高程序運行效率的改進思路提供參考依據(jù)。(詳細情況請參見在線幫助。) 7.提高程序執(zhí)行速度的原則

MATLAB是一種解釋性語言，它與C語言等編譯性語言有著相當大的區(qū)別，如果按照C語言的思路去編寫MATLAB程序，那么執(zhí)行效率肯定不是很理想。要提高編程的執(zhí)行效率，則一定要根據(jù)MATLAB的特點來編寫程序，具體有以下原則：首先，在編程中要盡量避免采用循環(huán)語句。利用向量化語句來代替循環(huán)語句可以大大提高程序運行速度。如果不得不采用多重循環(huán)，那么內(nèi)循環(huán)的次數(shù)應(yīng)該盡可能多于外循環(huán)的次數(shù)。其次，在使用大型數(shù)組或矩陣之前對其進行初始化，即采用指令“zeros”或“ones”對矩陣定維，這樣可以減少MATLAB在內(nèi)存分配過程中的耗時，大大提高速度。 第三，應(yīng)該優(yōu)先考慮使用MATLAB的內(nèi)在函數(shù)。MATLAB的內(nèi)在函數(shù)是采用C語言優(yōu)化構(gòu)造的，并固化在MATLAB的內(nèi)核中，其運行速度可以和C語言的速度等價，而用戶自行編制的M文件則是在MATLAB中解釋執(zhí)行的。另外，采用更先進更有效的算法也可以提供計算效率。快速傅立葉變換就是一個典型的例子。2.2MATLAB編程仿真 2.2.1時間連續(xù)信號與系統(tǒng)的計算機仿真問題 時間連續(xù)的確定信號在物理上是一個隨時間變化的（電壓或電流）波形，在數(shù)學上表示為一個時間連續(xù)的函數(shù)f(t)。時間連續(xù)信號也稱為模擬信號。而時間離散的確定信號在數(shù)學上可以表示為一個確定的序列{fn}。如果以滿足取樣定理的取樣速率對時間連續(xù)函數(shù)f(t)進行取樣，就可以得到對應(yīng)的時間離散序列{fn}，而將時間離散序列{fn}通過取樣定理所規(guī)定的理想低通濾波器后，就能夠恢復(fù)連續(xù)信號f(t)。 換句話說，在滿足取樣定理的前提下，f(t)與{fn}是一一對應(yīng)的關(guān)系。對時間離散序列{fn}進行幅度值上的離散化（這個過程稱為量化）就得到了數(shù)字信號{}。數(shù)字信號通常以二進制編碼的形式存儲在計算機中。因此，利用計算機所表達的時間連續(xù)信號f(t)實質(zhì)上是其所對應(yīng)的數(shù)字信號{}。本質(zhì)上，計算機不能直接處理模擬信號，計算機中的信號處理均是對數(shù)字信號的處理。而由取樣定理保證了它與模擬信號處理之間的一一對應(yīng)關(guān)系。 連續(xù)系統(tǒng)可以用微分方程來描述，離散系統(tǒng)可以用差分方程來描述。當系統(tǒng)為無記憶系統(tǒng)時，微分方程或差分方程退化為代數(shù)方程。因此，對系統(tǒng)進行仿真的過程就數(shù)學意義而言，就是求解這些微分方程、差分方程或代數(shù)方程的過程。 為了對連續(xù)系統(tǒng)進行仿真，首先需要建立其數(shù)學模型，然后利用計算機求解這些數(shù)學模型，從而得出數(shù)學模型的數(shù)值解。計算機對數(shù)學模型的求解過程就是對系統(tǒng)的仿真過程。由于連續(xù)系統(tǒng)是通過微分方程來建模的，因此對此的計算機仿真本質(zhì)上就是微分方程的數(shù)值求解問題。通常，對連續(xù)系統(tǒng)進行計算機仿真的過程是：

（1)建立數(shù)學模型：就是將物理模型轉(zhuǎn)換為數(shù)學模型的過程。對于電子電路而言，就是利用電路的拓撲、元件的屬性，列出網(wǎng)孔方程或節(jié)點方程，并簡化為標準形式的計算機可求解的一階微分方程組（即狀態(tài)方程）的過程。對于線性時不變系統(tǒng)來說，就是建立其傳遞函數(shù)模型，并簡化為計算機易于求解的線性一階微分方程組（即狀態(tài)方程）的過程。而對于一般的系統(tǒng)，通常可建模為一組非線性的微分方程組。 （2)選擇適合的計算機求解方法求解仿真模型。不同的微分方程求解算法的精度和適用范圍均有所不同。 （3)編寫MATLAB仿真程序或建立Simulink模塊方框圖，調(diào)試并運行程序。（4)得出數(shù)值解，即仿真結(jié)果，對仿真結(jié)果進行分析，以確定結(jié)果的可靠性和有效性。圖2-2一個二階電路系統(tǒng) 下面我們舉例說明。對于一個二階系統(tǒng)，電路如圖2-2所示，試利用MATLAB編程仿真求解：

（1)零狀態(tài)響應(yīng)：系統(tǒng)在t=-1秒時刻的初狀態(tài)為零,輸入信號為階躍信號f(t)=u(t)，輸出為電容電壓，試對時間t=［-1,10］秒的輸出信號波形進行仿真求解。（2)零輸入響應(yīng)：系統(tǒng)在t=-1秒時刻電容電壓為uC=1V，當t=-1秒時刻信號輸入二端閉合，試對時間t=［-1,10］秒的輸出信號波形進行仿真求解。（3)全響應(yīng)：系統(tǒng)在t=-1秒時刻電容電壓為uC=1V，輸入信號為階躍信號f(t)=u(t)，試對時間t=［-1,10］秒的輸出信號波形進行仿真求解。 ［解］（1)設(shè)電感電流為iL(t)，電容電壓為uC(t)，根據(jù)電路，列出KVL方程：

(2-2) 改寫為標準形式：

令電路的狀態(tài)變量iL(t)=x1(t)，uC(t)=x2(t)，輸入信號為f(t),則微分方程簡寫為 寫出該微分方程的MATLAB函數(shù)。注意，這一組微分方程組就是該二階電路的狀態(tài)方程，對于同一個物理系統(tǒng)，如果選取的狀態(tài)變量不同，所得到的狀態(tài)方程也就不同，但都表達了該物理系統(tǒng)的狀態(tài)信息。(2-3) 按照MATLAB的固定寫法，編寫出該微分方程組的MATLAB函數(shù)（函數(shù)文件名為funcforex123.m）如下： 程序2-3

functionxdot=funcforex123(t,x,flag,R,L,C)

xdot=zeros(2,1);%矩陣初始化

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);%方程1

xdot(2)=1/C*x(1);%方程2

functionin=f(t)%輸入信號

in=(t>0)*1;%階躍信號

然后，利用MATLAB提供的求解微分方程的指令對該微分方程組求解。MATLAB提供的求解微分方程的算法有多個，如“ode45”、“ode23”、“ode15s”等，不同的算法適用的場合稍有不同。例如，通過“ode45”函數(shù)求解，MATLAB程序（程序名為ex123.m）如下： 程序2-4

%filenameex123.m L=1;%電感值

C=0.1;%電容值

forR=［1.535］%仿真電阻值分別為1.5，3，5歐姆的情況 ［t,x］=ode45(′funcforex123′,［-1,10］,［0;0］,［］,R,L,C);

%也可采用ode23，ode15s等求解

figure(1);plot(t,x(:,1));holdon;xlabel(′timesec′); text(0.9,0.17,′＼leftarrowi-L(t)′);grid; figure(2);plot(t,x(:,2));holdon;xlabel(′timesec′); text(0.5,0.3,′＼leftarrowu-C(t)′);grid; end

運行程序后，得到的電感電流、電容電壓波形仿真結(jié)果如圖2-3（a）、(b）所示。圖2-3電感電流、電容電壓零狀態(tài)響應(yīng)波形仿真結(jié)果（a)電感電流單位階躍響應(yīng)（零狀態(tài)響應(yīng)）波形仿真結(jié)果；（b)電容電壓單位階躍響應(yīng)（零狀態(tài)響應(yīng)）波形仿真結(jié)果 （2)當t=-1秒時刻信號輸入二端閉合，相當于輸入信號為零，因此修改MATLAB函數(shù)funcforex123.m如下： 程序2-5

functionxdot=funcforex123(t,x,flag,R,L,C)

xdot=zeros(2,1);%矩陣初始化

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);%方程1

xdot(2)=1/C*x(1);%方程2

functionin=f(t)%輸入信號

in=0;%信號始終為零 并修改ex123.m的“ode45”語句，將系統(tǒng)狀態(tài)改為［0，1］，程序如下： 程序2-6 %filenameex123.m L=1;%電感值

C=0.1;%電容值

forR=［1.535］%仿真電阻值分別為1.5，3，5歐姆的情況［t,x］=ode45(′funcforex123′,［-1,10］,［0;1］,［］,R,L,C); %也可采用ode23，ode15s等求解

figure(1);plot(t,x(:,1));holdon;xlabel(′timesec′);

text(0.9,0.07,′＼leftarrowi-L(t)′);grid; figure(2);plot(t,x(:,2));holdon;xlabel(′timesec′); text(0.5,0.3,′＼leftarrowu-C(t)′);grid; end

然后運行，得到的波形仿真結(jié)果如圖2-4所示。圖2-4電感電流、電容電壓零輸入響應(yīng)波形仿真結(jié)果（a)電感電流零輸入響應(yīng)波形仿真結(jié)果；（b)電容電壓零輸入響應(yīng)波形仿真結(jié)果

（3)對于全響應(yīng)情況，可相應(yīng)修改MATLAB函數(shù)funcforex123.m如下： 程序2-7

functionxdot=funcforex123(t,x,flag,R,L,C)

xdot=zeros(2,1); %矩陣初始化

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);%方程1

xdot(2)=1/C*x(1);%方程2

functionin=f(t)%輸入信號

in=(t>0)*1;%信號階躍 而ex123.m的“ode45”語句中仍然將系統(tǒng)狀態(tài)改為［0，1］。運行后得到的波形仿真結(jié)果如圖2-5所示。圖2-5電感電流、電容電壓全響應(yīng)波形仿真結(jié)果（a)電感電流全響應(yīng)波形仿真結(jié)果；（b)電容電壓全響應(yīng)波形仿真結(jié)果 由此可見，一旦從物理系統(tǒng)中建立了其數(shù)學模型，也即得出系統(tǒng)的微分方程描述之后，給定其邊界條件，就可得出其數(shù)值解，求解的過程就是仿真的過程。在MATLAB中，我們可以方便地修改微分方程函數(shù)以及求解指令的參數(shù)，從而得出系統(tǒng)的零輸入響應(yīng)、零狀態(tài)響應(yīng)以及全響應(yīng)等各種狀態(tài)下的輸出信號。讀者可以將這個例子中的輸入信號修改為一個寬度為0.3s的矩形脈沖，然后觀察其響應(yīng)波形。(提示：只需修改funcforex123.m中語句“in=(t>0)*1;”為“in=(t>0)*1-(t>0.3)*1；”即可。） 本例演示了對于一個動態(tài)電路瞬態(tài)響應(yīng)的仿真過程。仿真的關(guān)鍵在于建立正確的數(shù)學模型。 對于任何電子電路，都可以根據(jù)電路理論的知識得出其數(shù)學模型,即一階微分方程組，也稱為狀態(tài)方程。接下來以MATLAB所規(guī)定的格式編寫狀態(tài)方程的描述函數(shù)，然后選擇合適的求解方式求解，也就完成了對電路的仿真。最后，對所得出的仿真波形進行分析。 需要強調(diào)的是，微分方程的MATLAB函數(shù)格式是固定的，必須按MATLAB所規(guī)定的格式來編寫。本例中，需要將方程參數(shù)“R;L;C”傳入。傳入?yún)?shù)前的標志變量“flag”是必需的，以符合求解函數(shù)“ode45”的要求。微分方程的MATLAB函數(shù)引導語句的格式為

functionxdot=方程函數(shù)名（t,x,flag,附加參數(shù)） 其中，“t”為時間變量矩陣，“x”為方程的狀態(tài)變量矩陣，“xdot”為狀態(tài)變量對時間的一階導數(shù)矩陣。在MATLAB中，對微分方程的求解有多種算法可供利用。通常有“ode23”、“ode45”以及“ode15s”等等，這些函數(shù)的調(diào)用方式完全相同，但針對的方程類型有所不同。其中，“ode23”和“ode45”分別采用二階三級RKF算法和四階五級RKF算法，適用于一般微分方程的求解；“ode15s”則適用于剛性微分方程的求解。所謂剛性方程，指的是方程的某些解（這些解是時間“t”的函數(shù)）變化緩慢，而另外一些解變化快速，兩者形成明顯對比的這樣一類方程?！皁de15s”同時也適合于一般微分方程的求解。關(guān)于其它解法和使用方法詳見MTALAB的在線手冊。 以“ode45”為例，其用法是：

［t,x］=ode45(方程函數(shù)名，仿真時間范圍，狀態(tài)變量初始值，算法選項，附加參數(shù))

例如，以默認算法選項計算微分方程組funcforex123.m，仿真時間從-1s到10s，狀態(tài)變量初始值為：電感電流x1(0)=0A，電容電壓x2(0)=1V，附加 參數(shù)為R=1.5Ω，L=1H，C=0.1F,則解法寫為 ［t,x］=ode45(′funcforex123′,［-1,10］,［0;1］,［］,1.5,1,0.1);

另外需要注意的是，為了達到算法快速有效，MATLAB中的狀態(tài)方程的求解函數(shù)是自適應(yīng)步長的，即其數(shù)值結(jié)果中時間自變量t的離散取值是非均勻的，在解變化劇烈的地方取值點密，而在解變化緩慢的曲線段相對取值會稀疏一些。在信號處理中，進一步對求解信號進行分析時，往往又需要在時間上均勻取樣間隔點上的波形函數(shù)值。我們可以采用一維插值的辦法來解決。讀者可以閱讀插值命令“interp1”的在線幫助進一步了解這方面的信息。 下面介紹求解系統(tǒng)的沖激響應(yīng)的方法。由于輸入沖激信號是超越函數(shù)，不適于作數(shù)值計算，因此首先求系統(tǒng)的階躍響應(yīng)，然后對其進行微分得到?jīng)_激響應(yīng)。為了應(yīng)用差分函數(shù)“diff”進行近似微分，首先需要通過插值將被微分的數(shù)據(jù)在時間上均勻化，然后再利用差分函數(shù)“diff”進行近似微分，從而得出系統(tǒng)的沖激響應(yīng)波形。 例如，對于電路如圖2-2所示的一個二階系統(tǒng)，輸出信號從電容電壓取出?，F(xiàn)在通過仿真計算求系統(tǒng)的沖激響應(yīng)，并與理論值相對比,分析仿真結(jié)果的正確性和仿真精度。 為了分析仿真結(jié)果的正確性和仿真精度，我們首先計算電路中電容電壓的沖激響應(yīng)的解析值。將狀態(tài)方程寫為二階微分方程形式，有(2-4) 進行拉氏變換，得出傳遞函數(shù)為 將R=3Ω,L=1H,C=0.1F代入，并將H(s)化簡為 利用拉氏變換對(2-5)(2-6)(2-7) 得出沖激響應(yīng)為 (2-8) 利用下面的程序求出上例電路中電容電壓的均勻取樣點的輸出波形值（階躍響應(yīng)），然后對該波形進行微分得出沖激響應(yīng)，并與理論值進行比較。程序運行結(jié)果如圖2-6所示。圖2-6沖激響應(yīng)的仿真結(jié)果與理論結(jié)果的比較 （1)狀態(tài)方程函數(shù)程序funcforex123.m如下： 程序2-8

functionxdot=funcforex123(t,x,flag,R,L,C)

xdot=zeros(2,1);%矩陣初始化

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);%方程1

xdot(2)=1/C*x(1);%方程2

functionin=f(t)%輸入信號是

in=(t>0)*1;%階躍信號 （2)主程序ex123b.m如下： 程序2-9 %filenameex123b.m ［t,x］=ode45(′funcforex123′,［-1,10］,［0;0］,［］,3,1,0.1); %求解上例電路的階躍響應(yīng)

ts=0.001; t1=-1:ts:10; %均勻取樣點間隔0.001s x1=interp1(t,x(:,2),t1,′spline′);%插值得出均勻時間樣點的電容電壓波形值;

plot(t1,x1,′k-.′);%階躍響應(yīng)波形繪制

holdon; x1dot=［diff(［x1］)/ts,0］;%近似微分得出沖激響應(yīng)波形

plot(t1,［diff(［x1］)/ts,0］,′k′);xlabel(′timesec′); %沖激響應(yīng)波形繪制

ht=10/sqrt(7.75).*exp(-1.5*t1).*sin(sqrt(7.75)*t1).*(t1>0); %沖激響應(yīng)的理論值

plot(t1(1:50:length(t1)),ht(1:50:length(t1)),′ko′); legend(′Stepresponse′,′Impluseresponse′,′Theoreticimpluseresponse′); 從圖2-6中可以得出，在插值精度達到足夠高的條件下，采用仿真數(shù)值計算方法所得出的系統(tǒng)沖激響應(yīng)與采用理論解析公式計算的曲線是吻合的，這也從一個側(cè)面說明了仿真數(shù)值計算的正確性。 通過以上仿真的實例，我們可總結(jié)出連續(xù)系統(tǒng)建模的一般規(guī)律。一般地，設(shè)線性時不變系統(tǒng)中有n個獨立的儲能元件（電容或電感），則系統(tǒng)中存在n個獨立的狀態(tài)，記為

x1(t),x2(t),…,xn(t)(2-9) 又設(shè)該系統(tǒng)有m個輸入信號，記為

f1(t),f2(t),…,fm(t)(2-10)

系統(tǒng)有r個輸出，記為

y1(t),y2(t),…,yr(t)(2-11)

由于所有輸入輸出信號、系統(tǒng)狀態(tài)都是時間t的函數(shù)，簡寫時可以省略(t)，于是狀態(tài)方程和輸出方程可表示為矩陣形式：(2-12) 對于更一般的（非線性）時不變系統(tǒng)，狀態(tài)方程和輸出方程可表示為一般函數(shù)式： 對于線性時不變系統(tǒng)，狀態(tài)方程為n元一階線性微分方程組，輸出方程是r元線性方程組,其中(2-13)(2-14)(2-15)(2-16)(2-17)(2-18)(2-19) 對連續(xù)系統(tǒng)的計算機仿真，實質(zhì)上就是利用求解常微分方程的數(shù)值算法求解的過程。那么，計算機是如何求解這樣一個一階常微分方程組的呢？ 不失一般性，設(shè)常微分方程組用矩陣形式表示為 式中(2-20)(2-21)(2-22)(2-23) 式,a≤t≤b的解x是時間t的連續(xù)函數(shù)，數(shù)值解法就是求解x在若干離散時間點a=t0<t1<…<tn=b的值

x(t0),x(t1),…,x(tn)。對式,a≤t≤b兩邊在某一區(qū)間［tk,tk+1］上積分，設(shè)區(qū)間長度為h=tk+1-tk，有 即

(2-24)(2-25) 移項并利用矩形法對積分近似，得 上式簡記為(2-26)(2-27) 利用初值x(a)=x(t0)=x0，由上式遞推，就可以求出各離散時間t1,t2,…,tn的狀態(tài)值xk。這就稱為歐拉算法。 由于歐拉算法利用了矩形法近似求解積分，計算精度較低。為了提高計算精度，可以減小步長h，但是這樣將增加計算量，減小計算速度，再考慮到計算機字長有限所引入的舍入誤差，通過減小計算步長來提高計算精度的方法效果是有限的。 根據(jù)上述原理，可以編制出歐拉算法的MATLAB程序sybeuler.m。

程序2-10

%filenamesybeuler.m %sybeuler實現(xiàn)了歐拉定步長微分方程求解 % %［tout,yout］=sybeuler(odefile,t0,h,th,y0,P) % %odefile為描述微分方程的M函數(shù)，t0,h和th分別為 %初始時間、步長和終止時間，y0為狀態(tài)變量初值列向量，返回 %的tout和yout分別為時間向量和狀態(tài)變量解矩陣，P為 %微分方程的參數(shù)矩陣。

與odefile函數(shù)的參數(shù)設(shè)計相匹配

function［tout,yout］=sybeuler(odefile,t0,h,th,y0,P) tout=［t0:h:th］′; %仿真時間段以及仿真步長

yout(length(tout),length(y0))=0; %分配并初始化存儲空間，以加快運算速度

kk=1;%輸出矩陣腳標變量

fort=tout′

yout(kk,:)=y0′;%保存當前計算值

kk=kk+1; k1=h*eval(［odefile′(t,y0,P)′］); y0=y0+k1;%計算歐拉公式

end 為便于結(jié)果比較，仍然采用圖2-2所示的電路作為計算對象。其狀態(tài)方程的MATLAB函數(shù)如下（mystateEQ.m）：

程序2-11

functionxdot=mystateEQ(t,x,P)

xdot=zeros(2,1); R=P(1);L=P(2);C=P(3);%RLC參數(shù)從P矩陣中讀入

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);

xdot(2)=1/C*x(1);

functionin=f(t)%輸入信號

in=(t>0)*1;%階躍信號

作為對比，調(diào)用“ode45”進行求解的狀態(tài)方程的MATLAB函數(shù)mystateEQforODE45.m如下： （注意，僅在函數(shù)變量中加上了“flag”以適應(yīng)“ode45”的調(diào)用規(guī) 程序2-12

functionxdot=mystateEQforODE45(t,x,flag,P)

xdot=zeros(2,1); R=P(1);L=P(2);C=P(3);%RLC參數(shù)從P矩陣中讀入

xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);

xdot(2)=1/C*x(1);

functionin=f(t)%輸入信號是

in=(t>0)*1;%階躍信號

對歐拉算法與“ode45”算法精度進行比較的主程序如下：(其中，歐拉算法采用步長分別為0.01，0.1，0.2。“ode45”算法以默認參數(shù)進行) 程序2-13

%filenamesybeulerTEST.m tic;%記時開始 ［t,x］=ode45(′mystateEQforODE45′,［0,10］,［0;0］,［］,［310.1］);

toc%記時結(jié)束，測定ode45計算耗時

plot(t,x,′k.′);holdon;

forh=［0.2,0.1,0.01］

tic; ［t,x］=sybeuler(′mystateEQ′,0,h,10,［0;0］,［310.1］); toc%記時結(jié)束，測定sybeuler計算耗時

plot(t,x,′k′); end legend(′u-C(t)ode45′,′i- L(t)ode45′,′sybeuler′);

程序運行結(jié)果如圖2-7所示。圖2-7歐拉算法與“ode45”算法的比較(歐拉算法精度隨著計算步長減小而提高) 計算耗時： “ode45”算法自適應(yīng)步長0.1s

歐拉算法步長0.2[DW]0.02s

歐拉算法步長0.1[DW]0.04s

歐拉算法步長0.01[DW]0.32s

可見,歐拉算法的計算精度是比較低的，當歐拉算法步長達到0.01時，其計算曲線才基本與“ode45”算法的曲線吻合。同時我們也看到，數(shù)值求解算法以及算法參數(shù)的選擇在系統(tǒng)仿真中對仿真結(jié)果的影響是相當大的。

MATLAB中的“ode45”算法是四階龍格-庫塔算法的改進算法，稱為四階五級的龍格-庫塔算法。與傳統(tǒng)的四階龍格-庫塔算法相比較，它具有更高的計算精度和穩(wěn)定性。 對于線性時不變連續(xù)時間系統(tǒng)，其數(shù)學模型為線性微分方程組（也即狀態(tài)方程），這種方程組在數(shù)學上已經(jīng)有了解析算法，可以得出解析解。然而，對于一般的連續(xù)時間信號與系統(tǒng)，其數(shù)學模型為一般微分方程組，通?？赡苁欠蔷€性的，對于非線性微分方程組的解析求解，目前數(shù)學上并無統(tǒng)一的方法，這時大量的要求求解非線性微分方程組的工程問題最為方便的辦法就是數(shù)值求解了。可以說，MATLAB仿真和數(shù)值計算擴展了傳統(tǒng)線性系統(tǒng)的研究范圍，為非線性系統(tǒng)的數(shù)值研究提供了一條捷徑。 總之，對連續(xù)時間信號與系統(tǒng)的計算機仿真問題就是建立連續(xù)系統(tǒng)的數(shù)學模型——狀態(tài)方程，然后利用適當?shù)臄?shù)值分析和求解算法計算其數(shù)值解的過程。類似地，對離散時間信號與系統(tǒng)的計算機仿真問題就是建立該離散時間系統(tǒng)的差分方程組（離散狀態(tài)方程）并遞推求解的過程。而對于混合系統(tǒng)的計算機仿真問題就成為了混合系統(tǒng)所對應(yīng)的微分方程和差分方程聯(lián)合方程組的求解問題。 對于無記憶系統(tǒng)，以上系統(tǒng)的建模方程組退化為代數(shù)方程組，仿真的問題也就相應(yīng)簡化為對這些代數(shù)方程的數(shù)值計算問題。

2.2.2基于數(shù)據(jù)流的仿真方法

MATLAB編程對數(shù)學方程的數(shù)值計算一般是基于數(shù)據(jù)流的，這稱為基于數(shù)據(jù)流的仿真方法。 所謂基于數(shù)據(jù)流的仿真，就是指在數(shù)值計算的一個階段完成之后，所得出的數(shù)據(jù)才進入下一個階段的數(shù)值計算中。例如，如果對由兩個傳遞函數(shù)H1(s)和H2(s)級聯(lián)的系統(tǒng)進行MATLAB編程仿真，第一步是求解在輸入信號激勵下H1(s)的輸出波形數(shù)據(jù)，第二步，再利用這些得出的波形數(shù)據(jù)作為系統(tǒng)H2(s)的激勵信號，去求解系統(tǒng)H2(s)的輸出波形數(shù)據(jù)?？梢姡谇蠼獾臅r間段上，我們無法“實時”地觀察到系統(tǒng)各個輸出點上的信號，只有當系統(tǒng)各個階段均已經(jīng)計算完畢之后，才能夠觀察到計算時間段上系統(tǒng)的狀態(tài)和各個輸出點的變化情況。2.3Simulink仿真 2.3.1Simulink仿真入門 1.Simulink的特點

Simulink是MATLAB中的一個建立系統(tǒng)方框圖和基于方框圖級的系統(tǒng)仿真環(huán)境，是一個對動態(tài)系統(tǒng)進行建模、仿真并對仿真結(jié)果進行分析的軟件包。使用Simulink可以更加方便地對系統(tǒng)進行可視化建模，并進行基于時間流的系統(tǒng)級仿真，使得仿真系統(tǒng)建模與工程中的方框圖統(tǒng)一起來。并且仿真結(jié)果可以近乎“實時”地通過可視化模塊，如示波器模塊、頻譜儀模塊以及數(shù)據(jù)輸入輸出模塊等顯示出來，使得系統(tǒng)仿真工作大為方便。

Simulink使得用戶可以用鼠標操作將一系列可視化模塊連接起來，從而建立直觀的功能上更為復(fù)雜的系統(tǒng)模型，避免了編寫MATLAB仿真程序，簡化了仿真建模過程，更加適用于大型系統(tǒng)的建模和仿真，如對IS-95CDMA通信系統(tǒng)全系統(tǒng)的建模仿真工作。 通過上節(jié)的討論我們知道，由于MATLAB具有微分方程的求解算法，利用MATLAB編程可以對微分方程進行求解，從而得到系統(tǒng)的仿真結(jié)果。然而，對于更為復(fù)雜的大系統(tǒng)，往往是一種連續(xù)、離散系統(tǒng)的混合，要通過編程語句的方式來建立整個系統(tǒng)的狀態(tài)方程模型將是一件不容易的事情，況且這樣的手工編程，直觀性不好，對于復(fù)雜系統(tǒng)的建模編程來說往往錯誤難以避免，因此仿真結(jié)果的可信度也就成了問題，而且手工編程的效率低下，出錯（往往是算法上的錯誤）檢查困難，編程程序的可重復(fù)利用程度不高。

這些原因使得MATLAB向易用性和可視化方向努力，Simulink就是這種努力的結(jié)果。1990年，Simulink首次作為MATLAB的軟件工具包出現(xiàn)，徹底改變了系統(tǒng)仿真界的軟件工具和建模模式，很快成為了系統(tǒng)建模和仿真的主流軟件工具，甚至將成為通用系統(tǒng)仿真的基準測試平臺。

Simulink仿真環(huán)境附帶了許多專業(yè)仿真模塊庫，利用這些模塊庫可以快速建立該專業(yè)領(lǐng)域的系統(tǒng)模型并進行仿真，而不需要用戶詳細了解各個模塊內(nèi)部的實現(xiàn)細節(jié)，大大方便了復(fù)雜的大系統(tǒng)的建模，而且，由于Simulink提供的這些專業(yè)模塊庫均通過各專業(yè)權(quán)威專家的評測，可信度和穩(wěn)定性都大為提高，從而也保證了系統(tǒng)建模整體的質(zhì)量和仿真精度。在通信工程和電子工程領(lǐng)域，Simulink提供的常用專業(yè)模塊庫有：CDMA參考模塊庫、通信系統(tǒng)模塊庫、DSP（數(shù)字信號處理器）模塊庫等等。隨著MATLAB版本升級，還將添加更多的專業(yè)模塊庫。

Simulink全方位地支持動態(tài)系統(tǒng)的建模仿真，它支持連續(xù)系統(tǒng)、離散系統(tǒng)、連續(xù)離散混合系統(tǒng)、線性系統(tǒng)、非線性系統(tǒng)、時不變系統(tǒng)、時變系統(tǒng)的建模仿真，也支持具有多采樣速率的多速率系統(tǒng)?？梢哉f，在通用系統(tǒng)仿真領(lǐng)域，Simulink是無所不包的。結(jié)合MATLAB編程和Simulink可視化建模仿真各自的特點，可以構(gòu)建更為復(fù)雜的系統(tǒng)模型，并進行自動化程度更高的仿真和仿真結(jié)果的數(shù)據(jù)分析，這是MATLAB的高級應(yīng)用方面。本書后面章節(jié)中的許多例子就是基于MATLAB編程和Simulink建模相融合的方式來構(gòu)建的。 2.使用Simulink建模和仿真的過程 我們首先通過實例來簡單介紹Simulink建模和仿真的全過程。啟動MATLAB之后，在命令窗口中輸入命令“simulink”或單擊MATLAB工具欄上的Simulink圖標，打開Simulink

模塊庫窗口（使用命令“simulink3”可以打開老版本的Simulink模塊庫界面）。 在Simulink模塊庫窗口中單擊菜單項“File|New|Model”，就可以新建一個Simulink模型文件，如圖2-8所示。圖2-8Simulink模塊庫界面和新建模型文件窗口 為了與MATLAB編程仿真進行對比，我們?nèi)匀灰?.3.1小節(jié)描述的二階系統(tǒng)為例。其輸入輸出的傳遞函數(shù)為 (2-28)

其中，R=3Ω,L=1H,C=0.1F。H(s)簡化為 利用鼠標單擊Simulink基礎(chǔ)庫中的Continuous子庫，選取傳遞函數(shù)模塊，將它拖動到新建模型窗口中的適當位置，如圖2-9所示。

(2-29)圖2-9利用模塊庫建立仿真模型 如果需要對模型模塊進行參數(shù)設(shè)置和修改，只需選中模型文件中的相應(yīng)模塊，單擊鼠標右鍵，彈出快捷菜單，從中選取相應(yīng)參數(shù)進行修改，如圖2-10所示。還可以在選中模塊之后通過鼠標拖動修改模塊的位置、大小和形狀。單擊模塊下方的“TransferFun”可以對其進行編輯，例如修改為“傳遞函數(shù)”字樣。從快捷菜單中選取“TransferFunparameters...”項修改傳遞函數(shù)參數(shù)，在彈出的對話框中的傳遞函數(shù)分子系數(shù)“Numerator：”欄填入［10］；在傳遞函數(shù)分母系數(shù)“Denominator”欄填入［1，3，10］，其余參數(shù)使用默認值。如果需要進了解該模塊的參數(shù)設(shè)置說明，可以單擊該對話框下方的“Help”按鈕，然后確認,就得到了所要進行仿真的傳遞函數(shù)。圖2-10修改仿真模型的參數(shù) 通過快捷菜單的其它選項還可以修改模型的顏色、旋轉(zhuǎn)、字體、陰影等屬性，也可對模型進行剪切、拷貝或刪除。

采用同樣的方法，在Simulink基礎(chǔ)庫中的Sources子庫中選取激勵信號源，例如我們選取階躍信號源，將之拖入建模窗口中。在Sinks子庫中選取示波器作為系統(tǒng)輸出波形顯示。接下來利用鼠標將這三個模塊連接起來。模塊外部的大于符號“>”分別表示信號的輸入輸出節(jié)點，為了連接兩個模塊的輸入輸出，可以將鼠標置于節(jié)點處，這時鼠標顯示為“十”字形狀，拖動鼠標到另一個模塊的端口，然后釋放鼠標按鈕，則形成了帶箭頭的連線，箭頭方向表示信號的流向。完成后的建模系統(tǒng)可以通過“File”菜單存盤為模型文件，擴展名為“mdl”，如“l(fā)izi1.mdl”，如圖2-11所示。圖2-11完成的建模方框圖 接下來，需要對輸入信號源（階躍）的參數(shù)進行設(shè)置。將鼠標指向階躍信號模塊雙擊或通過快捷菜單打開屬性設(shè)置對話框，設(shè)置階躍信號的參數(shù)，如圖2-12所示。圖中右邊的幫助窗口是通過單擊參數(shù)設(shè)置對話框下方的“Help”按鈕顯示的。通過閱讀幫助文檔可以了解參數(shù)的含義和設(shè)置情況。對于階躍信號源來說，其參數(shù)含義、默認值以及我們根據(jù)仿真需要修改后的參數(shù)值如表2-1所示。表2-1階躍信號源的參數(shù)含義、默認值及修改后的參數(shù)值圖2-12階躍信號源模塊的參數(shù)設(shè)置對話框及其幫助文檔 根據(jù)仿真要求，我們將仿真-2~10s時間區(qū)間內(nèi)的系統(tǒng)階躍響應(yīng)的輸出波形，在時間為0s時的躍變。 然后設(shè)置仿真參數(shù)，主要是仿真求解器的選擇和仿真步長等參數(shù)的選取。通過前面的分析可知，動態(tài)系統(tǒng)仿真的本質(zhì)就是求解其狀態(tài)方程，而對狀態(tài)方程的數(shù)值求解算法有多種，求解算法的步長也可以不同。不同的算法適用的范圍有所不同，而算法的步長也直接影響求解的精度。因此，對求解器的選擇以及其仿真步長等參數(shù)的設(shè)定對系統(tǒng)仿真來說就成為相當重要的事情。從系統(tǒng)建模窗口的狀態(tài)欄可以看到當前使用的求解器，如圖2-11中顯示的仿真求解器是“ode45”算法。從建模窗口菜單項“Simulation|SimulationParameters...”打開仿真參數(shù)設(shè)置對話框（快捷鍵為Ctrl+E），我們現(xiàn)在設(shè)置求解器標簽下的參數(shù)部分，如圖2-13所示。 設(shè)置仿真起始時間為-2秒，仿真結(jié)束時間為10秒，其余參數(shù)為默認值：求解器采用“ode45”算法，步長設(shè)定為自適應(yīng)變步長的，最大步長、最小步長以及初始步長均設(shè)為自動選取，相對求解精度為1e-3，絕對求解精度自動選取。圖2-13仿真參數(shù)設(shè)置對話框 最后，雙擊示波器模型圖標，打開示波器顯示窗口。在顯示窗口中單擊鼠標右鍵，通過快捷菜單設(shè)置顯示坐標范圍等屬性，這里我們設(shè)置為自動刻度，如圖2-14所示。圖2-14示波器顯示窗口以及參數(shù)設(shè)置 所有這些工作完成之后，就可以進行仿真了?？赏ㄟ^建模窗口菜單項“Simulation|Start”啟動仿真，也可以單擊工具欄上的小三角按鈕或使用快捷鍵Ctrl+T啟動仿真。仿真結(jié)果如圖2-15所示。讀者可以將結(jié)果與圖2-3所示的結(jié)果進行比較。圖2-15仿真結(jié)果

更換信號源為Sources子模塊庫中的SignalGenerator，并設(shè)置信號源為0.2Hz的方波，幅度為1，如圖2-16左邊對話框所示。設(shè)置示波器顯示窗口的屬性（Parameters），使之成為雙蹤顯示，然后將示波器第二輸入節(jié)點與信號源輸出相連，這樣我們就可以同時觀察系統(tǒng)的輸入輸出波形了。系統(tǒng)建模如圖2-16中間窗口所示。將仿真時間設(shè)定為0秒到20秒，其余參數(shù)使用默認參數(shù)。運行仿真后的結(jié)果如圖2-16右邊窗口顯示。讀者還可以進一步修改信號源參數(shù)，使用三角波、正弦波等作為激勵信號，觀察輸出信號的情況。

圖2-16更換信號源并使用雙蹤示波器之后的仿真結(jié)果

3.MATLAB命令窗口與Simulink之間的交互 從以上仿真過程中我們可以感受到Simulink的方便和快捷。事實上，我們還可以通過MATLAB命令來打開Simulink模型并進行仿真。 在MATLAB命令窗口中，使用“openlizi1.mdl”，然后使用“sim(′lizi1.mdl′)”就可以啟動對模型lizi1.mdl的仿真計算，從而實現(xiàn)Simulink仿真的自動化。而Simulink仿真的數(shù)據(jù)結(jié)果也可以送回MATLAB工作空間中作進一步數(shù)值分析。仍以上例加以說明。首先設(shè)置示波器，使得顯示波形數(shù)據(jù)能夠送回MATLAB工作空間。在示波器波形顯示窗口單擊“參數(shù)Parameters”工具圖標，打開顯示參數(shù)設(shè)置對話框。選中Datahistory標簽下的Savedatatoworkspace，并設(shè)傳遞變量名稱，例如設(shè)為ScopeData，格式選擇為Structurewithtime（帶時間的結(jié)構(gòu)型變量）。參數(shù)設(shè)置的情況如圖2-17所示。圖2-17示波器顯示窗口的參數(shù)設(shè)置

將建模存盤為lizi1.mdl，然后在MATLAB下運行如下命令（參見圖2-18）： >>clear;%工作空間初始化 >>openlizi1.mdl;%這時將看到建模模型文件被打開 >>sim(′lizi1.mdl′);%啟動模型仿真，顯示出仿真波形 >>whos

NameSizeBytesClas

s ScopeData1x13578structarray Grandtotalis307elementsusing3578bytes圖2-18通過命令啟動Simulink仿真

可見，仿真完成之后，工作空間中出現(xiàn)了“ScopeData”結(jié)構(gòu)變量，其中包含了示波器顯示的全部波形數(shù)據(jù)。通過“plot”命令可以作出這些數(shù)據(jù)對應(yīng)的波形，即

>>t=ScopeData.time;%仿真的時間變量 >>signal=ScopeData.signals; >>wave1=signal(1,1).values; >>wave2=signal(1,2).values; >>subplot(2,1,1);plot(t,wave1); >>subplot(2,1,2);plot(t,wave2);axis(［0,20,-2,2］);

作出的波形如圖2-19所示。讀者可以對比示波器上得到的波形（參見圖2-18）。我們在此說明了Simulink與MATLAB工作空間進行數(shù)據(jù)交換的一種方法。圖2-19示波器數(shù)據(jù)傳入工作空間后進行波形作圖顯示 事實上，MATLAB提供了許多途徑用于與Simulink的數(shù)據(jù)交互。通過Sources子模型庫中的FromWorkspace模塊可以從工作空間中讀入仿真所需要的輸入數(shù)據(jù)，而通過Sinks子模型庫中的ToWorkspace模塊可以將Simulink中產(chǎn)生的數(shù)據(jù)回送到工作空間，從而可以實現(xiàn)MATLAB編程與Simulink模型相結(jié)合的混合仿真，這樣就大大加強了仿真的功能。下面的例子中，我們將從工作空間中讀入輸入波形數(shù)據(jù)，通過Simulink建模仿真之后，將仿真結(jié)果再送回工作空間中進行處理。通過這種方式，我們可以通過MATLAB