MATLAB 7.0 采用行向量来表示多项式，将多项式的系数按降幂次序存放在行向量中。
例子：设计函数 poly2str()将行向量表示的多项式转换为常见的字符串表 多项式，
%poly2str.m 将多项式行向量表示转换为字符串表示
function Y=poly2str(X) %X 表示多项向量 %Y 多项字符串表示 如果%输入检查 X 退出不是一个向量
if isvector(X)==0,
disp输入错误：输入 X 不是一个向量，请输入一个代表多项式的向量!');
return; %函数返回
end; Y=''; %输出字符串
n=length(X);
for i=1:n, 将多项式的每一个功率转换为字符串
if(i~=1&&X(i)>0) 如果%是正系数，必须添加‘ ’字符
Y=[Y ' '];
end; %输出系数
if(X(i)==0), 如果这次的幂系数是% 0，不输出字符串
continue;
elseif(X(i)==1&&i~=n), 如果这次的幂系数是% 不输出系数，只输出 x^n
Y=Y;
else
Y=[Y num2str(X(i))]; 输出系数在其他情况下%
end; %输出 x^n
if(i==n-1), %1次力输出字符串x’
Y=[Y 'x'];
elseif(i==n), 百次不输出 x^n
Y=Y;
else
Y=[Y 'x^' num2str(n-i)]; 输出%其他情况 x^n
end; 如果不是最后一项，%输出’ ’
end; if(Y(1)==' ') %修正如果 0 次幂为 0 当时，字符串末尾有多余的字符串。 ’
Y(1)=[];
end;
实例：设计一个函数，实现从多项字符串表示到多项行向量表示的转换，需要输入 不码设置如下：
%str2poly.m %把多项式的字符串表示转换为行向量表示
function Y=str2poly(X) %X 多项式字符串形式 %Y 多项式是行向量形式 %输入格式检查
if(ischar(X)==0)
disp输入错误：输入 X 必须是字符串！
end; %用正则表达式寻找 或-的下标位置
index=regexp(X,'\ |\-'); 多项式项数%
L=length(index); 单元字符串矩阵用于存储多项式每个信息
term=cell(1,L 1);
term(1)=cellstr(X(1:(index(1)-1))；
for i=1:L-1,
term(i 1)=cellstr(X(index(i):(index(i 1)-1)));
end;
term(L 1)=cellstr(X(index(L):end)); 如果第一项为空，则删除该项目
if(isempty(char(term(1))))
term(1)=[]; 多项式项数减一%
L=L-1;
end; %多项系数矩阵
coefficient=[]; %多项式幂次矩阵与多项式系数矩阵一一对应
power=[];
for i=1:L 1, 单项多项字符串表达式
substring=char(term(i)); %用正则表达式寻找字符串x^’，由于’^正则表达式中的特殊字符大多使用\^’来表示
index2=regexp(substring,'x\^');
if(isempty(index2)==0), %如果匹配
if(index2(1)==1), %单项多项式字符串为’x^*’形式
coefficient=[coefficient 1];
power=[power str2num(substring((index2(1)+2):end))];
elseif(index2(1)==2), %单项多项式字符串为’+x^*’或者’-x^*’，’3x^*’形式
if(substring(1)=='+'),
coefficient=[coefficient 1];
power=[power str2num(substring((index2(1)+2):end))];
elseif(substring(1)=='-'),
coefficient=[coefficient -1];
power=[power str2num(substring((index2(1)+2):end))];
end;
else %单项多项式字符串为’2.2x^*’
coefficient=[coefficient str2num(substring(1:(index2(1)-1)))];
power=[power str2num(substring((index2+2):end))];
end;
else %单项多项式字符串不含’x^’ %用正则表达式，寻找字符串’x’
index2=regexp(substring,'x');
if(isempty(index2)==0), %如果匹配上’x’
if(index2(1)==1), %单项多项式字符串为’x*’形式
coefficient=[coefficient 1];
power=[power 1];
elseif(index2(1)==2), %单项多项式字符串为’+x*’或者’-x*’，’3x*’形式
if((substring(1)=='+')==1),
coefficient=[coefficient 1];
power=[power 1];
elseif(substring(1)=='-'),
coefficient=[coefficient -1];
power=[power 1];
else %单项多项式字符串为’2.2x*’
coefficient=[coefficient str2num(substring(1:(index2-1)))];
power=[power 1];
end;
else
coefficient=[coefficient str2num(substring(1:(index2-1)))];
power=[power 1];
end;
else %单项多项式字符串不含’x^’,’x’,则断定它是常数项
coefficient=[coefficient str2num(substring)];
power=[power 0];
end;
end;
end; %合并同类项
N=max(power)+1;
Y=zeros(1,N);
for i=1:N,
index3=find(power==(N-i));
Y(i)=sum(coefficient(index3));
end;
调用代码
X=’-10x^2+8x^2-x^3+4x^5+x^3-10+9-x^7’;
P=str2poly(X)
Str=poly2str(P)
多项式求值
在 MATLAB 7.0 中使用函数 polyval()来计算多项式的值
实例：求多项式的值和矩阵多项式的值，请注意这两者的区别，具体代码如下：
%polyval_example.m %多项式求值
str1='x^2+2x+3';
p1=str2poly(str1);
A=[1 0;
0 -1]; p1_A=polyval(p1,A) %求多项式值
p1_Am=polyvalm(p1,A) %求矩阵多项式的值
多项式乘法和多项式除法
多项式乘法和除法运算也就是多项式向量的卷积和解卷积运算。在 MATLAB 7.0 中，函 数 conv()和 deconv()用来实现这些运算。这两个函数的调用格式如下。 • w = conv(u,v)，实现向量 u,v 的卷积，在代数上相当于多项式 u 乘上多项式； • [q,r] = deconv(v,u)，实现解卷积运算，他们之间的关系为 v = conv(u,q)+r。在代数上相 当于实现多项式 u 除以 v，得到商 q 和余多项式 r。
实例：求多项式 x2+2x+3 与多项式 2x2+3x+4 的乘积
%poly_conv.m %多项式乘法
str1=’x^2+2x+3’;
str2=’2x^2+3x+4’;
p1=str2poly(str1);
p2=str2poly(str2);
p3=conv(p1,p2);
str3=poly2str(p3)
例如，求上个例子中的多项式乘积结果 2x4+7x3+16x2+17x+12 除以多项式 x2+2x+3 的商，
%poly_deconv.m %多项式除法
str3=’2x^4+7x^3+16x^2+17x+12’;
str1=’x^2+2x+3’;
p3=str2poly(str3);
p1=str2poly(str1);
[p2 r]=deconv(p3,p1);
str2=poly2str(p2)
r
多项式的导数和微分
1．多项式的导数
在 MATLAB 7.0 中使用函数 polyder()来计算多项式的导数，其调用方式如下： • k = polyder(p)，返回多项式 p 的导数； • k = polyder(a,b)，返回多项式 a 与多项式 b 乘积的导数； • [q,d] = polyder(b,a)，返回多项式 a 除以 b 的商的导数，并以 q/d 格式表示
实例：，对多项式求导，代码设置如下：
%polyder_example.m %多项式求导数的示例
str1=’x^2+2x+3’;
str2=’2x^2+3x+4’;
p1=str2poly(str1);
p2=str2poly(str2); %求多项式’x^2+2x+3’的导数
q1=polyder(p1);
str1_der=poly2str(q1) %求多项式’2x^2+3x+4’的导数
q2=polyder(p2);
str2_der=poly2str(q2) %求多项式’x^2+2x+3’与’2x^2+3x+4’乘积的导数
c=polyder(p1,p2);
str3_der=poly2str(c) %求多项式’x^2+2x+3’除以’2x^2+3x+4’的商的导数
[q d]=polyder(p1,p2);
str_der_q=poly2str(q)
str_der_q=poly2str(d)
．多项式的积分
在 MATLAB 7.0 中使用函数 polyint()来计算多项式的积分，其调用方式如下： • polyint(p,k)，返回多项式 p 的积分，设积分的常数项为 k； • polyint(p)，返回多项式 p 的积分，设积分的常数项为 0
实例：求多项式 3x2+4x+5 的积分，代码设置如下：
%polyint_example.m %多项式求积分的示例
str1=’3x^2+4x+5’;
p1=str2poly(str1); %求积分并设常数项为 5
p2=polyint(p1,5);
str2=poly2str(p2) %求积分并设常数项为 0
p3=polyint(p1);
str3=poly2str(p3)
多项式的根和由根创建多项式
1．多项式的根
在 MATLAB 7.0 中使用函数 roots()来求多项式的根，其调用方式如下： • r = roots(c)，返回多项式 c 的所有根 r，r 是向量，其长度等于根的个数。
实例：求多项式 2x3-2x2-8x+8 的根，代码设置如下
%roots_example.m %多项式求根
str1=’2x^3-2x^2-8x+8’;
p1=str2poly(str1); r=roots(p1) %求多项式根
2.由根创建多项式
与多项式求根相反的过程是由根创建多项式，它由函数 poly()来实现，其调用格式如下： • p = poly(r)，输入 r 是多项式所有根，返回值为代表多项式的行向量形式； • p = poly(A)，输入是 N×N 的方阵，返回值 p 是长度为 N+1 的行向量多项式，它是矩阵 A 的特征多项式，也就是说多项式 p 的根是矩阵 A 的特征值。
实例：，由根[-2 2 1]创建多项式
%poly_example.m %由根创建多项式
r=[-2 2 1];
p=poly(r);
str=poly2str(p)
实例：，求矩阵的特征多项式，代码设置如下
%poly_example2.m %求矩阵 A 的特征多项式
A =[1 2 3;
4 5 6;
7 8 0];
p=poly(A);
str=poly2str(p)；
r=roots(p) %求多项式 p 的根
A_eig=eig(A) %求矩阵 A 的特征值
多项式部分分式展开
函数 residue()可以将多项式之比用部分分式展开，也可以将一个部分分式表示为多项式 之比。
实例：，求多项式之比的部分分式展开，然后在把部分分式表示为多项式之比的形式，代 码设置如下：
%residue_example.m %求多项式之比的部分分式，然后把部分分式表示为多项式之比的形式
str1=’2x^3+3x^2-4x+1’;
str2=’x^2-3x+2’;
p1=str2poly(str1);
p2=str2poly(str2); [r,p,k]=residue(p1,p2) %求多项式之比的部分分式
[p1_res,p2_res]=residue(r,p,k); %把部分分式表示为多项式之比的形式 str1_res=poly2str(p1_res) %把多项式显示为字符串形式
str2_res=poly2str(p2_res)
多项式曲线拟合
函数 polyfit()采用最小二乘法对给定数据进行多项式拟合，最后给出多项式的系数。该 函数的调用方式如下： • p = polyfit(x,y,n)，采用 n 次多项式 p 来拟合数据 x 和 y，从而使得 p(x)与 y 最小均方差 最小
实例：下面例子说明函数 polyfit()的用法，并讨论采用不同多项式阶数对拟合结果的影 响
%polyfit_example.m %说明函数 polyfit()的用法，并讨论采用不同多项式阶数对拟合结果的影响
x=0:0.2:10;
y=0.25*x+20*sin(x); %5 阶多项式拟合
p5=polyfit(x,y,5);
y5=polyval(p5,x); %8 阶多项式拟合
p8=polyfit(x,y,8);
y8=polyval(p8,x); %60 阶多项式拟合
p60=polyfit(x,y,60);
y60=polyval(p60,x); %画图
hold on;
plot(x,y,’ro’);
plot(x,y5,’b--’);
plot(x,y8,’b:’);
plot(x,y60,’r-.’);
xlabel(’x’);
ylabel(’y’); legend(’原始数据’,’5 阶多项式拟合’,’8 阶多项式拟合’,’60 阶多项式拟合’);
%得到输出如图 6-1 所示。由图 6-1 可以看出：使用 5 次多项式拟合时，拟合得到的结果 比较差。而使用 8 次多项式拟合时，得到的结果与原始数据符合得很好。但使用 60 次多项式 拟合时，拟合的结果非常差。可见用多项式拟合必须选择适中的阶数，而不是阶数越高精度 越高。
曲线拟合图形用户接口
为了方便用户的使用，MATLAB 7.0 提供了支持曲线拟合的图形用户接口。它位于 “Figure”窗口的“Tools\Basic Fitting”菜单中。为了使用该工具，首先用待拟合的数据画图， 代码设置如下：
x=0:0.2:10;
y=0.25*x+20*sin(x);
plot(x,y,’ro’);
插值
插值是在已知数据之间寻找估计值的过程。在信号处理和图像处理中，插值是极其常用 的方法。MATLAB 7.0 提供大量的插值函数。这些函数在获得数据的平滑度、时间复杂度和 空间复杂度方面上有不同的性能。 插值函数保存在 MATLAB 7.0 工具箱的 polyfun 子目录下
函数名 功能描述 pchip 分段三次厄米多项式插值 interp1 一维插值 interp1q 一维快速插值 interpft 一维快速傅立叶插值方法 interp2 二维插值 interp3 三维插值 interpn N 维插值 griddata 栅格数据插值 griddata3 3 维栅格数据插值 griddatan N 维栅格数据插值 spline 三次样条插值 ppval 分段多项式求值
一维插值
一位插值就是对一维函数 y=f(x)进行插值，图 6-6 显示了插值的含义，实心点(x,y)代表的 是已知数据，空心点(xi,yi)的横坐标 xi 代表需要估计值的位置，纵坐标 yi 表示插值后的估计值。在 MATLAB 7.0 中，一维插值有基于多项式的插值和基于快速傅立叶的插值两种类型。
一维多项式插值可以通过函数 interp1()来实现
一维插值可以指定的方法如下： （1）最邻近插值（Nearest neighbor interpolation，method=’nearest’），这种插值方法在已知 数据的最邻近点设置插值点，对插值点的数进行四舍五入。对超出范围的点将返回一个 NaN （Not a Number）。 （2）线性插值（Linear interpolation，method=’linear’）,该方法是未指定插值方法时 MATLAB7.0 默认采用的方法。该方法采用直线连接相邻的两点。对超出范围的点将返回一 个 NaN（Not a Number）。 （3）三次样条插值（Cubic spline interpolation，method=’spline’），这样方法采用三次样条 函数来获得插值点。在已知点为端点情况下，插值函数至少具有相同的一阶和二阶导数。 （4）分段三次厄米多项式插值（Piecewise cubic Hermite interpolation，method=’pchip’）。 （5）三次多项式插值（method=’cubic’），与分段三次厄米多项式插值相同。 （6）MATLAB5 中使用的三次多项式插值（method=’v5cubic’），该方法使用一个 3 次多项 式函数对已知数据进行拟合。 当选择一个插值方法时应该考虑方法的执行速度、占用内存大小和获得数据的平滑度。 以上方法的特点如下。
（1）最邻近插值：最快的插值方法，但是数据平滑方面最差，其得到的数据是不连续的。 （2）线性插值：比邻近插值占用更多的内存，执行速度也稍慢，但其数据平滑方面优于 邻近插值。与邻近插值不同的是，线性插值的数据变化是连续的。 （3）三次样条插值：处理速度最慢，占用内存小于分段三次厄米多项式插值，可以产生 最光滑的结果。但是如果输入数据不均匀或者某些点靠得很近，会出现一些错误。样条插值 是极其有用的插值方法，因此除了提供三次样条插值函数外，MATLAB 7.0 还提供了一个样 条插值工具箱，它位于 toolbox\splines 下。 （4）分段三次厄米多项式插值：处理速度和消耗的内存比线性插值差，但插值得到的数 据和一阶导数都是连续的。 这些方法的相对优劣不仅适用于一维插值，而且适用于二维插值情况甚至高维插值 情况。
下面例子用不同插值方法对一维数据进行插值，并比较其不同，代码设置如下
%interp1_example.m %用不同插值方法对一维数据进行插值，并比较其不同
x = 0:1.2:10;
y = sin(x);
xi = 0:0.1:10; yi_nearest = interp1(x,y,xi,’nearset’); %最邻近插值 yi_linear = interp1(x,y,xi); %默认插值方法是线性插值 yi_spline = interp1(x,y,xi,’spline ’); %三次样条插值 yi_cubic = interp1(x,y,xi,’cubic’); %三次多项式插值 yi_v5cubic = interp1(x,y,xi,’v5cubic’); %MATLAB 5 中使用的三次多项式插值
hold on;
subplot(2,3,1);
plot(x,y,’ro’,xi,yi_nearest,’b-’); title(’最邻近插值’);
subplot(2,3,2);
plot(x,y,’ro’,xi,yi_linear,’b-’); title(’线性插值’);
subplot(2,3,3);
plot(x,y,’ro’,xi,yi_spline,’b-’); title(’三次样条插值’);
subplot(2,3,4);
plot(x,y,’ro’,xi,yi_cubic,’b-’); title(’三次多项式插值’);
subplot(2,3,5);
plot(x,y,’ro’,xi,yi_v5cubic,’b-’); title(’三次多项式插值(MATLAB5)’);
一维快速傅立叶插值
一维快速傅立叶插值通过函数 interpft()来实现，该函数用傅立叶变换把输入数据变换到 频域，然后用更多点的傅立叶逆变换，变换回时域，其结果是对数据进行增采样
实例：利用一维快速傅立叶插值实现数据增采样，代码设置如下：
%interpft_example.m %一维快速傅立叶插值实现数据增采样
x = 0:1.2:10;
y = sin(x); n = 2*length(x); %增采样 1 倍 yi = interpft(y,n); %一维快速傅立叶插值
xi = 0:0.6:10.4;
hold on;
plot(x,y,’ro’); %画图
plot(xi,yi,’b.-’);
title(’一维快速傅立叶插值’);
legend(’原始数据’,’插值结果’);
二维插值
二维插值主要应用于图像处理和数据的可视化，其基本思想与一维插值相同，它是对两 变量的函数 z=f(x,y)进行插值
MATLAB 7.0 中的二维插值函数为 interp2()，其调用格式如下。 • zi = interp2(x,y,z,xi,yi)，原始数据 x,y,z 决定插值函数 z=f(x,y)，返回值 zi 是(xi,yi)在函数 f(x,y)上的值； • zi = interp2(z,xi,yi)，若 z=n×m，则 x=1:n，y=1:m； • zi = interp2(z,ntimes)，在两点之间递归地插值 ntimes 次； • zi = interp2(x,y,z,xi,yi,method)，可采用的插值方法在本小节将详细介绍； • zi = interp2(...,method, extrapval)，当数据超过原始数据范围时，用输入 extrapval 指定 一种外推方法。 二维插值可以采用的插值方法如下： （1）最邻近插值（Nearest neighbor interpolation，method=’nearest’），这种插值方法在已知 数据的最邻近点设置插值点，对插值点的数进行四舍五入。对超出范围的点将返回一个 NaN （Not a Number）。 （2）双线性插值（Bilinear interpolation，method=’linear’）,该方法是未指定插值方法时 MATLAB7.0 默认采用的方法。插值点的值只决定于最邻近的 4 个点的值。 （3）三次样条插值（Cubic spline interpolation，method=’spline’），这样方法采用三次样条
函数来获得插值数据。
实例：）双三次多项式插值（method=’cubic’）。 例如，下面例子比较各种二维插值插值算法的不同，如下：
%interp2_example2 %采用二次插值对三维高斯型分布函数进行插值 [x,y] = meshgrid(-3:0.8:3); %原始数据
z = peaks(X,Y); [xi,yi] = meshgrid(-3:0.25:3); %插值点 zi_nearest = interp2(x,y,z,xi,yi,’nearset’); %最邻近插值 zi_linear = interp2(x,y,z,xi,yi); %默认插值方法是线性插值 zi_spline = interp2(x,y,z,xi,yi,’spline ’); %三次样条插值 zi_cubic = interp2(x,y,z,xi,yi,’cubic’); %三次多项式插值
hold on;
subplot(2,3,1);
surf(x,y,z); title(’原始数据’);
subplot(2,3,2);
surf(xi,yi,zi_nearest); title(’最邻近插值’);
subplot(2,3,3);
surf(xi,yi,zi_linear); title(’线性插值’);
subplot(2,3,4);
surf(xi,yi,zi_spline); title(’三次样条插值’);
subplot(2,3,5);
surf(xi,yi,zi_cubic); title(’三次多项式插值’); figure; %新开绘图窗口 subplot(2,2,1); %画插值结果的等高线
contour(xi,yi,zi_nearest); title(’最邻近插值’);
subplot(2,2,2);
contour(xi,yi,zi_linear); title(’线性插值’);
subplot(2,2,3);
contour(xi,yi,zi_spline); title(’三次样条插值’);
subplot(2,2,4);
contour(xi,yi,zi_cubic);
title(’三次多项式插值’);
数据分析和傅立叶变换
MATLAB7.0 中与数据分析和傅立叶变换相关的函数位于目录\toolboxs\MATLAB\datafun 下，本节将分类介绍这些函数。 MATLAB 7.0 在这些函数中的约定。 • 一维统计数据可以用行向量或者列向量来表示，不管输入数据是行向量或者是列向量， 运算是对整个向量进行的； • 二维统计数据可以采用多个向量表示，也可以采用二维矩阵来表示。对二维矩阵运算 时，运算总是按列进行的。 这两条约定不仅适用于本节提到的函数，而且适用于 MATLAB 7.0 各个工具箱的函数。
基本数据分析函数
MATLAB 7.0 提供的基本数据分析函数的功能和调用格式如表 6-3 所示。
表 6-3 基本数据分析函数 函数名 功能描述 基本调用格式
max 求最大值
C = max(A)，如果 A 是向量，返回向量中的最大值；如果 A 是矩阵 返回一个包含各列最大值的行向量 C = max(A,B)，返回矩阵 A 和 B 之中较大的元素，矩阵 A、B 必须 具有相同的大小 C = max(A,[],dim)，返回 dim 维上的最大值 [C,I] = max(...)，多返回最大值的下标 min 求最小值 与最大值函数 max()调用格式一致
mean 求平均值
M = mean(A)，如果 A 是向量，返回向量 A 的平均值，如果 A 是矩 阵，返回含有各列平均值的行向量 M = mean(A,dim)，返回 dim 维上的平均值 median 求中间值 与求平均值函数 mean()的调用格式一致
std 求标准方差
s = std(A)，如果 A 是向量，返回向量的标准方差；如果 A 是矩阵， 返回含有各列标准方差的行向量 s = std(A,flag)，用 flag 选择标准方差的定义式 s = std(A,flag,dim)，返回 dim 维上的标准方差
var 方差（标准方差的平反）
var(X)，如果 A 是向量，返回向量的方差；如果 A 是矩阵，返回含 有各列方差的行向量 var(X,1)，返回第二种定义的方差 var(X,w)，利用 w 做为权重计算反差 var(X,w,dim)，返回 dim 维上的方差
sort 数据排序
B = sort(A)，如果 A 是向量，升序排列向量；如果 A 是矩阵，升序 排列各个列 B = sort(A,dim)，升序排列矩阵 A 的 dim 维 B = sort(...,mode)，用 mode 选择排序方式：’ascend’为升序，’descend’ 为降序 [B,IX] = sort(...)，多返回数据 B 在原来矩阵中的下标 IX
sortrows 对矩阵的行排序
B = sortrows(A)，升序排序矩阵 A 的行 B = sortrows(A,column)，以 column 列数据作为标准，升序排序矩阵 A 的行 [B,index] = sortrows(A)，多返回数据 B 在原来矩阵 A 中的下标 IX
sum 求元素之和
B = sum(A)，如果 A 是向量，返回向量 A 的各元素之和；如果 A 是 矩阵，返回含有各列元素之和的行向量 B = sum(A, dim)，求 dim 维上的矩阵元素之和 B = sum(A, ’double’)，返回数据类型指定为双精度浮点数 B = sum(A, ’native’)，返回数据类型指定为与矩阵 A 的数据类型相同
prod 求元素的连乘积
B = prod(A)，如果 A 是向量，返回向量 A 的各元素连乘积；如果 A 是矩阵，返回含有各列元素连乘积的行向量 B = prod(A,dim)，返回 dim 维上的矩阵元素连乘积
hist 画直方图
n = hist(Y)，在 10 个等间距的区间统计矩阵 Y 属于该区间的元素个 数 n = hist(Y,x)，在 x 指定的区间统计矩阵 Y 属于该区间的元素个数 n = hist(Y,nbins)，用 nbins 个等间距的区间统计矩阵 Y 属于该区间的 元素个数 hist(...)，直接画出直方图
histc 直方图统计
n = histc(x,edges)，计算在 edges 区间内向量 x 属于该区间的元素个 数 n = histc(x,edges,dim)，在 dim 维上统计 x 出现的次数
trapz 梯形数值积分（等间距）
Z = trapz(Y)，返回 Y 的梯形数值积分 Z = trapz(X,Y)，计算以 X 为自变量时 Y 的梯形数值积分 Z = trapz(...,dim)，在 dim 维上计算梯形数值积分
cumsum 矩阵的累加
B = cumsum(A)，如果 A 是向量，计算向量 A 的累计和；如果 A 是 矩阵，计算矩阵 A 在列方向上的累计和 B = cumsum(A,dim)，在 dim 维上计算矩阵 A 的累计和 cumprod 矩阵的累积 函数调用格式与函数 cumsum()相同 cumtrapz 梯形积分累计 函数调用格式与函数 trapz()相同
实例：．最大值、最小值、平均值、中间值、元素求和
这几个函数的用法都比较简单，下面的例子用来说明这些函数的用法：
%max_example %最大值、最小值、平均值、中间值、元素求和函数使用的例子
x=1:40; y=randn(1,40); %正态分布随机数据
hold on;
plot(x,y); [y_max,I_max]=max(y); %求向量最大值和相应下标
plot(x(I_max),y_max,’o’); [y_min,I_min]=min(y); %求向量最小值和相应下标
plot(x(I_min),y_min,’*’); y_mean=mean(y); %求向量平均值
plot(x,y_mean*ones(1,length(x)),’:’); y_median=median(y); %求向量反差
plot(x,y_median*ones(1,length(x)),’-.’); y_sum=sum(y); %求向量元素之和
plot(x,y_sum*ones(1,length(x)),’--’);
xlabel(’x’);
ylabel(’y’); legend(’数据’,’最大值’,’最小值’,’平均值’,’中间值’,’向量元素之和’);
MATLAB 7.0 中默认使用（5-1）式来计算数据的标准方差。要使用（5-2）式来计算标 准方差可以使用函数调用格式 std(A,1)。 方差是标准方差的平反。对应标准方差的两种定义，方差也有两种定义。同样 MATLAB 7.0 中默认使用（5-1）式来计算数据的方差。要使用（5-2）式来计算方差可以使用函数调用 格式 var(A,1)。 例如，有一维随机变量 x，要求验证 2*x 的方差是 x 的 4 倍，2x 标准方差是 x 的 2 倍， 代码设置如下：
%std_example.m %验证 2*x 的方差是 x 的 4 倍，2x 标准方差是 x 的 2 倍
x=rand(1,1000);
x_std=std(x);
x_var=var(x);
x2_std=std(2*x);
x2_var=var(2*x); disp([’随机变量 x2 的标准方差与 x1 的标准方差之比 = ’ num2str(x2_std/x_std)]); disp([’随机变量 x2 的标准方差与 x1 的标准方差之比 = ’ num2str(x2_var/x_var)]);
．元素排序
MATLAB 7.0 可以对实数、复数和字符串进行排序。对复数矩阵进行排序时，先按复数 的模进行排序，如果模相等则按其在区间[-π,π]上的相角进行排序。MATLAB7.0 中实现排 序的函数为 sort()。
实例：对复数矩阵排列，代码设置如下：
%sort_complex.m %对复数进行排序
a = [0 -1 -0.5i 1 0.5 i -i -0.5 0.5i];
b = sort(a)
MATLAB 7.0 提供函数 sortrows()用于对矩阵的行进行排列。
实例：%minvar.m %x 是一个向量，k 是小于等于向量 x 长度的正标量 %返回值是 y 是 x 中反差最小的 k 个数据，i 是 y 的元素在 x 向量中的下标
function [y,i]=minvar(x,k)
if(isvector(x)==0), %输入数据检查
disp('输入数据 x 必须是向量!');
return;
end;
if(isscalar(k)==0) disp('输入数据 k 必须是标量!');
return;
end;
if(k<=0||k>length(x)),
disp('输入数据 k 必须是是小于等于向量 x 长度的正标量!');
return;
end; k=floor(k); %把 k 截断为整数
[x_sort index]=sort(x); %对数据 x 排序
n=length(x); %数据 x 的长度
for i=1:(n-k+1), %对排好序的数据求所有相邻 k 个值
std_x_sort(i)=std(x_sort(i:(i+k-1))); %的标准方差
end;
[min_std_x_sort j] = min(std_x_sort); %求最小的相邻 k 个值的标准方差和下标
i=index(j:(j+k-1)); %求最小标准方差的相邻 k 个值在向量
x %中的方差
y=x(i); %求最小标准方差的 k 个值
调用实例：
%minvar_example.m %验证已知 n 个数据中挑选出方差最小的 k 个数据的函数 minvar()
data = (1:10).^2; %数列 1 4 9 16 ...
index = randperm(10); %1:10的随机排列
x = data(index) %随机扰乱数据 data 的排列
[y i]=minvar(x,4) %求数据 x 中方差最小的 4 个值
协方差和相关系数矩阵
在数理统计中，协方差和相关函数是极其重要的随机变量的数字特征
%corrcoef_example.m %计算矩阵的相关系数
x = randn(10000,4); % 4个完全无关随机变量
y(:,1) = x(:,1);
y(:,2) = x(:,2);
y(:,3) = 0.3*x(:,3) + 0.7*x(:,4); %在第 3 个随机变量和第 4 个随机变量之间
y(:,4) = 0.2*x(:,3) + 0.8*x(:,4); %引入相关性
[r,p] = corrcoef(y) %计算相关矩阵和不相关概率矩阵
[i,j] = find(p<0.05); %寻找相关的随机变量对
[i,j] %显示相关的随机变量对的编号
有限差分和梯度
在 MATLAB7.0 中使用函数 diff()来计算差分，其调用格式如下。
实例：利用有限差分计算正弦函数的导数，代码设置如下：
%sin_diff.m %正弦函数的导数
x=0:0.1:10; %自变量
y=sin(x); %正弦函数
y_der=diff(y)./diff(x); %导数等于有限差分之比(dy/dx)
hold on;
x_der=x(1:(end-1)); %导数的自变量
plot(x,y,'b-'); %画图
plot(x_der,y_der,'b-.');
axis([0 10 -1.2 1.4]); %设定坐标轴范围
legend('正弦函数','正弦函数的导数'); %添加图例
在 MATLAB 7.0 中，用函数 gradient()来计算梯度
，计算二维高斯函数的梯度场，代码设置如下
%gradient_example.m %计算二维高斯函数的梯度场
v=-2:0.25:2;
[x,y]=meshgrid(v,v); %产生自变量 x,y
z= exp(-(x.^2+y.^2+0.5*x.*y)); %二维高斯函数
[px py]=gradient(z,0.25); %梯度场
contour(v,v,z,4); %画 4 条等高线
hold on;
quiver(v,v,px,py); %画出梯度场
图中椭圆是二维高斯函数的等高线，箭头的方 向代表梯度的方向，箭头的大小代表梯度的模。由图可见梯度方向总是垂直于等高线的，这 是梯度场的基本性质之一。
信号滤波和卷积
MATLAB 7.0 中有关信号滤波和卷积的函数，
函数名 功能描述 filter 一维数字滤波器 filter2 二维数字滤波器（将在本书的第 9 章详细介绍） conv 一维卷积和多项式乘法 conv2 二维卷积（将在本书的第 9 章详细介绍） convn N 维卷积 deconv 反卷积和多项式除法 detrend 去除信号中的直流或者线形成分，主要用于 FFT 运算中
一维数字滤波
MATLAB 7.0 的一维数字滤波可以用函数 filter()来实现，它是直接 II 型线性差分方程组 的解，可以用于实现 FIR 和 IIR 滤波
如何得到滤波器的系数 b 和 a 涉及到滤波器设计问题，具体将在本书的第 9 章作出详细 介绍。在有些情况下，可以很容易得到滤波器的系数 b 和 a。例如，一个 5 阶平均值滤波器： y(n)=1/5*x(n) + 1/5*x(n-1) + 1/5*x(n-2) + 1/5*x(n-3) + 1/5*x(n-4)，其系数 a 可以表示为 1，系
数 b 可以表示为[1/5 1/5 /15 1/5 1/5]。 下面通过一个示例来实现对带噪声的正弦信号进行平均值滤波，代码设置如下：
%filter_example.m %对带噪声的正弦信号进行平均值滤波
t=0:0.1:10; %时间
n = 6*randn(size(t)); %高斯白噪声
x = 40*sin(t)+n; %在正弦信号中添加噪声
a = 1; %频率值滤波器的系数
b = [1/5 1/5 1/5 1/5 1/5];
y=filter(b,a,x); %滤波
plot(t,x,'b-'); %画原始信号
hold on;
plot(t,y,'r:'); %画滤波后的信号
axis([0 10 -60 65]); %设置坐标轴范围
xlabel('时间/ts(s)');
legend('原始数据','滤波后数据'); %添加图例
．信号卷积
向量 u 和 v 的卷积为 w，则记为w u v = ⊗ 。如果 u，v 的长度分别为 m 和 n，则 w 的长 度为 m+n-1。
从信号中去除直流成分和线性成分，代码设置如下：
%detrend_example.m %从信号中去除直流成分和线性成分
t = 0:0.04:5; %时间
x = 2*t + 0.5*randn(size(t)); %带线性成分的随机信号
x_no_linear=detrend(x); %去除去除线性成分
subplot(2,1,1); %画图
hold on;
plot(t,x,'b-');
plot(t,x_no_linear,'b:'); axis([0 5 -2 14]); %设置坐标轴范围
title('从信号中去除线形成分');
legend('原始数据','去除线性成分的数据'); %添加图例
y = 3 + 0.5*randn(size(t)); %带直流成分的随机信号
y_no_constant = detrend(y,'constant'); %去除直流成分
subplot(2,1,2);
hold on;
plot(t,y,'b-');
plot(t,y_no_constant,'b:');
axis([0 5 -2 8]); %设置坐标轴范围
title('从信号中去直流成分');
legend('原始数据','去除直流成分的数据'); %添加图例
傅立叶变换
对信号进行分析通常可以在时域中进行，也可以在频域中进行，时域分析方法和频域分 析方法各有其优缺点。傅立叶变换是把信号从时域变换到频域，因此它在信号分析中占有极 其重要的地位，如滤波器设计、频谱分析等方面。 傅立叶变换既可以对连续信号进行分析，也可以对离散信号进行分析。连续信号的傅立
叶变换实际上要计算傅立叶积分，这部分内容可以参见本书第 7 章
MATLAB 7.0 提供的有关傅立叶变换的函数如表 6-5 所示。
函数名 功能描述 fft 一维离散快速傅立叶变换 fft2 二维离散快速傅立叶变换 fftn N 维离散快速傅立叶变换 ifft 一维离散快速傅立叶逆变换 ifft2 二维离散快速傅立叶逆变换 ifftn N 维离散快速傅立叶逆变换 fftshift 把频谱的直流成分移到中间 nextpow2 返回大于或等于参数值的二次幂，满足 2p>=abs(A)的最大整数
在 MATLAB 7.0 中一维离散傅立叶变换由函数 fft()来实现
在MATLAB 7.0中一维离散傅立叶逆变换由函数ifft()来实现，其调用方式基本与函数fft() 相同，只是添加一个选项，其调用方式如下。 • y = ifft(..., ’symmetric’)，把频谱 X(k)看成是共轭对称。该选项在 X(k)不是严格共轭对称 时有效； • y = ifft(..., ’nonsymmetric’)，该选项为默认选项
与傅立叶变换函数经常使用的函数是 fftshift()函数。如图 6-19 所示，信号 x(n)经过 FFT 后得到频谱 X(k)，X(k)的头部和尾部代表信号的低频分量大小，X(k)的中间则代表高频部分。 习惯上画频率响应曲线时把直流成分放在曲线的中间，而把高频部分放在曲线的头部和尾部。 因此频谱 X(k)不可以直接用于画频率响应曲线。为了画图时方便，MATLAB 7.0 提供函数 fftshift()用于把 FFT 变换的结果频谱 X(k)转换为适合画图显示的格式。对一维信号操作时， fftshift()函数把信号分为左半部分和右半部分，然后交换其左右部分。
，计算单位冲激信号经过一个带阻滤波器前后的频谱，代码设置如下：
%fliter_example.m %单位冲击信号通过带阻滤波器
t = 1:40; %信号的时间
x =zeros(size(t));
x(1) = 1; %产生单位冲击信号
>> [b,a] = butter(10,[0.3 0.7],'stop'); %设计带阻滤波器
y = filter(b,a,x); %滤波 hold on; %画图滤波前后的时域数据
stem(t,x,'marker','o');
stem(t,y,'marker','.');
xlabel('时间/ts(s)');
legend('原始数据','滤波后数据');
fx=fft(x); %对单位冲击信号进行傅立叶变换
fx=fftshift(fx); fy=fft(y); %滤波后的信号进行傅立叶变换
fy=fftshift(fy);
figure; subplot(2,1,1); %画图显示滤波前后的数据频谱
f=(t-20)/20; plot(f,abs(fx),'b-',f,abs(fy),'b-.'); %画信号的幅频曲线
xlabel('数字频率/\pi(rad)'); title('幅频曲线');
legend('原始数据的幅频曲线谱','滤波后数据的幅频曲线');
subplot(2,1,2);
plot(f,angle(fx),'b-',f,angle(fy),'b-.'); %画信号的相频曲线
xlabel('数字频率/\pi(rad)'); title('相频曲线'); legend('原始数据的相频曲线','滤波后数据的相频曲线');
由图 6-21 可见，单位冲激信号的频谱是一个常数，它通过一个带阻滤波器后的频谱等于 带通滤波器的频率响应曲线。这个结果验证了一个信号通过一个线性系统后的频谱等于信号 频谱乘以线性系统的频率响应的结论。
二维傅立叶变换
在图像处理中经常使用二维傅立叶变换来进行图像滤波等操作。在 MATLAB7.0 中二维 傅立叶变换用函数 fft2()来实现。函数 fft2()可以看成是对信号进行两次一维傅立叶变换，即 fft2(X)=fft(fft(X).’).’。函数 fft2()的调用格式如下：
实例：分析图 6-22 的频谱，代码设置如下
%fft2_example.m %分析图像的频谱
img = imread('coins.png'); %读图像文件
f_img = fft2(double(img)); %二维 fft 变换
f_img = fftshift(f_img); %把直流成分移动频谱的中间
imshow(img); %显示图像
figure; f_img_abs = abs(f_img); %得到频谱幅度
f_img_abs = (f_img_abs-min(min(f_img_abs)))./(max(max(f_img_abs))-min(min(f_img_abs)))*255; %把频谱幅度变换到[0,255]范围内
imshow(f_img_abs); %显示频谱幅度
title('图像的幅频分布');
功能函数 自写
功能函数就是可以用其他函数作为输入变量的函数，例如一个可以画函数图像的函数 fplot()。其调用格式为 fplot(@fun,[-pi,pi])。输入变量@fun 为需要画成图的函数的函数句柄。 MATLAB 7.0 的功能函数都在目录“\toolbox\matlab\funfun”下。
在 MATLAB 7.0 中函数可以用 M 文件、匿名函数和 inline 对象来表示。其中匿名函数是
MATLAB 7.0 的新特性。
%customized_fun.m %用户自己定义的函数
function y=customized_fun(x)
y=2./(1+exp(-x))+3./(1+exp(-2*x));
把 customized_fun.m 文件放在 MATLAB 7.0 的路径下，例如“\work”目录下，就可以使 用该函数了。例如在命令行中输入如下代码：
customized_fun(magic(4))
．一元函数最小值
一元函数在给定区间内的最小值问题可以用函数 fminbnd()来实现
实例：求正弦函数在[0,10]内的最小值，代码设置如下：
%fminbnd_example1.m %求正弦函数在[0 10]内的最小值
x=fminbnd(@sin,0,10)
如果要显示计算函数最小值的过程，则可以使用 options 参数来设定，代码设置如下
%fminbnd_example2.m %求正弦函数在[0 10]内的最小值
x=fminbnd(@sin,0,10,optimset(’Display’,’iter’))
求多元函数的最小值
求多元函数的最小值用函数 fminsearch()来实现。使用该函数时必须指定开始矢量 X0， 此函数返回一个在矢量 X0 附近的局部最小值
实例：求二维函数 2 2 2 2 1 1 f( ) 100( ) (1 ) x x x x = − + − 的局部最小值，代码设置如下
%fminsearch_example.m %寻找二维最小值
banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;
[x,fval] = fminsearch(banana,[-1.2, 1])
求一元函数的零点
可以使用函数 fzero()来求一元函数的零点。寻找一元函数零点时，可以指定一个开始点， 或者指定一个区间。当指定一个开始点时，此函数在开始点附近寻找一个使函数值变号的区 间，如果没有找到这样的区间，则函数返回 NaN。如果知道函数零点所在的区间，则可以使 用一个包含两个元素的向量来指定此区间。fzero()函数的调用格式如下。
用函数 fzero()来解一元非线性方程
2 2 1 1 1 ( 4) 1 ( 4) 1 2 x x + = + + − +
。该问题等价
于求函数
2 2 1 1 1 f( ) ( 4) 1 ( 4) 1 2 x x x = + + + − +
的零点，函数 f(x)的曲线图形如图 6-25 所示
%fzero_example.m %求解一元非线性方程
f = @(x) 1./((x+4).^2+1)+1./((x-4).^2+1)-0.5; %用匿名函数来表示函数
fplot(f,[-10 10]); %画出函数的图
xlabel('x');
ylabel('f(x)');
x1 = fzero(f,-5) %求某个点附近的零点
x2 = fzero(f,-3)
x3 = fzero(f,3)
x4 = fzero(f,5)
x1_region = fzero(f,[-10,-3]) %求某个区间内的零点
．一元函数的数值积分
MATLAB 7.0 提供了两个函数 quad()和 quadl()计算一元函数的积分。函数 quad()采用低 阶的自适应递归 Simpson 方法，函数 quadl()则采用高阶的自适应 Lobatto 方法
实例：，求归一化高斯函数的在区间[-1 1]上的定积分，代码设置如下：
%quad_exam.m %求归一化高斯函数的在区间[-1 1]上的定积分,并得到积分过程的中间结点 y=@(x)1/sqrt(pi)*exp(-x.^2); %归一化高斯函数 quad(y,-1,1,2e-6,1) %求定积分，并显示中间迭代过程 fplot(y,[-1 1],’b’); %画出函数
hold on; %跟踪数据（运行完上面程序后，可以在命令行中复制这些数据）
trace = [ 9 -1.0000000000 5.43160000e-001 0.1804679399;
11 -1.0000000000 2.71580000e-001 0.0728222057;
13 -0.7284200000 2.71580000e-001 0.1076454255;
15 -0.4568400000 9.13680000e-001 0.4817487615;
17 -0.4568400000 4.56840000e-001 0.2408826755;
19 -0.4568400000 2.28420000e-001 0.1142172651;
21 -0.2284200000 2.28420000e-001 0.1266655031;
23 0.0000000000 4.56840000e-001 0.2408826755;
25 0.0000000000 2.28420000e-001 0.1266655031;
27 0.2284200000 2.28420000e-001 0.1142172651;
29 0.4568400000 5.43160000e-001 0.1804679399;
31 0.4568400000&nb