在大量的应用领域中，人们经常面临用一个解析函 数描述数据(通常是测量值)的任务。对这个问题有两种方法。在插 值法里，数据假定是正确的，要求以某种方法描述数据点之间所发生 的情况。这种方法在下一节讨论。这里讨论的方法是曲线拟合或回归 。人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任 何数据点。图11.1说明了这两种方法。标有'o'的是数据点； 连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合。
11.1 曲线拟合
曲线拟合涉及回答两个基本问题：最 佳拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义 最佳拟合，并存在无穷数目的曲线。所以，从这里开始，我们走向何 方？正如它证实的那样，当最佳拟合被解释为在数据点的最小误差平 方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。 数学上，称为多项式的最小二乘曲线拟合。如果这种描述使你混淆， 再研究图11.1。虚线和标志的数据点之间的垂直距离是在该点的 误差。对各数据点距离求平方，并把平方距离全加起来，就是误差平 方和。这条虚线是使误差平方和尽可能小的曲线，即是最佳拟合。最 小二乘这个术语仅仅是使误差平方和最小的省略说法。
< br>图11.12阶曲线拟合
在MATLAB中，函数p olyfit求解最小二乘曲线拟合问题。为了阐述这个函数的用法 ，让我们以上面图11.1中的数据开始。
»
x=[0.1.2.3.4.5.6.7.8.91];
»
y=[-.4471.9783.286.167 .087.347.669.569.489.3011.2];
< br>为了用polyfit，我们必须给函数赋予上面的数据和我 们希望最佳拟合数据的多项式的阶次或度。如果我们选择n=1作为 阶次，得到最简单的线性近似。通常称为线性回归。相反，如果我们 选择n=2作为阶次，得到一个2阶多项式。现在，我们选择一个2 阶多项式。
»
n=2;
%polynomi al order
»
p=polyfit( x, y, n)
p =
-9.810820.1 293-0.0317
polyfit 的输出是一个多项 式系数的行向量。其解是y = －9.8108x2 ＋20.1 293x－0.0317。为了将曲线拟合解与数据点比较，让我们 把二者都绘成图。
»
xi=linspac e(0, 1, 100);
%x-axis data for plotting
»
z=polyval( p, xi);
为了计算在xi数据点的多项式值，调用M ATLAB的函数polyval。
»
pl ot(x, y, ' o ' , x, y, xi, z, ' : ' )
画出了原始数据x和y，用'o'标出该数 据点，在数据点之间，再用直线重画原始数据，并用点' : '线 ，画出多项式数据xi和z。
»
xlabe l(' x '), ylabel(' y=f(x) '), title(' Second Order Curve Fit ting ')
将图作标志。这些步骤的结果表示于前面的 图11.1中。
多项式阶次的选择是有点任意的。两点决定 一直线或一阶多项式。三点决定一个平方或2阶多项式。按此进行， n+1数据点唯一地确定n阶多项式。于是，在上面的情况下，有1 1个数据点，我们可选一个高达10阶的多项式。然而，高阶多项式 给出很差的数值特性，人们不应选择比所需的阶次高的多项式。此外 ，随着多项式阶次的提高，近似变得不够光滑，因为较高阶次多项式 在变零前，可多次求导。例如，选一个10阶多项式
&ra quo;
pp=polyfit(x, y, 10) ;
»
format short e%change display format
»
pp. '%display polynomial coefficie nts as a column
ans =
-4 .6436e+005
2.2965e+006
- 4.8773e+006
5.8233e+006
-4.2948e+006
2.0211e+006
-6.0322e+005
1.0896e+005-1.0626e+004
4.3599e+002< br>-4.4700e-001
要注意在现在情况下，多 项式系数的规模与前面的2阶拟合的比较。还要注意在最小(-4. 4700e-001)和最大(5.8233e+006)系数之间 有7个数量级的幅度差。将这个解作图，并把此图与原始数据及2阶 曲线拟合相比较，结果如何呢？
»
zz=p olyval(pp, xi);
%evaluate 10th order polynomial
»
pl ot(x, y, ' o ' , xi, z, ' : ' , xi, zz)%plot data
»
xlabel(' x '),ylabel(' y=f(x) '),title(' 2nd and 10th Order curve Fitting ')
在下面的图11.2 中，原始数据标以'o'，2阶曲线拟合是虚线，10阶拟合是实线 。注意，在10阶拟合中，在左边和右边的极值处，数据点之间出现 大的纹波。当企图进行高阶曲线拟合时，这种纹波现象经常发生。根 据图11.2，显然，‘ 越多就越好 ’的观念在这里不适用。< br>
图11.22阶和10阶曲线拟合
11.2 一维插值
正如在前一节对曲线拟合所描述的那样 ，插值定义为对数据点之间函数的估值方法，这些数据点是由某些集 合给定。当人们不能很快地求出所需中间点的函数值时，插值是一个 有价值的工具。例如，当数据点是某些实验测量的结果或是过长的计 算过程时，就有这种情况。
或许最简单插值的例子是MAT LAB的作图。按缺省，MATLAB用直线连接所用的数据点以作 图。这个线性插值猜测中间值落在数据点之间的直线上。当然，当数 据点个数的增加和它们之间距离的减小时，线性插值就更精确。例如 ，
»
x1=linspace(0, 2* pi, 60);
»
x2=linspac e(0, 2*pi, 6);
»
plot (x1, sin(x1), x2, sin(x2), ' - ')
»
xlabel(' x '),y label(' sin(x) '),title(' Line ar Interpolation ')
图1 1.3线性插值
图11.3是sine函数的两个图，一个 在数据点之间用60个点，它比另一个只用6个点更光滑和更精确。
如曲线拟合一样，插值要作决策。根据所作的假设，有多种 插值。而且，可以在一维以上空间中进行插值。即如果有反映两个变 量函数的插值，z=f(x, y)，那么就可在x之间和在y之间 ，找出z的中间值进行插值。MATLAB在一维函数interp 1和在二维函数interp2中，提供了许多的插值选择。其中的 每个函数将在下面阐述。
为了说明一维插值，考虑下列问题 ，12小时内，一小时测量一次室外温度。数据存储在两个MATL AB变量中。
»
hours=1:12;
% index for hour data was record ed
»
temps=[5 8915252 9313022252724];
%recorded temp eratures
»
plot(hours , temps, hours, temps,' + ')%v iew temperatures
»
ti tle(' Temperature ')
» ;
xlabel(' Hour '),ylabel(' De grees Celsius ')
图11.4 在线性插值下室外温度曲线
正如图11.4看到的，MAT LAB画出了数据点线性插值的直线。为了计算在任意给定时间的温 度，人们可试着对可视的图作解释。另外一种方法，可用函数int erp1。
»
t=interp1(hou rs, temps, 9.3)%estimate tempe rature at hour=9.3
t =
2 2.9000
»
t=interp1(ho urs, temps, 4.7)%estimate temp erature at hour=4.7
t =
22
»
t=interp1(hours, temps, [3.26.57.111.7])%find temp at many points!
t =
10.2000
30.0000
30.9000
24.9000
interp1的缺省用法是由i nterp1(x, y, xo)来描述，这里x是独立变量(横 坐标)，y是应变量(纵坐标)，xo是进行插值的一个数值数组。 另外，该缺省的使用假定为线性插值。
若不采用直线连接数 据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法是 用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模 ，每个3次多项式的头两个导数与该数据点相一致。这种类型的插值 被称为3次样条或简称为样条。函数interp1也能执行3次样 条插值。
»
t=interp1(hour s, temps, 9.3, ' spline ')%est imate temperature at hour=9.3< br>t =
21.8577
»
t =interp1(hours, temps, 4.7, ' spline ')%estimate temperature at hour=4.7
t =
22.3143
»
t=interp1(hours, t emps, [3.26.57.111.7], ' splin e ')
t =
9.6734
30.0 427
31.1755
25.3820
注 意，样条插值得到的结果，与上面所示的线性插值的结果不同。因为 插值是一个估计或猜测的过程，其意义在于，应用不同的估计规则导 致不同的结果。
一个最常用的样条插值是对数据平滑。也就 是，给定一组数据，使用样条插值在更细的间隔求值。例如，
»
h=1:0.1:12;
%estimate temperature every 1/10 hour
»
t=interp1(hours, temp s, h, ' spline ') ;
»
plot(hours, temps, ' - ' , ho urs, temps, ' + ' , h, t)%plot comparative results
» ;
title(' Springfield Temperat ure ')
»
xlabel(' Hou r '),ylabel(' Degrees Celsius ')
在图11.5中，虚线是线性插值，实线是平滑的样条 插值，标有' + '的是原始数据。如要求在时间轴上有更细的分 辨率，并使用样条插值，我们有一个更平滑、但不一定更精确地对温 度的估计。尤其应注意，在数据点，样条解的斜率不突然改变。作为 这个平滑插值的回报，3次样条插值要求更大量的计算，因为必须找 到3次多项式以描述给定数据之间的特征。关于样条的更详细信息可 见下一章。
图11.5在不同插值下室外温度曲 线
在讨论二维插值之前，了解interp1所强制的二个 强约束是很重要的。首先，人们不能要求有独立变量范围以外的结果 ，例如，interp1(hours, temps, 13.5 )导致一个错误，因为hours在1到12之间变化。其次，独立 变量必须是单调的。即独立变量在值上必须总是增加的或总是减小的 。在我们的例子里，hours是单调的。然而，如果我们已经定义 独立变量为一天的实际时间，
»
time_ of_day=[7:121:6]%start at 7AM, end at 6PM
time_of_day =
789101112123456
则独立变量将不是单调 的，因为time_of_day增加到12，然后跌到1，再然后 增加。如果用time_of_day代替interp1中的ho urs，将会返回一个错误。同样的理由，人们不能对temps插 值来找出产生某温度的时间(小时)，因为temps不是单调的。
11.3 二维插值
二维插值是基于与一维插值同 样的基本思想。然而，正如名字所隐含的，二维插值是对两变量的函 数z=f(x, y)进行插值。为了说明这个附加的维数，考虑一 个问题。设人们对平板上的温度分布估计感兴趣，给定的温度值取自 平板表面均匀分布的格栅。
采集了下列的数据：
& raquo;
width=1:5;
%index for wi dth of plate (i.e.,the x-dimen sion)
»
depth=1:3;
%in dex for depth of plate (i,e,,t he y-dimension)
»
tem ps=[8281808284;
7963616581;
84 84828586]%temperature data
temps =
8281808284
79636 16581
8484828586
如同在标引点上 测量一样，矩阵temps表示整个平板的温度分布。temps的 列与下标depth或y-维相联系，行与下标width或x-维 相联系(见图11.6)。为了估计在中间点的温度，我们必须对它 们进行辨识。
»
wi=1:0.2:5;
% estimate across width of plate
»
d=2;
%at a depth of 2
»
zlinear=interp2( width, depth, temps, wi, d) ;
% linear interpolation
» ;
zcubic=interp2(width, depth, temps, wi,d, ' cubic ') ;
%cub ic interpolation
»
pl ot(wi, zlinear, ' - ' , wi, zc ubic)%plot results
»
xlabel(' Width of Plate '),yla bel(' Degrees Celsius ')
&r aquo;
title( [' Temperature at Depth ='num2str(d) ] )
另一种 方法，我们可以在两个方向插值。先在三维坐标画出原始数据，看一 下该数据的粗糙程度(见图11.7)。
»
mesh(width, depth, temps)%use mesh plot
»
xlabel(' Width of Plate '),ylabel(' Dep th of Plate ')
»
zlab el(' Degrees Celsius '),axis(' ij '),grid
图11.6在 深度d=2处的平板温度
图11.7平 板温度
然后在两个方向上插值，以平滑数据。
&r aquo;
di=1:0.2:3;
%choose highe r resolution for depth
&raq uo;
wi=1:0.2:5;
%choose higher resolution for width
» ;
zcubic=interp2(width, depth, temps, wi, di, ' cubic ') ;
%c ubic
»
mesh(wi, di, z cubic)
»
xlabel(' Wid th of Plate '),ylabel(' Depth of Plate ')
»
zlabel( ' Degrees Celsius '),axis(' ij '),grid
上面的例子清楚地证明了，二维插值更为 复杂，只是因为有更多的量要保持跟踪。interp2的基本形式 是interp2(x, y, z, xi, yi, meth od)。这里x和y是两个独立变量，z是一个应变量矩阵。x和y 对z的关系是
z(i, :) = f(x, y(i)) 和z(:, j) = f(x(j), y).
也就是， 当x变化时，z的第i行与y的第i个元素y(i)相关，当y变化 时，z的第j列与x的第j个元素x(j)相关，。xi是沿x-轴 插值的一个数值数组；yi是沿y-轴插值的一个数值数组。
图11.8二维插值后的平板温度
可选的参数 method可以是 'linear'，'cubic'或'ne arest'。在这种情况下，cubic不意味着3次样条，而是 使用3次多项式的另一种算法。linear方法是线性插值，仅用 作连接图上数据点。nearest方法只选择最接近各估计点的粗 略数据点。在所有的情况下，假定独立变量x和y是线性间隔和单调 的。关于这些方法的更多的信息，可请求在线帮助，例如，&raq uo;
help interp2，或参阅MATLAB参考手册 。
11.4 M文件举例
虽然对于许多应用，函数 interp1和interp2是很有用的，但它们限制为对单调 向量进行插值。在某些情况，这个限制太严格。例如，考虑下面的插 值：
»
x=linspace(0, 5) ;
»
y=1-exp(-x).*sin( 2*pi*x);
»
plot(x, y)
图11.9函数1-exp(-x). *sin(2*pi*x)的曲线
函数interp1可用 来在任何值或x的值上估计y值。
»
yi= interp1(x, y, 1.8)
yi =
1.1556
然而，interp1不能找出对应于某些y 值的x值。例如，如在图11.9上所示，考虑寻找y=1.1处的 x值：
图11.10给y值在函数曲线上求x的 值
»
plot(x, y, [0, 5] , [1.11.1] )
从图11.10上，我们看到有 四个交点。使用interp1，我们得到：
» ;
xi=interp1(y, x, 1.1)
??? Error using ==>
table1
First column of the table must be monotonic.
这个函数interp1失 败，由于y不是单调的。
在本章精通MATLAB工具箱所 说明的M文件例子，消除了单调性的要求。
»
table=[x;
y].' ;
%create colum n oriented table from data
»
xi=mminterp(table, 2, 1.1)
xi =
0.52811.10000.95801.1000
1.58251.1000
1.88471.1000
这里使用了线性插值， 函数mminterp估计了y=1.1处的四个点。由于函数mm interp的一般性质，要插值的数据是由面向列矩阵给出，在上 面的例子中称作为表(table)。第二个输入参量是被搜索矩阵 table的列，第三个参量是要找的值。
这个精通MAT LAB工具箱函数的主体由下面给出：
function y=mminterp(tab, col, val)
% MMINTERP 1-D Table Search by L inear Interpolation.
%Y=MMI NTERP(TAB,COL,VAL) linearly in terpolates the table
%TAB s earching for the scalar value VAL in the column COL.
%All crossings are found and TAB(: ,COL) need not be monotonic.%Each crossing is returned a s a separate row in Y and Y ha s as
%many columns as TAB.N aturally,the column COL of Y c ontains
%the value VAL. If VAL is not found in the table, Y=[].
%Copyright (c) 1996 b y Prentice-Hall,Inc.
[rt, c t]=size(tab);
if length(val ) >
1,error(' VAL must be a scalar. '),end
if col>
c t|col <
1,error(' Chosen co lumn outside table width. '),e nd
if rt <
2,error(' Tab le too small or not oriented i n columns. '),end
above=tab (: , col) >
val;
%True where >
VAL
below=tab(: , col ) <
val;
%True where <
VA L
equal=tab(: , col) = = va l;
%True where = VAL
if all( above = = 0) | all(below = = 0 ),%handle simplest case
y=t ab(find(equal), : );
return
end
pslope=find(below(1:rt- 1)&
above(2:rt));
%indices w here slope is +
nslope=find (below(2:rt)&
above(1:rt-1) );
%indices where slope is -
ib=sort([pslope;
nslope+1]);
% put indices below in order
ia=sort([nslope;
pslope+1]);
%p ut indices above in order
i e=find(equal);
%indices where e qual to val
[tmp,ix]=sort( [ib, ie] );
%find where equals fit in result
ieq=ix >
l ength(ib);
%True where equals v alues fit
ry=length(tmp);
% # of rows in result y
y=zer os(ry, ct);
%poke data into a z ero matrix
alpha=(val-tab(i b,col))./(tab(ia,col)-tab(ib,c ol));
alpha=alpha(: , ones( 1, ct));
%duplicate for all col umns
y(~ieq, : )=alpha.*tab (ia, : )+(1-alpha).*tab(ib, : );
%interpolated values
y(ie q, : )=tab(ie, : );
%equal valu es
y( : , col)=val*ones(ry, 1);
%remove roundoff error
正如所见的，mminterp利用了find和sort 函数、逻辑数组和数组操作技术。没有For循环和While循环 。不论用其中哪一种技术来实现将使运行变慢，尤其对大的表。注意 mminterp与含有大于或等于2的任意数列的表一起工作，如 同函数interp1一样。而且，在这种情况下，插值变量可以是 任意的列。例如，
»
z=sin(pi*x );
%add more data to table
& raquo;
table=[x;
y;
z].' ;
»
t=mminterp(table, 2, 1 .1)%same interpolation as earl ier
t =
0.52811.10000.99 30
0.95801.10000.1314
1. 58251.1000-0.9639
1.88471.1 000-0.3533
»
t=mminte rp(table, 3, -.5)%second third column now
t =
1.16690. 7316-0.5000
1.83291.1377-0. 5000
3.16710.9639-0.5000
3.83311.0187-0.5000
这些最后的结 果估计了x和y在z= -.5处的值。
尽管逐条地对函数 mminterp解释如何工作是很有帮助的，但这样做要求有更多 的篇幅和时间。解释mminterp如何工作最容易的方法是创建 一个小表格，然后，在重要的语句末尾删除分号以后，调用函数。这 样，中间值将帮助用户理解函数是如何找到与所需值相符的数据值以 及如何执行插值。
前面已阐述了interp1的用法。当 用于线性插值时，只要所要求的插值点的个数少，interp1工 作很好。在要求许多插值点情况下，由于所用的算法，interp 1工作较慢。为了克服这个问题，精通MATLAB工具箱包括了函 数mmtable，它的帮助文本是：
»
he lp table
MMTABLE 1-D Monoto nic Table Search by Linear Int erpolation.
YI=MMTABLE(TAB, COL,VALS) linearly interpolate s the table TAB
searching f or values VALS in the column C OL.
TAB(:,COL) must be mono tonic, but need NOT be equally spaced.
YI has as many row s as VALS and as many columns TAB
NaNs are returned where VALS are outside the range of TAB(:,COL).
YI=MMTABLE(TAB ,VALS) interpolates using COL= 1 and does not return
TAB( :,1) in Y.This matches the usa ge of TABLE1(TAB,X0).
YI=MM TABLE(X,Y,XI) interpolates the vector X to find YI associate d
with XI. This match the usage of INTERP1(X,Y,XI)
Th is routine is 10X faster than TABLE1which is called by INTER P1.
MMTABLE由线性插值实现一维单调表搜索YI=MMTABLE(TAB,COL,VALS) 线性地 对表TAB进行插值，在列COL中搜索值为VALS
TA B(:,COL)必须是单调的，但不必等价地生成空间。
YI与VALS有同样的行和与TAB有同样的列。
当VA LS超出TAB(:,COL)的范围，返回NaNs.
Y I=MMTABLE(TAB,VALS) 使用COL=1进行插 值，不返回在Y中的TAB(:,1)
这和TABLE1( TAB,X0)的用法匹配。
YI=MMTABLE(X, Y,XI) 为了找出YI和XI的关系，对向量X进行插值。这和INTERP1(X,Y,XI)的用法匹配。
这 个例程比由INTERP1调用TABLE1快10倍。
正 如前面描述的，可以用几种方式调用mmtable。此外，要插值 的列或向量不需要线性间隔。由于这个原因，mmtable比il inear函数更普遍。在MATLAB版本5中，interp1 将用ilinear来实现线性插值。
11.5 小结下面的表11.1总结了在MATLAB中所具有的曲线拟合和 插值函数。
表11.1
曲 线 拟 合 和 插 值 函 数
polyfit(x, y, n)< br> 对描述n阶多项式y=f(x)的数据 进行最小二乘曲线拟合
inte rp1(x, y, xo)
1维线性插 值
interp1(x, y, xo, ' spline ')
1维3次样条插值
interp 1(x, y, xo, ' cubic ')
1维3次插值
interp2(x, y, Z, x i, yi)
2维线性插值
inte rp2(x, y, Z, xi, yi, ' cubic ' )
2维3次插值
< br>
interp2(x , y, Z, xi, yi, ' nearest ') 2维最近邻插值
|