1、2.5MATLAB程序MATLAB是解釋性語言，輸入一行語句后回車，就會立即執(zhí)行得出結果。如果要實現(xiàn)比較復雜的功能，單靠一條一條地在命令窗口中輸入指令執(zhí)行，效率是很低的。如何解決這個問題呢？為此MATLAB提供了擴展名為“.m”的文本文件，在文件中事先寫入一行行的MATLAB命令，存盤后從MATLAB的命令窗口調入執(zhí)行（類似于DOS下的批處理），這種文件稱為底稿文件或MATLAB腳本文件。,在文本編輯器中輸入語句，然后將文件存盤，例如命名文件名為“my1stprg.m”，保存于MATLAB的默認工作路徑中。然后回到MATLAB命令窗口，在提示符“”下鍵入文件名(可以省略擴展名)后回車，即可運行

2、程序得到結果。,2.5.1M文件簡介MATLAB除了如前所述的在命令窗口進行的直接交互的指令操作方式外，另外一種更為重要的工作方式就是M文件的編程工作方式。,M文件有兩種形式：一種是腳本文件，另一種是函數文件。,腳本文件（命令文件）實際上是一串指令的集合，與在命令窗口逐行執(zhí)行文件中的所有指令，其結果是一樣的。如果M文件的第一行是以關鍵字“function”開頭的，則就是函數文件。,M文件,二者在內存使用上存在重要的區(qū)別：,函數文件內部所定義的變量僅僅在該函數文件內部有效，函數返回后這些內部變量將自動被清除，也就是說，函數內部所定義的變量僅在函數內部起作用，是局部變量。,而腳本文件中所定義的變量

3、，在使用“clear”命令清除之前，始終存在于工作空間當中，是全局變量。,注意:不要和MATLAB已經定義的系統(tǒng)函數和其它自定義的函數同名。,另外，函數文件的文件名必須和函數名相同.,（1）函數定義行（關鍵字function）functionout1,out2,.=filename(in1,in2,.)輸入和輸出（返回）的參數個數分別由nargin和nargout兩個MATLAB保留的變量來給出。（2）第一行幫助行，即H1行以（%）開頭，作為lookfor指令搜索的行（3）函數體說明及有關注解以（%）開頭，用以說明函數的作用及有關內容（4）函數體語句函數體內使用的除返回和輸入變量這些在func

4、tion語句中直接引用的變量以外的所有變量都是局部變量，即在該函數返回之后，這些變量會自動在MATLAB的工作空間中清除掉。如果希望這些中間變量成為在整個程序中都起作用的變量，則可以將它們設置為全局變量。,函數M文件,1.M文件的建立M文件是一個文本文件，它可以用任何編輯程序來建立和編輯，而一般常用且最為方便的是使用MATLAB提供的文本編輯器,有3種方法：(1)菜單操作。File菜單中選擇New菜單項，再選擇M-file命令。(2)命令操作。在MATLAB命令窗口輸入命令edit。(3)命令按鈕操作。單擊MATLAB主窗口工具欄上的NewM-File命令按鈕。,只是一個簡單的ASCII碼文本

5、文件，執(zhí)行程序時逐行解釋運行程序，matlab是解釋性的編程語言。,M文件的打開(1)菜單操作。從File菜單中選擇Open命令，在Open對話框中選中所需打開的M文件。(2)命令操作。在MATLAB命令窗口輸入命令：edit文件名。(3)命令按鈕操作。單擊MATLAB主窗口工具欄上的OpenFile命令按鈕，再從彈出的對話框中選擇所需打開的M文件。,3.數據和文件的輸入輸出1)inputinput指令提示用戶從鍵盤輸入數據、字符串或表達式，并接受該輸入。指令的調用格式有兩種：A=input(prompt)%輸入數據或表達式A=input(prompt,s)%輸入字符串例如：xm=input(

6、Whatsyourname?,s);,2)數據的輸出MATLAB提供的命令窗口輸出函數主要有disp函數，其調用格式為disp(輸出項)其中輸出項既可以為字符串，也可以為矩陣。,例輸入x,y的值，并將它們的值互換后輸出。程序如下：x=input(Inputxplease.);y=input(Inputyplease.);z=x;x=y;y=z;disp(x)disp(y),數據文件.matmat文件是matlab以標準二進制格式保存的數據文件，可將工作空間中有用的數據變量保存下來。mat文件的生成和調用是由函數save和load完成的。,(3)matlab文件的類型,m文件,m文件有兩類:獨立

