曲线拟合与插值

新闻摘要:在大量的应用领域中，人们经常面临用一个解析函数描述数据（通常是测量值）的任务在插值法 里，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况人们设法找出某条光滑曲线，它 最佳地拟合数据，但不必要经过任何数据点标有'o'的是数据点在大量的应用领域中， 们经常面临用一个解析函数描述数据（通常是［无忧机械电子］测量值）的任务对这［无忧机摭电子］个问题有两种方法在插值法里，数据假定是正确的，要求以某种方法描述数据点 之间所发生的情况这种方法在下一节讨论这里讨论的方法是曲线拟合或回归人们设法找出某条光滑 曲线，它最佳地拟合数据，但不必要经过任何数据点图11.1说明了这两种方法标有'o'的是数据点； 连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合11.1 曲线拟合曲线拟合涉及回答两个基本问题：最佳拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义 最佳拟合，并存在无穷数目的曲线所以，从这里开始，我们走向何方？正如它证实的那样，当最佳拟合 被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的数学 上，称为多项式的最小二乘曲线拟合如果这种描述使你混淆，再研究图11.1。

虚线和标志的数据点之间 的垂直距离是在该点的误差对各数据点距离求平方，并把平方距离全加起来，就是误差平方和这条虚 线是使误差平方和尽可能小的曲线，即是最佳拟合最小二乘这个术语仅仅是使误差平方和最小的省略说 法图11.12阶曲线拟合在MATLAB中，函数polyfit求解最小二乘曲线拟合问题为了阐述这个函数的用法，让我们以上面图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];为了用polyfit,我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度如果我 们选择n=1作为阶次，得到最简单的线性近似通常称为线性回归相反，如果我们选择n=2作为阶次， 得到一个2阶多项式现在，我们选择一个2阶多项式» n=2;%polynomial order» p=polyfit(x, y, n)P =-9.810820.1293-0.0317polyfit的输出是一个多项式系数的行向量其解是y = —9.8108x^20.1293x—0.0317为了将曲线拟 合解与数据点比较，让我们把二者都绘成图。

» xi=linspace(0, 1, 100);%x-axis data for plotting» z=polyval(p, xi);为了计算在xi数据点的多项式值，调用MATLAB的函数polyval» plot(x, y, ' o ' , x, y, xi, z,':')画出了原始数据x和y,用'o'标出该数据点，在数据点之间，再用直线重画原始数据，并用点'：'线， 画出多项式数据xi和z» xlabel(' x '), ylabel(' y=f(x) '), title(' Second Order Curve Fitting ')将图作标志这些步骤的结果表示于前面的图11.1中多项式阶次的选择是有点任意的两点决定一直线或一阶多项式三点决定一个平方或2阶多项式按此 进行，n+1数据点唯一地确定n阶多项式于是，在上面的情况下，有11个数据点，我们可选一个高达1 0阶的多项式然而，高阶多项式给出很差的数值特性，人们不应选择比所需的阶次高的多项式此外， 随着多项式阶次的提高，近似变得不够光滑，因为较高阶次多项式在变零前，可多次求导例如，选一个 10阶多项式» pp=polyfit(x, y, 10);» format short e%change display format» pp.'%display polynomial coefficients as a column ans =-4.6436e+0052.2965e+006-4.8773e+0065.8233e+006-4.2948e+0062.0211e+006-6.0322e+0051.0896e+005-1.0626e+0044.3599e+002-4.4700e-001要注意在现在情况下，多项式系数的规模与前面的2阶拟合的比较。

