第11章 曲线拟合与插值 在大量的应用领域中，人们经常面临用一个解析函数描述数据(通常是测量值)的任务。对这个问题有两种方法。在插值法里，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。这种方法在下一节讨论。这里讨论的方法是曲线拟合或回归。人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。图11.1说明了这两种方法。标有'o'的是数据点；连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合。 11.1 曲线拟合 曲线拟合涉及回答两个基本问题：最佳拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义最佳拟合，并存在无穷数目的曲线。所以，从这里开始，我们走向何方？正如它证实的那样，当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。数学上，称为多项式的最小二乘曲线拟合。如果这种描述使你混淆，再研究图11.1。虚线和标志的数据点之间的垂直距离是在该点的误差。对各数据点距离求平方，并把平方距离全加起来，就是误差平方和。这条虚线是使误差平方和尽可能小的曲线，即是最佳拟合。最小二乘这个术语仅仅是使误差平方和最小的省略说法。  图11.1 2阶曲线拟合 在MATLAB中，函数polyfit求解最小二乘曲线拟合问题。为了阐述这个函数的用法，让我们以上面图11.1中的数据开始。 ? x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1]; ? y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; 为了用polyfit，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度。如果我们选择n=1作为阶次，得到最简单的线性近似。通常称为线性回归。相反，如果我们选择n=2作为阶次，得到一个2阶多项式。现在，我们选择一个2阶多项式。 ? n=2; % polynomial order ? p=polyfit(x, y, n) p = -9.8108 20.1293 -0.0317 polyfit 的输出是一个多项式系数的行向量。其解是y = －9.8108x2 ＋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个数据点，我们可选一个高达10阶的多项式。然而，高阶多项式给出很差的数值特性，人们不应选择比所需的阶次高的多项式。此外，随着多项式阶次的提高，近似变得不够光滑，因为较高阶次多项式在变零前，可多次求导。例如，选一个10阶多项式 ? pp=polyfit(x, y, 10) ; ? format short e % change display format ? pp.' % display polynomial coefficients 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 -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，显然，‘ 越多就越好 ’的观念在这里不适用。  图11.2 2阶和10阶曲线拟合 11.2 一维插值 正如在前一节对曲线拟合所描述的那样，插值定义为对数据点之间函数的估值方法，这些数据点是由某些集合给定。当人们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况。 或许最简单插值的例子是MATLAB的作图。按缺省，MATLAB用直线连接所用的数据点以作图。这个线性插值猜测中间值落在数据点之间的直线上。当然，当数据点个数的增加和它们之间距离的减小时，线性插值就更精确。例如， ? x1=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 ')  图11.3 线性插值 图11.3是sine函数的两个图，一个在数据点之间用60个点，它比另一个只用6个点更光滑和更精确。 如曲线拟合一样，插值要作决策。根据所作的假设，有多种插值。而且，可以在一维以上空间中进行插值。即如果有反映两个变量函数的插值，z=f(x, y)，那么就可在x之间和在y之间，找出z的中间值进行插值。MATLAB在一维函数interp1和在二维函数interp2中，提供了许多的插值选择。其中的每个函数将在下面阐述。 为了说明一维插值，考虑下列问题，12小时内，一小时测量一次室外温度。数据存储在两个MATLAB变量中。 ? hours=1:12; % index for hour data was recorded ? temps=[5 8 9 15 25 29 31 30 22 25 27 24]; % recorded temperatures ? plot(hours, temps, hours, temps,' + ') % view temperatures ? title(' Temperature ') ? xlabel(' Hour '), ylabel(' Degrees Celsius ')  图11.4 在线性插值下室外温度曲线 正如图11.4看到的，MATLAB画出了数据点线性插值的直线。为了计算在任意给定时间的温度，人们可试着对可视的图作解释。另外一种方法，可用函数interp1。 ? t=interp1(hours, temps, 9.3) % estimate temperature at hour=9.3 t = 22.9000 ? t=interp1(hours, temps, 4.7) % estimate temperature at hour=4.7 t = 22 ? t=interp1(hours, temps, [3.2 6.5 7.1 11.7]) % find temp at many points! t = 10.2000 30.0000 30.9000 24.9000 interp1的缺省用法是由interp1(x, y, xo)来描述，这里x是独立变量(横坐标)，y是应变量(纵坐标)，xo是进行插值的一个数值数组。另外，该缺省的使用假定为线性插值。 若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的头两个导数与该数据点相一致。这种类型的插值被称为3次样条或简称为样条。函数interp1也能执行3次样条插值。 ? t=interp1(hours, temps, 9.3, ' spline ') % estimate temperature at hour=9.3 t = 21.8577 ? t=interp1(hours, temps, 4.7, ' spline ') % estimate temperature at hour=4.7 t = 22.3143 ? t=interp1(hours, temps, [3.2 6.5 7.1 11.7], ' spline ') t = 9.6734 30.0427 31.1755 25.3820 注意，样条插值得到的结果，与上面所示的线性插值的结果不同。因为插值是一个估计或猜测的过程，其意义在于，应用不同的估计规则导致不同的结果。 一个最常用的样条插值是对数据平滑。也就是，给定一组数据，使用样条插值在更细的间隔求值。例如， ? h=1: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次多项式以描述给定数据之间的特征。关于样条的更详细信息可见下一章。  图11.5 在不同插值下室外温度曲线 在讨论二维插值之前，了解interp1所强制的二个强约束是很重要的。首先，人们不能要求有独立变量范围以外的结果，例如，interp1(hours, temps, 13.5)导致一个错误，因为hours在1到12之间变化。其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减小的。在我们的例子里，hours是单调的。然而，如果我们已经定义独立变量为一天的实际时间， ? time_of_day=[7:12 1:6] % start at 7AM,end at 6PM time_of_day = 7 8 9 10 11 12 1 2 3 4 5 6 则独立变量将不是单调的，因为time_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours，将会返回一个错误。同样的理由，人们不能对temps插值来找出产生某温度的时间(小时)，因为temps不是单调的。 11.3 二维插值 二维插值是基于与一维插值同样的基本思想。然而，正如名字所隐含的，二维插值是对两变量的函数z=f(x, y) 进行插值。为了说明这个附加的维数，考虑一个问题。设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅。 采集了下列的数据： ? width=1: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=[82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86] % temperature data temps = 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 如同在标引点上测量一样，矩阵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 ') ; % 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 '), grid  图11.6 在深度d=2处的平板温度  图11.7 平板温度 然后在两个方向上插值，以平滑数据。 ? di=1:0.2:3; % choose higher resolution for depth ? wi=1: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的基本形式是interp2(x, y, z, xi, yi, method)。这里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'或'nearest'。在这种情况下，cubic不意味着3次样条，而是使用3次多项式的另一种算法。linear方法是线性插值，仅用作连接图上数据点。nearest方法只选择最接近各估计点的粗略数据点。在所有的情况下，假定独立变量x和y是线性间隔和单调的。关于这些方法的更多的信息，可请求在线帮助，例如，? 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.1 1.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 column oriented table from data ? xi=mminterp(table, 2, 1.1) xi = 0.5281 1.1000 0.9580 1.1000 1.5825 1.1000 1.8847 1.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. '), end if 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 > VAL below=tab(: , col) < val; % True where < VAL equal=tab(: , col) = = val; % True where = VAL if all(above = = 0) | all(below = = 0), % handle simplest case y=tab(find(equal), : ); return end pslope=find(below(1:rt-1)&above(2:rt)); % indices where 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]); % put indices above in order ie=find(equal); % indices where equal to val [tmp,ix]=sort( [ib, ie] ); % find where equals fit in result ieq=ix > length(ib); % True where equals values fit ry=length(tmp); % # of rows in result y y=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 columns y(~ieq, : )=alpha.*tab(ia, : )+(1-alpha).*tab(ib, : ); % interpolated values y(ieq, : )=tab(ie, : ); % equal values 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 ? table=[x; y; z].' ; ? t=mminterp(table, 2, 1.1) % same interpolation as earlier t = 0.5281 1.1000 0.9930 0.9580 1.1000 0.1314 1.5825 1.1000 -0.9639 1.8847 1.1000 -0.3533 ? t=mminterp(table, 3, -.5) % second third column now t = 1.1669 0.7316 -0.5000 1.8329 1.1377 -0.5000 3.1671 0.9639 -0.5000 3.8331 1.0187 -0.5000 这些最后的结果估计了x和y在z= -.5处的值。 尽管逐条地对函数mminterp解释如何工作是很有帮助的，但这样做要求有更多的篇幅和时间。解释mminterp如何工作最容易的方法是创建一个小表格，然后，在重要的语句末尾删除分号以后，调用函数。这样，中间值将帮助用户理解函数是如何找到与所需值相符的数据值以及如何执行插值。 前面已阐述了interp1的用法。当用于线性插值时，只要所要求的插值点的个数少，interp1工作很好。在要求许多插值点情况下，由于所用的算法，interp1工作较慢。为了克服这个问题，精通MATLAB工具箱包括了函数mmtable，它的帮助文本是： ?help table MMTABLE 1-D Monotonic Table Search by Linear Interpolation. YI=MMTABLE(TAB,COL,VALS) linearly interpolates the table TAB searching 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 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 usage of TABLE1(TAB,X0). YI=MMTABLE(X,Y,XI) interpolates the vector X to find YI associated with 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中搜索值为VALS TAB(:,COL)必须是单调的，但不必等价地生成空间。 YI与VALS有同样的行和与TAB有同样的列。 当VALS超出TAB(:,COL)的范围，返回NaNs. YI=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比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维最近邻插值