7、的m文件稱命令文件(腳本文件)可調用m文件稱函數文件,4)利用文件輸入輸出數據指令“save”和“l(fā)oad”用于MATLAB與磁盤的數據交換?！皊ave”指令將MATLAB工作空間的數據存入磁盤；“l(fā)oad”指令則將數據從磁盤送入到MATLAB工作空間。舉例說明如下：A=1,2;3,4;b=81;save(Cmydate.mat,A,b)%將變量A和b保存在文件C:mydate.mat,load(Cmydate.mat,A,b);%調入數據文件whos%顯示調入的變量NameSizeBytesClassA2x232doublearrayb11x18doublearrayGrandtotalis

8、5elementsusing40bytes,2.5.2程序控制任何計算機語言，只要利用順序結構、循環(huán)結構以及分支結構，就可以完成任何程序功能。由于MATLAB語言矩陣計算功能十分強大，常常僅使用順序結構借以矩陣的邏輯運算就可以完成計算任務。由于循環(huán)結構和分支結構在MATLAB語言中的運行速度相對較慢，因此在算法優(yōu)化的編程中應當盡可能避免使用，而代之以矩陣運算，從而提高程序運行速度，簡化程序代碼，并使得程序代碼更加接近于數學表達。,1）.for循環(huán)語句基本格式for循環(huán)變量起始值：步長：終止值循環(huán)體end步長缺省值為1，可以在正實數或負實數范圍內任意指定。對于正數，循環(huán)變量的值大于終止值時，循環(huán)

9、結束；對于負數，循環(huán)變量的值小于終止值時，循環(huán)結束。循環(huán)結構可以嵌套使用。,1.循環(huán)結構,例如，最簡單的forend循環(huán)：fork=1:10x(k)=k.2;end事實上，采用矩陣思想也可以獲得相同的結果，但編程更加簡單明了：k=1:10;x=k.2,絕大部分循環(huán)都可以遵從這樣的方法變成向量化(矩陣化)的算法，避免采用循環(huán)語句，從而大大提高程序的執(zhí)行效率。,2）.while循環(huán)語句基本格式while表達式循環(huán)體end若表達式為真，則執(zhí)行循環(huán)體的內容，執(zhí)行后再判斷表達式是否為真，若不為真，則跳出循環(huán)體，向下繼續(xù)執(zhí)行。,While循環(huán)和for循環(huán)的區(qū)別在于，while循環(huán)結構的循環(huán)體被執(zhí)行的次數不

10、是確定的，而for結構中循環(huán)體的執(zhí)行次數是確定的。,與循環(huán)結構相關的語句還有break語句和continue語句,它們一般與if語句配合使用。break語句用于終止循環(huán)的執(zhí)行。當在循環(huán)體內執(zhí)行到該語句時，程序將跳出循環(huán)，繼續(xù)執(zhí)行循環(huán)語句的下一語句。continue語句控制跳過循環(huán)體中的某些語句。當在循環(huán)體內執(zhí)行到該語句時，程序將跳過循環(huán)體中所有剩下的語句，繼續(xù)下一次循環(huán)。,2.條件分支結構if，switch,if表達式語句1else語句2end,if表達式執(zhí)行語句end,(1)ifelseend語句有3種形式,if表達式1多分支語句1elseif表達式2語句2elseif表達式3語句3else

