讓你的MATLAB代碼飛起來
MATLAB語言是一種被稱為是“演算紙”式的語言，因此追求的是方便性、靈活性以及交互性。在快速性上要比C語言這種性能強勁著稱的稍遜一籌。然而，通過一些手段，我們也能讓MATLAB語言快起來，甚至和C差不多了!
首先聲明：本文是一個初級教程，因此很多知識是假定你已經很熟悉了的;雖然我在討論讓代碼飛起來，但從來不會說最快有多快，究竟有多快你要自己感覺;作者水平不是很高，難免誤導你，小心甄別。
在正式討論之前，先看看這些好習慣你有沒有?
1. 使用 M-Lint
M-Lint是一個代碼分析檢查工具，它在你寫代碼的過程中實時交互，發現你代碼的問題，按照最佳性能和最可維護性給出修改建議。
注意：我可沒說是最正確!
如果沒有激活這個功能，依次使用File > Preferences > M-Lint，勾選Enable integrated M-Lint warning and error messages 。同時，還可以設定你的偏好。
激活后，在你寫代碼時就會實時交互了，錯誤的或者不推薦的部分會以紅色下劃線標出，鼠標經過紅色下劃線的語句或單詞，M-Lint給出提示信息。想一下子看遍全部提示信息。使用Tools >M-Lint > (Save and) Show M-Lint Report2。
注：首次“觀看”先提示先保存一下。
2. 組織
給每一個項目(project)建立一個單獨的文件夾。同屬于一個項目的文件保存在哪兒的都有，你找的時候就不費勁嗎!
寫頭部注釋，尤其是H1。第一行就是H1。MATLAB中的內置函數的 help的內容其實就是讀取的這個函數的頭部注釋。怎么寫，參照MATLAB內置函數。
將經常用到的控制臺命令存儲為腳本(script)。如果有些命令反復使用，還是存為腳本吧，沒別的意思，你要少敲多少次鍵盤啊!
3. 避免數據丟失
不要在腳本中使用 clear all。不幸的是這是一個大家常用的命令，有些書上還作為一條規則確立起來，建議必須使用!要知道這個命令一執行，工作空間的數據可就不可逆轉的全沒了啊!
警告：注意呦!
小心同名覆蓋。如果你在一個文件中，本來你的意思是兩個變量，你卻給他們起了相同的名字，那么第一次的數據可就沒了。比如：
- result=max(a,b); %想求 a和 b之較大者
- result=max(c,d); %想求 c 和d之較大者
result結果是什么?恐怕不是你想要的。不妨將其改為result1和 result2。類似的，也要小心文件重名的覆蓋，這個后果貌似更嚴重些。
下條內容請重視！
如何讓 MATLAB崩潰。
盡管 MATLAB是很穩定的，但是我們仍然可以讓它崩潰!使用第三方的MEX函數或者耗內存的操作比如視頻處理或者超大規模矩陣都可能會造成MATLAB崩潰。
#p#
如果你已經有這些好習慣，那么恭喜，你要是還有其他好習慣麻煩也告訴我一聲!如果沒有，相信你看完之后總該有了吧?好了，我們開始!
1.使用profile
profile，Longman 給出的解釋是：a short description that gives important details about a person, a group of people, or a place。
MATLAB中內置了一個叫做profile的工具，來協助評估程序,也就是對程序運行過程的一個short description吧。主要命令有：
profile on 開啟
profile off 關閉
profile clear 清空數據
profile viewer 在profiler中看結果
下面我們評估一下下面這個函數：
- function result =example1(count)
- for k = 1:count
- result(k) = sin(k/50);
- if result(k)<-0.9
- result(k) = gammaln(k);
- end
- end
為了分析這個函數的效率，首先開啟并清空 profiler，然后運行這個函數，接下來看結果報告。即依次輸入:
- >> profile on, profile clear
- >> example1(5000);
- >> profile viewer
這就是 profile 的基本語法。也有使用鼠標操作的方法，這里就不介紹了，那樣雖然直觀單遠不及使用，命令方便。
由于系統的不同，報告的結果一般是不一樣的。以下是我的系統得出的結果。
1.先看profile summary:
2.點擊example1鏈接，進入具體各小項的評估。
（1）調用函數(children)、被調用函數(parents)。本例中都沒有。如果被 profile 的對象有調用函數或者被調用函數的話，會給出相應的數據。
（2）時間在哪些行被消耗(Lines where the most time was spent):
從數據中我們可以看出哪些行消耗了多少時間(總時間和相對時間)，被調用了多少次，以及直觀的柱形圖。
（3）另一個有用的項目是 M-Lint 結果，給出了錯誤(警告、提示)所在的行，以及對應的建議修改信息，這些建議對代碼的改進是很有價值的信息：
（4）最下面還有一個函數列表，是（2）的另一種形式。看圖：
最右側是函數代碼，前有行號、每一行調用的次數和小號的時間。消耗時間最多的行被標示了出來。最紅的消耗時間最多。
profiler工具的時間分辨率不是很高，因此，如果你的代碼運行的時間很短，有時候恐怕不能感知到。這時候不妨人為的加入幾個循環，讓程序所運行幾次，然后進行分析。
必須指出，profile工具的作用主要是分析程序，獲得程序運行的信息。如果想要知道程序運行的精確時間，使用計時器 tic/toc。以上面程序為例,在命令行輸入：
- >> tic;example1(5000);toc
輸出是：
- Elapsed time is 0.058522 seconds.
為了獲得更為精準的結果，你最好把瀏覽器、殺毒軟件、防火墻等等占用CPU時間片的程序先關了，只剩下不能關掉的系統進程。
注意：profile在新版本中不斷被加強，可使用的參數也越來越多，不過大多數根本用不著，如果你覺得那些參數很有用，我相信你根本用不找看我這個小冊子了，要真是這樣，麻煩您不吝賜教，分享一些經驗。更詳細的內容，您還是去看文檔去吧!
#p#
2. 預分配矩陣
MATLAB中的矩陣變量可以動態增長行和列。比如：
- >>x=2
- x=
- 2
- >>x(2,3)=1
- x=
- 2 0 0
- 0 0 1
看到沒?MATLAB自動調整了矩陣的大小!從內部實現上看，矩陣數據存儲單元被重新分配了更大的單元。如果矩陣的大小被反復的調整(比如在循環中)，重新分配存儲空間帶來的額外開銷會是很顯著的。為了避免反復的矩陣存儲重新分配，預分配矩陣的存儲單元是一個不錯的選擇。一個推薦的方法是使用 zeros 函數命令。看下面的代碼：
- a(1) = 1;
- b(1) = 0;
- for k = 2:8000
- a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);
- b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);
- end
- tic/toc計時運行得到：
- Elapsed time is 0.013306 seconds.
簡單分析上面的代碼，知道，每一次 for，矩陣 a 和 b 的大小都要被重新分配，最終的大小事 8000 的列向量。如果我們提前就給它們分配好大小為 8000的存儲空間，看看結果怎么樣：
- a=zeros(1,8000); %預分配矩陣存儲單元
- b=zeros(1,8000);
- a(1) = 1;
- b(1) = 0;
- for k = 2:8000
- a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);
- b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);
- end
- 及時運行得到：
- Elapsed time is 0.000753 seconds.
看出來沒?速度提高了近 18 倍!像這種只需添加幾行代碼就能做到的情況是很多的。這個例子也有特殊性，就是最后的結果大小已知，如果結果的大小可變、未知呢?沒關系，我們可以估計一下，最終結果最大能是多少?比估計到的最大再留出一些余量就成了!如果你估計的還是不夠大，那超出的部分還要反復重新分配，不過這樣節省下來的時間也是很可觀的，畢竟可以少分配很多次了! 最后呢，還要處理一下后事，比如你分配給變量 a 有 1000 個單元，但最終它只占了300個，那你還要將那700個給收回來。看下面的代碼：
- a = zeros(1,10000); % 預分配
- count = 0;
- for k = 1:10000
- v = exp(rand*rand);
- if v > 0.5 % 增長結果不確定的來源
- count = count + 1; a(count) = v;
- end
- end
- a = a(1:count); %調整矩陣大小
- 未預分配時：Elapsed time is 0.052395 seconds.
- 預分配后：Elapsed time is 0.008935 seconds.
感慨：些微時間的意義在哪呢?背后是你對 MATLAB 的理解深度。哥玩的不是時間，是技術。
#p#
3. 向量化
很多情況下，程序中的某些代碼可以被向量化，向量化前后的速度往往在10 倍以上!向量化是最基本和最有效的讓代碼快起來的技巧，我都不愿意在后面叫“之一”了。
（1）向量化的計算
很多常規函數都是向量化的，它們作用于數組時，就好像是作用于數組中的每一個元素。例如：
- >> sqrt([1,4,9,16])
- ans =
- 1 2 3 4
- 考慮下面的函數：
- function d = minDistance(x,y,z) %尋找點集中距離遠點最近點
- nPoints = length(x);
- d = zeros(nPoints,1); % 預分配
- for k = 1:nPoints % 計算每一個點的距離
- d(k) = sqrt(x(k)^2 + y(k)^2 + z(k)^2);
- end
- d = min(d); % 得到最小距離
- 取 x=[1 2 3 4 5 6]; y=[2 3 5 2 1 4];z=[9 2 3 2 1 5];
- 計時運行：Elapsed time is 0.008006 seconds.
如果你寫出上面類似的代碼，說明你認真看了前面的內容。為d預分配空間確實為本例節省了不少時間。如果采用向量化計算，我們可以去掉for循環，直接計算向量。這里要隆重推出“.”運算符，它表示的是對應元素進行運算。有.*和./和.\和.'和.^等。分別表示不帶.運算的對應元素運算。假設A是方陣，A^2是矩陣的 2 次乘冪，而 A.^2 表示矩陣 A 中的元素各自求平方組成新的矩陣。考慮下面的代碼：
- function d = minDistance(x,y,z)
- d = sqrt(x.^2 + y.^2 + z.^2); % 計算每一點的距離
- d = min(d);
- 計時運行：
- Elapsed time is 0.005326 seconds.
貌似差別不大?這就對了，別忘了，咱可就計算了6個值啊!這么幾個值就有了這樣的差距，那x、y、z向量要是大一點，結果的差異就可想而知了!
更進一步的，我們可以使用d = sqrt(min(x.^2 + y.^2 + z.^2))取代后兩行語句，讓程序更加簡潔。
一下函數使用向量化的計算會更為節省時間：min, max, repmat, meshgrid,sum, diff, prod等等。
（2）向量化邏輯
上面討論了計算的向量化，其實MATLAB的邏輯運算也是向量化的。比如：
- >> [1 4 2]>[2 3 1]
- ans =
- 0 1 1
兩個數組“按元素”進行比較。向量的邏輯操作返回二進制的邏輯結果向量，即用0代表假，用1代表真。這為什么有用呢?因為MATLAB中有一些強勁的針對邏輯向量的函數。例如：
- >> find([1,5,3] < [2,2,4])
- ans =
- 1 3
- >> any([1,5,3] < [2,2,4])
- ans =
- 1
- >> all([1,5,3] < [2,2,4])
- ans =
- 0
其實，對一般向量(非邏輯向量)也是適用的!
- >> find(eye(4)==1)
- ans =
- 1
- 6
- 11
- 16
以上函數的用法請自己查閱函數說明。
#p#
4. 示例
（1）向量歸一標準化
將一個向量v歸一標準化，我們可是使用v = v/norm(v)，norm函數的作用是求模(范數)。
如果對一組向量 v(:,1), v(:,2),…進行歸一標準化，可以使用一個循環計算v(:,k)/norm(v(:,k))。更好的策略是向量化計算：
- vMag = sqrt(sum(v.ˆ2));
- v = v./vMag(ones(1,size(v,1)),:);
（2）剔除元素
有時候，我們需要將矩陣中的符合某些條件的元素剔除，當然可以使用條件判斷加循環。我們使用向量化剔除矩陣中的NaN和無窮兩類數：
- i = find(isnan(x) | isinf(x)); %在x中找到符合條件的數的位置
- x(i) = []; %剔除它
- 或者，同樣的功能：
- i = find(˜isnan(x) & ˜isinf(x)); %找到不符合的數
- x = x(i); %保留它
- 進一步的，我們可以更加簡化，省略中間變量：
- x(isnan(x) | isinf(x)) = [];
- 以及
- x = x(˜isnan(x) & ˜isinf(x));
（3）分段函數
信號分析中十分重要的 sinc(x)函數是分段的：x=0 時的值是 1，x!=0 時，sinc(x)=sin(x)/x。下面的代碼使用向量化方法處理分段：
- function y = sinc(x)
- y = ones(size(x)); % 先設所有的y都是1
- i = find(x ˜= 0); % 找到非零x值
- y(i) = sin(x(i)) ./ x(i); % 計算非零值處的函數值
- 更簡潔的，可以寫成：
- y = (sin(x) + (x == 0))./(x + (x == 0))
能看出來嗎?里面用到了邏輯運算，實在是巧妙的很!
（4）其他
還有些不常用的，算了，知道也八輩子用不著，珍惜腦細胞吧!
感慨：向量、矢量、相量、復數、數組、矩陣，這些名詞能分清楚么?能分清楚知道內涵也就是為什么要這樣規定么?不會也別問我!
#p#
5. 內嵌簡單函數
內嵌函數的意思就是將函數調用的函數的代碼直接寫到這個函數里面來。由于函數調用要做保護現場以及恢復現場等工作，也會額外增加一些時間消耗。如果調用的次數不是很多，這些時間是可以忽略的，但是當調用次數很多的時候(比如500次)，這個時間就很可觀了!
什么樣的被調用函數適合內嵌呢?正如標題所說，是簡單的函數，特征呢就是這個函數只有幾行代碼。如果這個函數很復雜，代碼很長，還是死了這個心吧，內嵌是內嵌了，可是你看不懂代碼了，得不償失。程序的可讀性是非常重要的!
注意：必須是 M-File 實現的函數才能內嵌!
下面的代碼演示一個反復調用median函數的內嵌方法。原代碼：
- y = zeros(size(x)); % 預分配
- for k = 3:length(x)-2
- y(k) = median(x(k-2:k+2));
- end
- 取 x=rand(1,2500);
- 計時運行：Elapsed time is 0.030949 seconds.
下面我們試試內嵌。首先，要研究一下你要內嵌的函數，本例中就是median。在命令行中輸入：edit median，發現它是使用sort進行工作的。將核心代碼內嵌：
- y = zeros(size(x));
- for k = 3:length(x)-2
- tmp = sort(x(k-2:k+2));
- y(k) = tmp(3); ;
- end
- 仍取x=rand(1,2500);
- 計時運行：Elapsed time is 0.011379 seconds.
以上就是一個演示，可見時間確實省去了不少。為了確認你想內嵌的函數是否是用M-File實現的，你可以使用“edit 函數名”命令試試看。
【編輯推薦】
