第2章 MATLAB编程基础
2.1 M-文件
MATLAB中的M-文件可以是简单执行一系列MATLAB语句的源文件，也可以是接收自变量并产生一个或多个输出的函数。
M-文件由文本编辑器创建，并以filename.m形式的文件名存储，比如average.m以及filter.m。M-文件的组成部分如下：
2.1.1 函数定义行
函数定义行的形式为：
function [outputs]=name(inputs)
例如，某个计算两幅图像的求和与乘积（两个不同的输出）的函数应该具有如下形式：
function [s,p]=sumprod(f,g)
其中，f和g是输入图像，s是求和图像，p是乘积图像。名称sumprod可以任意定义，但function总是在左侧，注意，输出参量必须位于方括号内，而输入参量位于圆括号内。如果函数只有单个输出参量，可不用方括号而直接列出。如果函数没有输出，只需要使用function, 不需要使用括号或等号。函数名必须以字母开头，后面可以跟字母、数字、下划线，但不允许有空格
函数可以在命令提示符中调用，例如：
[s,p]=sumprod(f,g);
也可以被用作其他函数的元素，在这种情况下，这些函数就成为子函数。如果输出只有单个变量，也可以不写括号：
y=sum(x)
2.1.2 H1语句
H1语句是第一个文本行，也就是函数定义行后面的单独注释行。函数定义行和H1语句之间无空行或空格，H1语句的示例如下：
% SUMPROD Computes the sum and product of two images
当用户在MATLAB提示符处输入
>> help function_name
H1语句是最先出现的文本。输入lookfor keyword就会显示出所有含有字符串keyword的H1语句。
2.1.3 帮助文本
帮助文本是紧跟在H1语句后面的文本块，二者之间无空行。帮助文本用来为函数提供注释或在线帮助。当用户在提示符后键入help function_name时，MATLAB会显示函数定义行和第一个非注释行（执行语句或空白语句）之间的全部注释行。但帮助系统会忽略帮助文本块后面的任何注释行。
2.1.4 函数主体
函数主体包含了执行计算并给出输出变量赋值的所有MATLAB代码。
符号"%"后面的非H1语句或帮助文本的所有行都被认为是函数注释行，他们不是帮助文本的一部分。代码行的末尾也可附加注释。
M-文件可以在任何文本编辑器中进行创建和编辑，并以扩展名.m保存到指定目录下，通常保存在MATLAB搜索路径中。创建和编辑M-文件的另→种方法是在提示符处使用edit函数。例如，如果文件存在于MATLAB搜索路径的目录中或者在当前目录下，可键入:
>> edit sumprod
就会打开文件sumprod.m并进行编辑。如果找不到该文件，MATLAB会为用户提供用于创建该文件的选项.
2.2 算子
MATLAB有两种不同类型的算子。矩阵运算由线性代数的规则来定义，而数组运算可以逐个元素地执行。字符.用来区分数组运算与矩阵运算，如(A*B)表示传统意义的矩阵乘法，而(A.*B)则表示数组乘法，这种乘法是指两个大小相同的数组对应元素的乘积。换句话说，假如(C=A.*B)，就有(C(i,j)=A(i,j)*B(i,j))。当书写(B=A)这样的表达式时，MATLAB将做(B=A)的"记录"，但并不将(A)的数据复制到(B)中，除非在后面的程序中，(A)的内容有了变化。这一点很重要，因为使用不同的v变量来“存储”相同的内容有时候可以增强代码的清晰性和可读性。下表列出了MATLAB的算子：
| 运算符 |
名称 |
描述 |
| + |
数组和矩阵加法 |
a+b、A+B 或 a+B |
| - |
数组和矩阵减法 |
a-b、A-B、A-a 或 a-A |
| .* |
数组乘法 |
Cv=A.B、C(i,j)=A(i,j)B(i,j) |
| * |
矩阵乘法 |
A*B，即标准矩阵乘 |
| ./ |
数组右除 |
C=A./B、C(i,j)=A(i,j)/B(i,j) |
| . |
数组左除 |
C=A.B、C(i,j)=B(i,j)/A(i,j) |
| / |
矩阵右除 |
A/B是计算A*inv(B) |
|
矩阵左除 |
AB是计算inv(A)*B |
| .^ |
数组乘幂 |
如果(C=A.^B)，那么C(i,j)=A(i,j)^B(i,j) |
| ^ |
矩阵乘幂 |
查阅帮助 |
| .' |
向量和矩阵转置 |
A.'标准向量和矩阵转置 |
| ' |
向量和矩阵复共轭转置 |
A'，标准向量和矩阵复共轭转置 |
2.3 关系算子
| 算子 |
名称 |
| < |
小于 |
| <= |
小于等于 |
| --- |
-------- |
| > |
大于 |
| >= |
大于等于 |
| == |
等于 |
| ~= |
不等于 |
2.4 逻辑算子
| 算子 |
名称 |
| & |
与 |
| | |
或 |
| ~ |
非 |
| && |
标量"与" |
| || |
标量"或" |
算子&和|针对数组操作，他们分别针对输入元素执行"与"和"或"运算；算子&&和||仅用于标量操作，主要与各种形式的if、while和for循环结合使用
2.5 流程控制
| 语句 |
描述 |
| if |
与else和elseif结合使用，执行一组基于指定逻辑条件的语句 |
| for |
执行规定次数的一组语句 |
| while |
根据规定的逻辑条件，执行不确定次数的一组语句 |
| break |
终止执行for或while循环 |
| continue |
将控制传递给for或while循环的下一次迭代，跳出循环体中的剩余部分 |
| switch |
与case和otherwise结合使用，根据指定的值或字符串执行不同的语句组 |
| return |
返回调用函数 |
| try...catch |
如果在执行期间检测到错误，就改变流程控制 |
2.6 数组索引
(1 imes N)的数组被称为行向量，这种向量的元素可以使用单一索引值来访问。例如(v(1))是向量(v)的第一个元素，(v(2))是第二个元素：
>> v=[1 3 5 7 9]
v =
1 3 5 7 9
>> v(2)
ans =
3
使用转置算子(.')可将行向量转换为列向量：
>> w=v.'
w =
1
3
5
7
9
为访问元素块，使用冒号。例如，为访问向量(v)的前3个元素，可使用语句：
>> v(1:3)
ans =
1 3 5
类似的，我们可以访问向量(v)的第3个元素到最后一个元素：
>> v(3:end)
ans =
5 7 9
其中，end表示向量中的最后一个元素。还可以将向量用作索引以进入另一个向量，例如：
>> v([1 4 5])
ans =
1 7 9
此外，索引并不限于连续的元素，例如：
>> v(1:2:end)
ans =
1 5 9
其中，符号1:2:end表示索引从1开始计数，步长为2，当计数达到最后一个元素时停止：
>> x=[1 2 3 4 5 6 7 8]
x =
1 2 3 4 5 6 7 8
>> x(1:2:end)
ans =
1 3 5 7
在MATLAB中，矩阵可以很方便地用一列被方括号括起并用分号隔开地行向量来表示：
>> A=[1 2 3;4 5 6;7 8 9]
A =
1 2 3
4 5 6
7 8 9
从矩阵中选择元素和从向量中选择元素是一样地，但需要两个索引：一个确定行的位置，另一个对应于列。我们也可以选择整行和整列，或使用冒号作为索引来选择整个矩阵：
>> A(2,:)
ans =
4 5 6
>> sum(A(:))
ans =
45
>> A(:)
ans =
1
4
7
2
5
8
3
6
9
函数sum计算参量每一列的和，单冒号索引把A转换为列向量，将结果传给sum。
另一种非常有用的索引形式是逻辑索引，逻辑索引表达式的形式是A(D)，其中A是数组，D是与A大小相同的逻辑数组，表达式A(D)提取A中与D中1值对应的所有元素：
>> D=logical([1 0 0;0 0 1;0 0 0])
D =
3×3 logical 数组
1 0 0
0 0 1
0 0 0
>> A(D)
ans =
1
6
>> D=[1 0 0;0 0 1;0 0 0]
D =
1 0 0
0 0 1
0 0 0
>> A(D)
下标索引必须为正整数类型或逻辑类型。
对图像处理很有用的最后一种索引是线性索引。线性索引表达式使用单个下标来编制矩阵或高维数组的索引。对于(M imes N)的矩阵，元素(r,c)可以用单一的下标r+M(c-1)。这样一来，A(2,3)就可以用A(8)来选择：
>> A(8)
ans =
6
2.7 函数句柄、单元数组与结构
函数句柄是MATLAB数据类型，包含用于引用函数的信息。使用函数句柄的主要优点之一是可以在调用中把函数句柄作为参数传递给另一个函数。
函数句柄有两种不同类型，这两种类型都用函数句柄算子@来创建。第一种函数句柄类型是命名的函数句柄，为了创建命名的函数句柄，在算子@后边写上所需的函数：
>> f=@sin
f =
包含以下值的 function_handle:
@sin
可以通过调用函数句柄f来间接调用函数sin:
>> f(pi/4)
ans =
0.7071
>> sin(pi/4)
ans =
0.7071
第二种函数句柄类型是匿名的函数句柄，由MATLAB表达式而不是函数名构成。构建匿名函数的通用格式是：
@(inoput-argument-list) expression
例如下列匿名函数句柄计算输入值的平方值：
>> g=@(x) x.^2;
>> g(2)
ans =
4
下列函数句柄计算两个变量平方之和的平方根：
>> r=@(x,y) sqrt(x.^2+y.^2)
r =
包含以下值的 function_handle:
@(x,y)sqrt(x.^2+y.^2)
>> r(1,2)
ans =
2.2361
单元数组(cell array)提供了一种在变量名下组合混合的一组对象（例如数字、字符、矩阵以及其他单元数组）的方法。例如，假设使用：（1）大小为(512 imes512)像素的uint8图像f；(2)188*2数组行形式的二维坐标序列b；（3）包含两个字符名的单元数组char_array={'area','centroid'}（花括号用来包含单元数组的内容）。可以用单元数组将这三种不同的实体组织成单个变量C。
C={f,b,char_array}
在提示符处键入C，将输出下列结果：
>>c
C=
[512*512 uint8] [188*2 double] {1*2 cell}
换句话说，显示的输出不是各种变量的值，而是对它们的某些特性的描述。为了看到单元素的全部内容，可以把单元元素的数值位置附加在花括号中。例如，要查看char_array的内容，键入：
>>c{3}
ans=
'area' 'centroid'
在C的元素中用圆括号代替花括号，给出变量的描述：
>>c(3)
ans =
{1*2 cell}
最后需要指出，单元数组包括参数的副本，而不是参数的指针。在前述的例子中，如果C的任何参数在C创建后改变了，那么改变在C中不会反映过来。
结构与单元数组类似，他们都允许一组不同的数据集成到单个变量中。但与单元数组不同的是，单元数组中的单元由数字寻址，而结构元素由所谓的“字段”来寻址，例如，如果f是一幅输入图像：
function s=image_stats(f)
s.dm=size(f);
s.AI=mens2(f);
s.AIrows=mean(f,2);
s.AIcols=mean(f,1);
2.8 优化代码
MATLAB是专门为数组操作而设计的编程语言。有两种重要的优化MATLAB代码的方法：预分配数组和向量化循环。
预分配是在进入计算数组元素的for循环之前初始化数组。假设要创建一个MATLAB函数来计算：
[f(x)=sin(x/100pi),x=0,1,2,...,M-1
]
下面是这个函数的第一种形式：
function y=sinfun1(M)
x=0:M-1;
for k=1:numel(x)
y(k)=sin(x(k)/(100*pi));
end
(M=5)时的输出是：
>> sinfun1(5)
ans =
0 0.0032 0.0064 0.0095 0.0127
MATLAB函数tic和toc可用来测量函数的执行时间，先调用tic，然后调用这个函数，之后再调用toc:
>> tic;sinfun1(100);toc
时间已过 0.001218 秒。
正如刚才介绍那样，调用计时函数再测量时间中可产生较大的变化，在命令提示符测量时尤其明显，例如，重复前边的调用将得到不同的结果：
>> tic;sinfun1(100);toc
时间已过 0.000356 秒。
函数timeit可用于得到函数调用的可靠且可重复的时间测量，其调用语法是：
s=timeit(f)
其中，f是用于对函数定时的函数句柄，s是调用所需要的秒数。调用函数句柄f时不使用输入参量。例如，在(M=100)时对sinfun1使用timeit进行计时：
>> M=100;
>> f=@() sinfun1(M);
>> timeit(f)
ans =
2.6024e-05
这个timeit函数调用例子很好地证明了前面介绍地函数句柄地强大功能。因为能够接受没有输入的函数句柄，函数timeit与我们希望计时的函数参数无关。用timeit测量sinfun1的时间，取M=500,1000,1500,...,20000:
M=500:500:20000;
>> for k=1:numel(M)
f=@() sinfun1(M(k));
t(k)=timeit(f);
end
预分配意味着在循环开始之前把它初始化为希望的输出大小。通常，采用调用函数zeros来做预分配：
function y=sinfun2(M)
x=0:M-1;
y=zeros(1,numel(x));
for k=1:numel(x)
y(k)=sin(x(k)/(100*pi));
end
比较sinfun1(20000)和sinfun2(20000)所需的时间：
>> timeit(@() sinfun1(20000))
ans =
0.0021
>> timeit(@() sinfun2(20000))
ans =
5.2661e-04
MATLAB中的向量化是使用矩阵/向量算子的组合、索引技术和现有的MATLAB或工具箱函数来完全消除循环。作为示例，sinfun的第三种形式：
function y=sinfun3(M)
x=0:M-1;
y=sin(x./(100*pi));
函数sinfun3没有for循环。在MATLAB旧版本中，用矩阵和向量算子消除循环几乎总能得到有意义的加速。然而，MATLAB新版本可自动编译简单的for循环，例如sinfun2中的那个，可加快机器代码。许多for循环在MATLAB的旧版本中很慢，但在向量化版本中不再慢：
>> timeit(@() sinfun2(20000))
ans =
0.0014
>> timeit(@() sinfun3(20000))
ans =
4.6300e-04
下面，我们写两种MATLAB版本的函数，创建一幅以下面等式为基础的合成图像：$$f(x,y)=Asin(u_0x+v_0y)$$第一个函数twodsin1使用两个嵌套的for循环计算(f)：
function f=twodsin1(A,u0,v0,M,N)
f=zeros(M,N);
for c=1:N
v0y=v0*(c-1);
for r=1:M
u0x=u0*(r-1)
f(r,c)=A*sin(u0x+v0y);
end
end
在for循环之前，预分配步骤f=zeros(M,N)。使用timeit创建一幅大小为512*512像素的图像，看看这个函数用了多长时间：
>> timeit(@() twodsin1(1,1/(4*pi),1/(4*pi),512,512))
ans =
0.0152
可以用imshow的自动范围语法[]显示结果图像：
>> f=twodsin1(1,1/(4*pi),1/(4*pi),512,512);
>> imshow(f,[])
下面我们写一个该函数的向量化版本，格式语法如下：
[C,R]=meshgrid(c,r)
输入参量c和r分别是水平（行）和垂直（列）坐标（首先写出列）。函数meshgrid把坐标向量变换为两个数组C和R，它们可以来计算两个变量的函数求值结果。例如，下面的命令用meshgrid对整数范围为1到3的x和范围10到14的y计算函数(z=x+y)：
>> [X,Y]=meshgrid(1:3,10:14)
X =
1 2 3
1 2 3
1 2 3
1 2 3
1 2 3
Y =
10 10 10
11 11 11
12 12 12
13 13 13
14 14 14
>> Z=X+Y
Z =
11 12 13
12 13 14
13 14 15
14 15 16
15 16 17
因此，我们用meshgrid重写2D且没有循环的sine函数：
function f=twodsin2(A,u0,v0,M,N)
r=0:M-1;%Row
c=0:N-1;%Column
[C,R]=meshgrid(c,r);
f=A*sin(u0*R+v0*C);
用timeit测量执行速度
>> timeit(@() twodsin2(1,1/(4*pi),1/(4*pi),512,512))
ans =
0.0055
向量化版本的运行时间要少三分之二。
因为MATLAB每一个新版本对运行循环的速度都倾向于有所改进，所以在向量化MATLAB代码时，给出通用的指导是困难的。但是，向量化的代码常常比基于循环的代码更易读。