11、語句nend,這種結構實現(xiàn)了一種多路選擇，比較復雜。可替代一般高級語言中的switch-case-end語句,在條件表達式中，通常都是由關系操作符,=,=,=等邏輯操作符%清除MATLAB內存空間,這一命令常用于腳本文件的首句%以避免前面命令在內存空間形成的變量的意外影響F=inline(t,t);%用inline函數建立積分的被積函數forx=1,3,5;f_x(k)=x3+x+log(x)*sin(x)+quadl(F,0,x);%計算表達式值k=k+1;%數組下標加1end,f-x%顯示計算結果文件編輯存盤后,在MATLAB命令窗口執(zhí)行:,helpcalcfx這是開頭的注釋行,可以用he

12、lp命令顯示程序calcfx.m的功能是計算表達式calcfxf-x=2.500034.6550140.9567,下面是函數編程的例子，函數的文件名與函數名相同，為stat.m。程序2-2functionmean,stdev=stat(x)n=length(x);%求輸入向量x的長度mean=sum(x)/n;%求平均值stdev=sqrt(sum(x-mean).2/n);%求均方根值,m-函數實例:利用級數展開近似計算sin,sin的級數展開為,1:functions=powersin(x);2:%POWERSIN.Powerseriesforsin(x).3:%POWERSIN(x)tr

13、iestocomputesin(x)fromapowerseries,4:s=0;5:t=x;6:n=1:,7:whileabs(t)eps8:s=s+t;9:t=-x.2/(n+1)*(n+2).*t;10:n=n+2;11:end,利用“tic”和“toc”指令可以對程序段的執(zhí)行時間進行測定，從而估計出程序執(zhí)行效率，并找出改進程序、提高效率的方法?！皌ic”用于計時開始，而“toc”用于計時結束并顯示計時結果。,6.測定程序執(zhí)行時間和時間分配,7.提高程序執(zhí)行速度的原則要提高編程的執(zhí)行效率，則一定要根據MATLAB的特點來編寫程序，具體有以下原則：,其次，在使用大型數組或矩陣之前對其進行初

14、始化，即采用指令“zeros”或“ones”對矩陣定維，這樣可以減少MATLAB在內存分配過程中的耗時，大大提高速度。,首先，在編程中要盡量避免采用循環(huán)語句。利用向量化語句來代替循環(huán)語句可以大大提高程序運行速度。如果不得不采用多重循環(huán)，那么內循環(huán)的次數應該盡可能多于外循環(huán)的次數。,第三，應該優(yōu)先考慮使用MATLAB的內在函數。MATLAB的內在函數是采用C語言優(yōu)化構造的，并固化在MATLAB的內核中，其運行速度可以和C語言的速度等價，而用戶自行編制的M文件則是在MATLAB中解釋執(zhí)行的。另外，采用更先進更有效的算法也可以提供計算效率?？焖俑盗⑷~變換就是一個典型的例子。,2.5.4MATLAB編

15、程仿真,仿真也稱為模擬，在本質上，系統(tǒng)的計算機仿真就是根據物理系統(tǒng)的運行原理建立相應的數學描述并進行計算機數值求解的過程。,建模和仿真過程,建立系統(tǒng)數學模型或仿真模型。編寫系統(tǒng)的計算機仿真模型。求解。分析。,求解算法算法的改進簡化系統(tǒng)模型,關鍵問題：,1.系統(tǒng)仿真的數學基礎,2.計算機仿真的一般過程,建立計算機程序（層次化的）。執(zhí)行仿真。對仿真模型和仿真結果的檢驗。,對仿真數學模型有效性的驗證。對計算機仿真模型（程序）的驗證。對仿真算法的驗證。仿真結果置信度分析。,仿真驗證包含以下方面內容：,例1：試對空氣中在重力作用下不同質量物體的下落過程進行建模和仿真。已知重力加速度g=9.8m/s2，在

