【转载】让你的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。不幸的是这是一个大家常用的命令，有些书上还作为一条规则确立起来，建议必须使用!要知道这个命令一执行，工作空间的数据可就不可逆转的全没了啊!

警告：注意呦! 

小心同名覆盖。如果你在一个文件中，本来你的意思是两个变量，你却给他们起了相同的名字，那么第一次的数据可就没了。比如：

1. result=max(a,b); %想求 a和 b之较大者  
2. result=max(c,d); %想求 c 和d之较大者 

result结果是什么?恐怕不是你想要的。不妨将其改为result1和 result2。类似的，也要小心文件重名的覆盖，这个后果貌似更严重些。

下条内容请重视！

如何让 MATLAB崩溃。

尽管 MATLAB是很稳定的，但是我们仍然可以让它崩溃!使用第三方的MEX函数或者耗内存的操作比如视频处理或者超大规模矩阵都可能会造成MATLAB崩溃。

如果你已经有这些好习惯，那么恭喜，你要是还有其他好习惯麻烦也告诉我一声!如果没有，相信你看完之后总该有了吧?好了，我们开始!

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中看结果

下面我们评估一下下面这个函数：

1. function result =example1(count)   
2. for k = 1:count   
3.         result(k) = sin(k/50);   
4.        
5.         if result(k)<-0.9   
6.            result(k) = gammaln(k);   
7.         end   
8. end  

为了分析这个函数的效率，首先开启并清空 profiler，然后运行这个函数，接下来看结果报告。即依次输入:

1. >> profile on, profile clear   
2. >> example1(5000);   
3. >> 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。以上面程序为例,在命令行输入：

1. >> tic;example1(5000);toc 

输出是：

1. Elapsed time is 0.058522 seconds. 

为了获得更为精准的结果，你最好把浏览器、杀毒软件、防火墙等等占用CPU时间片的程序先关了，只剩下不能关掉的系统进程。

注意：profile在新版本中不断被加强，可使用的参数也越来越多，不过大多数根本用不着，如果你觉得那些参数很有用，我相信你根本用不找看我这个小册子了，要真是这样，麻烦您不吝赐教，分享一些经验。更详细的内容，您还是去看文档去吧!

2. 预分配矩阵

MATLAB中的矩阵变量可以动态增长行和列。比如：

1. >>x=2   
2. x=   
3. 2   
4. >>x(2,3)=1   
5. x=   
6. 2   0   0   
7. 0   0   1  

看到没?MATLAB自动调整了矩阵的大小!从内部实现上看，矩阵数据存储单元被重新分配了更大的单元。如果矩阵的大小被反复的调整(比如在循环中)，重新分配存储空间带来的额外开销会是很显著的。为了避免反复的矩阵存储重新分配，预分配矩阵的存储单元是一个不错的选择。一个推荐的方法是使用 zeros 函数命令。看下面的代码：

1. a(1) = 1;   
2. b(1) = 0;   
3. for k = 2:8000   
4. a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);   
5. b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);   
6. end   
7. tic/toc计时运行得到：   
8. Elapsed time is 0.013306 seconds.  

简单分析上面的代码，知道，每一次 for，矩阵 a 和 b 的大小都要被重新分配，最终的大小事 8000 的列向量。如果我们提前就给它们分配好大小为 8000的存储空间，看看结果怎么样：

1. a=zeros(1,8000);    %预分配矩阵存储单元   
2. b=zeros(1,8000);   
3. a(1) = 1;   
4. b(1) = 0;   
5. for k = 2:8000   
6. a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);   
7. b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);   
8. end   
9. 及时运行得到：   
10. Elapsed time is 0.000753 seconds.  

看出来没?速度提高了近 18 倍!像这种只需添加几行代码就能做到的情况是很多的。这个例子也有特殊性，就是最后的结果大小已知，如果结果的大小可变、未知呢?没关系，我们可以估计一下，最终结果最大能是多少?比估计到的最大再留出一些余量就成了!如果你估计的还是不够大，那超出的部分还要反复重新分配，不过这样节省下来的时间也是很可观的，毕竟可以少分配很多次了! 最后呢，还要处理一下后事，比如你分配给变量 a 有 1000 个单元，但最终它只占了300个，那你还要将那700个给收回来。看下面的代码：

1. a = zeros(1,10000);     %  预分配   
2. count = 0;   
3. for k = 1:10000   
4. v = exp(rand*rand);   
5. if v > 0.5        %  增长结果不确定的来源   
6. count = count + 1; a(count) = v;   
7. end   
8. end   
9. a = a(1:count);    %调整矩阵大小   
10. 未预分配时：Elapsed time is 0.052395 seconds.   
11. 预分配后：Elapsed time is 0.008935 seconds. 

感慨：些微时间的意义在哪呢?背后是你对 MATLAB 的理解深度。哥玩的不是时间，是技术。

3. 向量化

很多情况下，程序中的某些代码可以被向量化，向量化前后的速度往往在10 倍以上!向量化是最基本和最有效的让代码快起来的技巧，我都不愿意在后面叫“之一”了。

（1）向量化的计算

很多常规函数都是向量化的，它们作用于数组时，就好像是作用于数组中的每一个元素。例如：

1. >> sqrt([1,4,9,16])   
2. ans =   
3.       1     2     3     4   
4. 考虑下面的函数：   
5. function d = minDistance(x,y,z)    %寻找点集中距离远点最近点   
6. nPoints = length(x);   
7. d = zeros(nPoints,1);    %  预分配   
8. for k = 1:nPoints    %  计算每一个点的距离   
9. d(k) = sqrt(x(k)^2 + y(k)^2 + z(k)^2);   
10. end   
11. d = min(d);    %  得到最小距离   
12. 取  x=[1 2 3 4 5 6]; y=[2 3 5 2 1 4];z=[9 2 3 2 1 5];   
13. 计时运行：Elapsed time is 0.008006 seconds. 