还要注意在最小(-4.4700e-001)和最 大(5.8233e+006)系数之间有7个数量级的幅度差将这个解作图，并把此图与原始数据及2阶曲线拟合相 比较，结果如何呢？» zz=polyval(pp, xi);%evaluate 10th order polynomial» plot(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，显然，'越多就越好'的观念在这里不适用2nd and 10th Order curve Fitting141210它8t60 0.2 0.4 0.6 0.8 1图11.22阶和10阶曲线拟合11.2 一维插值正如在前一节对曲线拟合所描述的那样，插值定义为对数据点之间函数的估值方法，这些数据点是由某些 集合给定。

当人们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具例如，当数据点是 某些实验测量的结果或是过长的计算过程时，就有这种情况或许最简单插值的例子是MATLAB的作图按缺省，MATLAB用直线连接所用的数据点以作图这个线性插 值猜测中间值落在数据点之间的直线上当然，当数据点个数的增加和它们之间距离的减小时，线性插值 就更精确例如，» xl=linspace(0, 2*pi, 60);» x2=linspace(0, 2*pi, 6);» plot(x1, sin(x1), x2, sin(x2),' - ') » xlabel(' x '),ylabel(' sin(x) '),title(' Linear Interpolation ')Linear Interpolation0 1 2 3 4 5 6图11.3线性插值图11.3是sine函数的两个图，一个在数据点之间用60个点，它比另一个只用6个点更光滑和更精确 如曲线拟合一样，插值要作决策根据所作的假设，有多种插值而且，可以在一维以上空间中进行插值即如果有反映两个变量函数的插值，z=f(x, y),那么就可在x之间和在y之间，找出z的中间值进行插值。

MATLAB在一维函数interp1和在二维函数interp2中，提供了许多的插值选择其中的每个函数将在下面阐述为了说明一维插值，考虑下列问题，12小时内，一小时测量一次室外温度数据存储在两个MATLAB变量 中» hours=l:12;%index for hour data was recorded» temps二[5 89152529313022252724]; %recorded temperatures» plot(hours, temps, hours, temps,' + ' )%view temperatures» title(' Temperature ') » xlabel(' Hour ' ),ylabel(' Degrees Celsius ')Temperature35 1 1 1 30-+■10 - + -5 1 1 1 1 1 1 0 2 4 6 0 10 12Hour图11.4在线性插值下室外温度曲线人们可试正如图11.4看到的，MATLAB画出了数据点线性插值的直线为了计算在任意给定时间的温度, 着对可视的图作解释另外一种方法，可用函数interp1» t=interpl(hours, temps, 9.3)%estimate temperature at hour=9.3:000» t=interp1(hours, temps, 4.7)%estimate temperature at hour=4.7：2=» t=interp1(hours, temps, [3.26.57.111.7])%find temp at many points!t =10.200030.000030.900024.9000 interpl的缺省用法是由interpl(x, y, xo)来描述，这里x是独立变量(横坐标)，y是应变量(纵坐标)， xo是进行插值的一个数值数组。

另外，该缺省的使用假定为线性插值若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点最常用的方法是用一个3阶多项 式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的头两个导数与该数据点相一致 这种类型的插值被称为3次样条或简称为样条函数interpl也能执行3次样条插值» t=interpl(hours, temps, 9.3, ' spline ')%estimate temperature at hour=9.3t =21.8577» t=interp1(hours, temps, 4.7, ' spline ')%estimate temperature at hour=4.7t =22.3143» t=interpl(hours, temps, [3.26.57.111.7], ' spline ')t =9.673430.042731.175525.3820注意，样条插值得到的结果，与上面所示的线性插值的结果不同因为插值是一个估计或猜测的过程，其 意义在于，应用不同的估计规则导致不同的结果一个最常用的样条插值是对数据平滑也就是，给定一组数据，使用样条插值在更细的间隔求值。

例如，» h=l:0.1:12;%estimate temperature every 1/10 hour» t=interp1(hours, temps, h, ' spline ');» plot(hours, temps, ' - ' , hours, temps, ' + ' , h, t)%plot comparative results» title(' Springfield Temperature ') » xlabel(' Hour '),ylabel(' Degrees Celsius ')在图11.5中，虚线是线性插值，实线是平滑的样条插值，标有'+ '的是原始数据如要求在时间轴上有 更细的分辨率，并使用样条插值，我们有一个更平滑、但不一定更精确地对温度的估计尤其应注意，在 数据点，样条解的斜率不突然改变作为这个平滑插值的回报，3次样条插值要求更大量的计算，因为必 须找到3次多项式以描述给定数据之间的特征关于样条的更详细信息可见下一章Springfield Temperature35 1 1 1 130252010 - -XE I I I I I 0 2 4 6 8 10 12Hour图11.5在不同插值下室外温度曲线在讨论二维插值之前，了解interpl所强制的二个强约束是很重要的。

首先，人们不能要求有独立变量范 围以外的结果，例如，interp1(hours, temps, 13.5)导致一个错误，因为hours在1到12之间变化其 次，独立变量必须是单调的即独立变量在值上必须总是增加的或总是减小的在我们的例子里，hours 是单调的然而，如果我们已经定义独立变量为一天的实际时间，» time_of_day=[7:121:6]%start at 7AM,end at 6PMtime_of_day =789101112123456则独立变量将不是单调的，因为time_of_day增加到12,然后跌到1,再然后增加如果用time_of_day代替interp1中的hours，将会返回一个错误同样的理由, 人们不能对temps插值来找出产生某温度的时间(小时)，因为temps不是单调的11.3 二维插值 二维插值是基于与一维插值同样的基本思想然而，正如名字所隐含的，二维插值是对两变量的函数Z=f(x, y)进行插值为了说明这个附加的维数，考虑一个问题设人们对平板上的温度分布估计感兴趣，给定的 温度值取自平板表面均匀分布的格栅采集了下列的数据：» width=l:5;%index for width of plate (i.e.,the x-dimension)» depth=1:3;%index for depth of plate (i,e,,the y-dimension)» temps=[8281808284; 7963616581; 8484828586]%temperature datatemps =828180828479636165818484828586如同在标引点上测量一样，矩阵temps表示整个平板的温度分布。

temps的列与下标dep th或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 ') ;%cubic interpolation» plot(wi, zlinear, ' - ' , wi, zcubic)%plot results» xlabel(' Width of Plate '),ylabel(' Degrees Celsius ')» title( [' Temperature at Depth ='num2str(d)])另一种方法，我们可以在两个方向插值先在三维坐标画出原始数据，看一下该数据的粗糙程度(见图11. 7)» mesh(width, depth, temps)%use mesh plot» xlabel(' Width of Plate '),ylabel(' Depth of Plate ') » zlabel(' Degrees Celsius '),axis(' ij '),grid434Temperature at Depth = 280End-20岳 auEOJcl656075702 3 4Width of Plate图11.6在深度d=2处的平板温度畧一空OJOEaJQJ-JEiQJclDepth of Plate2Width of Plate图11.7平板温度然后在两个方向上插值，以平滑数据。

» di=l:0.2:3;%choose higher resolution for depth» wi=l:0.2:5;%choose higher resolution for width» zcubic=interp2(width, depth, temps, wi, di, ' cubic ') ;%cubic» mesh(wi, di, zcubic) » xlabel( Width of Plate ),ylabel( Depth of Plate )» zlabel(' Degrees Celsius '),axis(' ij '),grid上面的例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持跟踪interp2的基本形式是i nt erp2(x, y, z, xi, yi, met hod)这里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二维插值后的平板温度可选的参数met hod可以是'linear'，'cubic'或'neares t'在这种情况下，cubic不意味着3次样条，而 是使用3次多项式的另一种算法linear方法是线性插值，仅用作连接图上数据点nearest方法只选择 最接近各估计点的粗略数据点在所有的情况下，假定独立变量x和y是线性间隔和单调的关于这些方 法的更多的信息，可请求在线帮助，例如，》help interp2,或参阅MATLAB参考手册11・4 M文件举例虽然对于许多应用，函数interp1和interp2是很有用的，但它们限制为对单调向量进行插值在某些情 况，这个限制太严格例如，考虑下面的插值：10.80.60.40.212345yi =» x=linspace(0, 5);» y=l-exp(-x).*sin(2*pi*x);» plot(x, y)图 11.9 函数 l-exp(-x).*sin(2*pi*x)的曲线函数interpl可用来在任何值或x的值上估计y值» yi=interp1(x, y, 1.8)1.1556然而，interp1不能找出对应于某些y值的x值例如，如在图11.9上所示，考虑寻找y=1.1处的x值:10.80.60.40.212345Xi =图11.10给y值在函数曲线上求x的值» plot(x, y, [0, 5], [1.11.1])从图11.10上，我们看到有四个交点。

使用interpl,我们得到:» xi=interp1(y, x, 1.1)??? Error using ==> table1First column of the table must be monotonic.这个函数interp1失败，由于y不是单调的在本章精通MATLAB工具箱所说明的M文件例子，消除了单调性的要求» table=[x; y].' ;%create column oriented table from data» xi=mminterp(table, 2, 1.1)0.52811.10000.95801.10001.58251.10001.88471.1000这里使用了线性插值，函数mminterp估计了 y=1.1处的四个点由于函数mminterp的一般性质，要插值 的数据是由面向列矩阵给出，在上面的例子中称作为表(table)第二个输入参量是被搜索矩阵table的列, 第三个参量是要找的值这个精通MATLAB工具箱函数的主体由下面给出：function y=mminterp(tab, col, val) %MMINTERP 1-D Table Search by Linear Interpolation. %Y=MMINTERP(TAB,COL,VAL) linearly interpolates the table %TAB searching for the scalar value VAL in the column COL.%All crossings are found and TAB(:,COL) need not be monotonic.%Each crossing is returned as a separate row in Y and Y has as%many columns as TAB.Naturally,the column COL of Y contains%the value VAL. If VAL is not found in the table,Y=[]. %Copyright (c) 1996 by Prentice-Hall,Inc.[rt, ct]=size(tab);if length(val) > 1,error(' VAL must be a scalar. '),endif col>ct|col < 1,error(' Chosen column outside table width. '),end if rt < 2,error(' Table too small or not oriented in columns. '),end above=tab(: , col) > val;%True where > VALbelow=tab(: , col) < val;%True where < VALequal=tab(: , col) = = val;%True where = VALif all(abol(below = = 0),%handle simplest case y=tab(find(equal), : );returnendpslope=find(below(1:rt-1)&above(2:rt));%indices where slope is + nslope=find(below(2:rt )& above(1:rtT));%indices where slope is 一ib=sort([pslope; nslope+l]);%put indices below in order ia=sort([nslope; pslope+l]);%put indices above in order ie=find(equal);%indices where equal to val[tmp,ix]=sort( [ib, ie] );%find where equals fit in resultieq=ix > length(ib);%True where equals values fitry=length(tmp);%# of rows in result yy=zeros(ry, ct);%poke data into a zero matrix alpha=(val-tab(ib,col))./(tab(ia,col)-tab(ib,col));alpha=alpha(: , ones(1, ct));%duplicate for all columnsy(~ieq, : )=alpha.*tab(ia, : ) + (l-alpha).*tab(ib, : );%interpolated valuesy(ieq, : )=tab(ie, : );%equal valuesy( : , col)=val*ones(ry, 1);%remove roundoff error正如所见的，mminterp利用了 find和sort函数、逻辑数组和数组操作技术。

没有For循环和While循环不论用其中哪一种技术来实现将使运行变慢，尤其对大的表注意mminterp与含有大于或等于2的任意数列的表一起工作，如同函数interpl 一样而且，在这种情况下，插值变量可以是任意的列例如，» z=sin(pi*x);%add more data to table» table=[x; y; z].';» t=mminterp(table, 2, 1.1)%same interpolation as earliert =0.52811.10000.99300.95801.10000.13141.58251.1000- 0.96391.88471.1000- 0.3533ve = = 0) | al » t=mminterp(table, 3, -.5)%second third column now1.16690.7316-0.50001.83291.1377-0.50003.16710.9639-0.50003.83311.0187-0.5000这些最后的结果估计了 x和y在z= -.5处的值尽管逐条地对函数mminterp解释如何工作是很有帮助的，但这样做要求有更多的篇幅和时间。

解释mmint erp如何工作最容易的方法是创建一个小表格，然后，在重要的语句末尾删除分号以后，调用函数这样, 中间值将帮助用户理解函数是如何找到与所需值相符的数据值以及如何执行插值前面已阐述了 int erp1的用法当用于线性插值时，只要所要求的插值点的个数少，in terpl工作很好 在要求许多插值点情况下，由于所用的算法，interp1工作较慢为了克服这个问题，精通MATLAB工具箱 包括了函数mmtable，它的帮助文本是：»help tableMMTABLE 1-D Monotonic Table Search by Linear Interpolation.YI=MMTABLE(TAB,COL,VALS) linearly interpolates the table TABsearching for values VALS in the column COL.TAB(:,COL) must be monotonic, but need NOT be equally spaced.YI has as many rows as VALS and as many columns TABNaNs are returned where VALS are outside the range of TAB(:,COL).YI=MMTABLE(TAB,VALS) interpolates using COL=1 and does not returnTAB(:,1) in Y.This matches the usage of TABLE1(TAB,X0).YI=MMTABLE(X,Y,XI) interpolates the vector X to find YI associatedwith XI. This match the usage of INTERP1(X,Y,XI)This routine is 10X faster than TABLE1which is called by INTERP1.MMTABLE由线性插值实现一维单调表搜索YI=MMTABLE(TAB,COL,VALS)线性地对表TAB进行插值，在列COL中搜索值为VALSTAB (:, COL)必须是单调的，但不必等价地生成空间。

YI与VALS有同样的行和与TAB有同样的列当VALS超出TAB (:, COL)的范围，返回NaNs.YI=MMTABLE(TAB,VALS)使用COL=1进行插值，不返回在Y中的TAB(:,1)这和TABLE1(TAB,XO)的用法匹配YI=MMTABLE(X,Y,XI)为了找出YI和XI的关系，对向量X进行插值这和INTERP1(X,Y,XI)的用法匹配这个例程比由INTERP1调用TABLE1快10倍正如前面描述的，可以用几种方式调用mmtable此外，要插值的列或向量不需要线性间隔由于这个原 因，mmtable比ilinear函数更普遍在MATLAB版本5中，interp1将用ilinear来实现线性插值11.5小结下面的表11.1总结了在MATLAB中所具有的曲线拟合和插值函数表 11.1曲线拟合和插值函数polyfit(x, y, n)对描述n阶多项式y=f(x)的数据进行最小二乘曲线拟合interp1(x, y, xo)1维线性插值interp1(x, y, xo, ' spline ')1维3次样条插值interp1(x, y, xo, ' cubic ') 1维3次插值interp2(x, y, Z, xi, yi)2维线性插值interp2(x, y, Z, xi, yi, ' cubic ')2维3次插值interp2(x, y, Z, xi, yi, ' nearest ')2维最近邻插值。