16、初始時刻t0=0s時物體由靜止開始墜落?？諝鈱β潴w的影響可以忽略不計。,質量為m的物體在自由墜落過程中受到豎直向下的恒定重力的作用，由牛頓第二定律，我們知道，重力F，加速度a以及物體質量m之間的關系是：,F=ma其中a=g,初始時刻為t0=0，物體的初始速度為v(t0)=0，并設物體下落的瞬時速度為v(t)。,（1）建立數學模型,設物體在t時刻的位移為s(t),并設初始位移為零，即s(t0)=0,根據加速度、速度、位移三者之間的微積分關系，我們得到一組數學方程：,a=dt/dvv=dt/dsF=ma,以及初始條件（也稱為方程的邊界條件）,v(t0)=0s(t0)=0,（2）數學模型的解析分析,

17、將方程轉換為一種在自變量（時間）上的“遞推”表達式,（3）根據數學模型建立計算機仿真模型（編程）,v(t+dt)=v(t)+dv=v(t)+adts(t+dt)=s(t)+ds=s(t)+v(t)dt,（4）執(zhí)行仿真和結果分析,%ch1example1prg1.mg=9.8;%重力加速度v=0;%設定初始速度條件s=0;%設定初始位移條件t=0;%設定起始時間dt=0.1;%設置計算步長N=20;%設置仿真遞推次數.仿真時間等于N與dt的乘積fork=1:Nv=v+g*dt;%計算新時刻的速度s(k+1)=s(k)+v*dt;%新位移t(k+1)=t(k)+dt;%時間更新end,%理論計算,

18、以便與仿真結果對照t_theory=0:0.01:N*dt;%設置解析計算的時間點v_theory=g*t_theory;%解析計算的瞬時速度s_theory=1/2*g*t_theory.2;%解析計算的瞬時位移%作圖:仿真結果與解析結果對比t=0:dt:N*dt;plot(t,s,o,t_theory,s_theory,-);xlabel(時間t);ylabel(位移s);legend(仿真結果,理論結果);,5)仿真程序的功能擴展-以動態(tài)方式來觀察物體墜落的過程,將作圖語句放在遞推計算循環(huán)內，并設置即時作圖刷新方式，從而得到這種“動畫”仿真的效果。,%ch1example1prg2.mg