如果你写出上面类似的代码，说明你认真看了前面的内容。为d预分配空间确实为本例节省了不少时间。如果采用向量化计算，我们可以去掉for循环，直接计算向量。这里要隆重推出“.”运算符，它表示的是对应元素进行运算。有.*和./和.\和.\'和.^等。分别表示不带.运算的对应元素运算。假设A是方阵，A^2是矩阵的 2 次乘幂，而 A.^2 表示矩阵 A 中的元素各自求平方组成新的矩阵。考虑下面的代码：

1. function d = minDistance(x,y,z)   
2. d = sqrt(x.^2 + y.^2 + z.^2);    %  计算每一点的距离   
3. d = min(d);     
4. 计时运行：   
5. 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. >> [1 4 2]>[2 3 1]   
2. ans =   
3.       0     1     1 

两个数组“按元素”进行比较。向量的逻辑操作返回二进制的逻辑结果向量，即用0代表假，用1代表真。这为什么有用呢?因为MATLAB中有一些强劲的针对逻辑向量的函数。例如：

1. >> find([1,5,3] < [2,2,4])     
2. ans =   
3. 1      3   
4. >> any([1,5,3] < [2,2,4])     
5. ans =   
6. 1   
7. >> all([1,5,3] < [2,2,4])     
8. ans =   
9. 0 

其实，对一般向量(非逻辑向量)也是适用的!

1. >> find(eye(4)==1)   
2. ans =   
3.       1   
4.       6   
5.      11   
6.      16 

以上函数的用法请自己查阅函数说明。

4. 示例

（1）向量归一标准化

将一个向量v归一标准化，我们可是使用v = v/norm(v)，norm函数的作用是求模(范数)。

如果对一组向量 v(:,1), v(:,2),…进行归一标准化，可以使用一个循环计算v(:,k)/norm(v(:,k))。更好的策略是向量化计算：

1. vMag = sqrt(sum(v.ˆ2));   
2. v = v./vMag(ones(1,size(v,1)),:); 

（2）剔除元素

有时候，我们需要将矩阵中的符合某些条件的元素剔除，当然可以使用条件判断加循环。我们使用向量化剔除矩阵中的NaN和无穷两类数：

1. i = find(isnan(x) | isinf(x));   %在x中找到符合条件的数的位置   
2. x(i) = [];    %剔除它   
3. 或者，同样的功能：   
4. i = find(˜isnan(x) & ˜isinf(x));     %找到不符合的数   
5. x = x(i);     %保留它   
6. 进一步的，我们可以更加简化，省略中间变量：   
7. x(isnan(x) | isinf(x)) = [];     
8. 以及   
9. x = x(˜isnan(x) & ˜isinf(x)); 

（3）分段函数

信号分析中十分重要的 sinc(x)函数是分段的：x=0 时的值是 1，x!=0 时，sinc(x)=sin(x)/x。下面的代码使用向量化方法处理分段：

1. function y = sinc(x)   
2. y = ones(size(x));          %  先设所有的y都是1   
3. i = find(x ˜= 0);          %  找到非零x值   
4. y(i) = sin(x(i)) ./ x(i);      %  计算非零值处的函数值   
5. 更简洁的，可以写成：   
6. y = (sin(x) + (x == 0))./(x + (x == 0)) 

能看出来吗?里面用到了逻辑运算，实在是巧妙的很!

（4）其他

还有些不常用的，算了，知道也八辈子用不着，珍惜脑细胞吧!

感慨：向量、矢量、相量、复数、数组、矩阵，这些名词能分清楚么?能分清楚知道内涵也就是为什么要这样规定么?不会也别问我!

5. 内嵌简单函数

内嵌函数的意思就是将函数调用的函数的代码直接写到这个函数里面来。由于函数调用要做保护现场以及恢复现场等工作，也会额外增加一些时间消耗。如果调用的次数不是很多，这些时间是可以忽略的，但是当调用次数很多的时候(比如500次)，这个时间就很可观了!

什么样的被调用函数适合内嵌呢?正如标题所说，是简单的函数，特征呢就是这个函数只有几行代码。如果这个函数很复杂，代码很长，还是死了这个心吧，内嵌是内嵌了，可是你看不懂代码了，得不偿失。程序的可读性是非常重要的!

注意：必须是 M-File 实现的函数才能内嵌!

下面的代码演示一个反复调用median函数的内嵌方法。原代码：

1. y = zeros(size(x));      %  预分配   
2. for k = 3:length(x)-2   
3. y(k) = median(x(k-2:k+2));   
4. end   
5. 取  x=rand(1,2500);   
6. 计时运行：Elapsed time is 0.030949 seconds. 

下面我们试试内嵌。首先，要研究一下你要内嵌的函数，本例中就是median。在命令行中输入：edit median，发现它是使用sort进行工作的。将核心代码内嵌：

1. y = zeros(size(x));     
2. for k = 3:length(x)-2   
3. tmp = sort(x(k-2:k+2));   
4. y(k) = tmp(3); ;   
5. end   
6. 仍取x=rand(1,2500);   
7. 计时运行：Elapsed time is 0.011379 seconds. 

以上就是一个演示，可见时间确实省去了不少。为了确认你想内嵌的函数是否是用M-File实现的，你可以使用“edit 函数名”命令试试看。