19、=9.8;%重力加速度forL=1:5%仿真重復5次以便于觀察v=0;%初始速度s=0;%初始位置t=0;dt=0.01;%計算步長fork=1:200v=v+g*dt;%速度s=s+v*dt;%位移t=t+dt;%時間,plot(0,-s,o);axis(-22-200);%坐標范圍固定text(0.5,-1,當前時間：t=,num2str(t);text(0.5,-2,當前速度：v=,num2str(v);text(0.5,-3,當前位置：s=,num2str(s);set(gcf,DoubleBuffer,on);%雙緩沖避免作圖閃爍drawnow;%立即作圖endend,例2對乒乓球的

20、彈跳過程進行仿真。忽略空氣對球的影響，乒乓球垂直下落，落點為光滑的水平面，乒乓球接觸落點立即反彈。如果不考慮彈跳中的能量損耗，則反彈前后的瞬時速率不變，但方向相反。如果考慮撞擊損耗，則反彈速率有所降低。我們希望通過仿真得出乒乓球位移隨時間變化的關系曲線，并進行彈跳過程的“實時”動畫顯示。,設球是剛性的，質量為m，垂直下落。碰擊面為水平光滑平面。在理想情況下碰擊無能量損耗。如果考慮碰擊面損耗，則碰擊前后速度方向相反，大小按比例系數K(0K1)下降。在t時刻的速度設為v=v(t),位移設為y=y(t),并以碰擊點為坐標原點，水平方向為坐標橫軸建立直角坐標系。球體的速度以豎直向上方向為正方向。重力加

21、速度為g=9.8m/s2。,首先對乒乓球彈跳過程進行一些理想化假設。,設初始時刻t0=0球體的初始速度為v0=v(t0),初始位移為y0=y(t0)。,v(t+dt)=v(t)-gdt在t+dt時刻小球的位移為y(t+dt)=y(t)+v(t)dt,在小球撞擊水平面的瞬間，即y(t)=0的時刻，它的速度方向改變，大小按比例K衰減。當K=1時，就是無損耗彈跳情況。因此，小球反彈瞬間（t+dt時刻）的速度為,v(t+dt)=-Kv(t)-gdt;0K1,初始條件假設：,在空中時小球受重力F=mg作用，其中，g=-(dt/dv)。則在t+dt時刻小球的速度為:,受力分析：,反彈瞬間的位移為,y(t+

22、dt)=y(t)-Kv(t)dt=-Kv(t)dt,(2)仿真模型設計（程序）從數學模型中可見，小球在空中自由運動時刻與撞擊時刻的動力方程不同。通過小球所處位置（位移）是否為零可判定小球處于何種狀態(tài)。,%ch1example2prg1.mg=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球質量t0=0;%起始時間K=0.85;%彈跳的損耗系數N=5000;%仿真的總步進數dt=0.001;%仿真步長v=v0;%初狀態(tài)y=y0;,fork=1:Nif(y0)|(v0)%小球在空中的(含剛剛彈起瞬間)動力方程計算v=v-g*dt;y=y+v*dt;else%碰擊瞬間的計算

23、y=y-K.*v*dt;v=-K.*v-g*dt;ends(k)=y;%將當前位移記錄到s數組中以便作圖end,t=t0:dt:dt*(N-1);%仿真時間長度plot(t,s);xlabel(時間t);ylabel(位移y(t);axis(0501.1);,%ch1example2prg2.mg=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球質量t0=0;%起始時間K=0.85;%彈跳的損耗系數N=5000;%仿真的總步進數dt=0.005;%仿真步長v=v0;%初狀態(tài)y=y0;,5)仿真程序的功能擴展-以動態(tài)方式來觀察物體墜落的過程,fork=1:Nify0%小

24、球在空中的動力方程計算v=v-g*dt;y=y+v*dt;else%碰擊瞬間的計算y=-K.*v*dt;v=-K.*v-g*dt;endplot(0,y,o);axis(-2201);set(gcf,doubleBuffer,on)%雙緩沖避免作圖閃爍drawnow;end,實例3試仿真得出一個幅度調制系統(tǒng)的輸入輸出波形。設輸入被調制信號是一個幅度為2v，頻率為1000Hz的余弦波，調制度為0.5,調制載波信號是一個幅度為5v，頻率為10KHz的余弦波。所有余弦波的初相位為0。,（1）數學模型,該調幅系統(tǒng)的輸入輸出關系表達式為：,y(t)=(M+maMcos2fmt)cos2fct,其中，M=

25、2是被調信號的振幅，fm=1000是其頻率。A=5是載波信號的幅度，fc=104是其頻率。ma=0.5是調制度。,（2）編程實現(xiàn)在這個系統(tǒng)中的信號最高工作頻率為(fm+fc)=11KHz，根據取樣定理，只要離散取樣率高于該頻率的2倍即可無失真。在計算量和數據存儲量許可的條件下，取樣率可以設置更高，以使仿真計算的結果波形圖顯示更加光滑。我們設置取樣率為105，即在一個載波周期上取樣10次，相應的取樣間隔為t=10-5s。本例中，取樣間隔也作為仿真步進。程序代碼如下:,%ch2example1prg1.mdt=1e-5;%仿真采樣間隔T=3*1e-3;%仿真終止時間t=0:dt:T;input=2

26、*cos(2*pi*1000*t);%輸入被調信號carrier=5*cos(2*pi*1e4*t);%載波output=(2+0.5*input).*carrier;%調制輸出,%作圖:觀察輸入信號,載波,以及調制輸出subplot(3,1,1);plot(t,input);xlabel(時間t);ylabel(被調信號);subplot(3,1,2);plot(t,carrier);xlabel(時間t);ylabel(載波);subplot(3,1,3);plot(t,output);xlabel(時間t);ylabel(調幅輸出);,在實際中用雙蹤示波器來測量調幅的基帶輸入波形和調制輸

27、出波形時，只要基帶輸入余弦波的頻率與調幅載波頻率不是整數倍關系，示波器就只能夠對其中一個信號進行同步顯示，而另外一個信號在示波器上是滑動的。,如果調整示波器的掃描周期使之與兩個測試信號的周期都不呈整數倍關系，那么示波器屏幕上顯示的兩個信號都是滑動。能否通過仿真程序來再現(xiàn)這種波形顯示“滑動”現(xiàn)象呢？,（3）讓輸出波形“動”起來,考慮一個更接近真實物理系統(tǒng)的調幅模型?？紤]了相位噪聲n(t)的載波表達為：,v(t)=Acos(2fct+n(t)）,設信道是無衰減的，其加性噪聲為n(t)，那么在接收機所收到的調幅信號r(t)為：,r(t)=y(t)+n(t),程序中，我們將仿真時間區(qū)間劃分為若干段，每

28、段稱為一幀，示波器的掃描周期等于幀周期，即我們每仿真得出一幀的數據就讓顯示刷新一次。程序中故意將輸入被調信號的頻率設置為1005Hz，這樣其信號周期就與仿真的幀周期不是整數倍關系，運行后將看到“不斷”滑動的被調信號，載波顯示出相位抖動現(xiàn)象，而接收信號則沾染了噪聲。,%ch2example1prg3.mdt=1e-6;%仿真采樣間隔T=2*1e-3;%仿真的幀周期forN=0:500%總共仿真的幀數t=N*T+(0:dt:T);%幀中的取樣時刻input=2*cos(2*pi*1005*t);%輸入被調信號carrier=5*cos(2*pi*(1e4)*t+0.1*randn);%載波outp

29、ut=(2+0.5*input).*carrier;%調制輸出noise=randn(size(t);%噪聲r=output+noise;%調制信號通過加性噪聲信道,%作圖:觀察輸入信號,載波,以及調制輸出subplot(3,1,1);plot(0:dt:T,input);xlabel(時間t);ylabel(被調信號);text(T*2/3,1.5,當前幀數：N=,num2str(N);subplot(3,1,2);plot(0:dt:T,carrier);xlabel(時間t);ylabel(載波);subplot(3,1,3);plot(0:dt:T,r);xlabel(時間t);yla

30、bel(調幅輸出);set(gcf,DoubleBuffer,on);%雙緩沖避免作圖閃爍drawnow;end,實例4仿真電容的充電過程。,一個電壓源通過電阻與電容串聯(lián)的網絡對電容充電。設t=0為初始時刻（初始時刻之前電路斷開，不工作），電壓源輸出電壓x(t)為單位階躍函數，電容兩端的電壓為y(t)，回路電流為i(t)，并將電壓源視為系統(tǒng)輸入，電容上的電壓視為系統(tǒng)輸出。電路的初始狀態(tài)為y(0)。如圖：,電容的充電電路以及等價系統(tǒng),（1）數學分析首先根據網絡拓撲和元件伏安特性建立該電路方程組,并化簡得：,該方程也稱為系統(tǒng)的狀態(tài)方程。在方程中，變量y代表電容兩端的電壓，是電容儲能的函數。本例中它既是系統(tǒng)的狀態(tài)變量，又是系統(tǒng)的輸出變量。,（2）數值求解最直接的求解方法是將上式轉換為以時間向前遞進的計算遞推式，并以微小仿真時