这里使用了线性插值,函数mminterp估计了y=1.1处的四个点。由于函数mminterp的一般性质,要插值的数据是由面向列矩阵给出,在上面的例子中称作为表(table)。第二个输入参量是被搜索矩阵table的列,第三个参量是要找的值。 这个精通MATLAB工具箱函数的主体由下面给出:
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 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 - above=tab(: , col) > val; % True where > VAL below=tab(: , col) < val; % True where < VAL equal=tab(: , col) = = val; % True where = VAL [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 % Copyright (c) 1996 by Prentice-Hall,Inc. 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=[].
正如所见的,mminterp利用了find和sort函数、逻辑数组和数组操作技术。没有For循环和While循环。不论用其中哪一种技术来实现将使运行变慢,尤其对大的表。注意mminterp与含有大于或等于2的任意数列的表一起工作,如同函数interp1一样。而且,在这种情况下,插值变量可以是任意的列。例如,
? t=mminterp(table, 3, -.5) % second third column now t =
1.1669 0.7316 -0.5000
? 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
? table=[x; y; z].' ;
? z=sin(pi*x); % add more data to table y(ieq, : )=tab(ie, : ); % equal values
y( : , col)=val*ones(ry, 1); % remove roundoff error 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=zeros(ry, ct); % poke data into a zero matrix
[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
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) interp1(x, y, xo, ' spline ') interp1(x, y, xo, ' cubic ') interp2(x, y, Z, xi, yi) interp2(x, y, Z, xi, yi, ' cubic ') interp2(x, y, Z, xi, yi, ' nearest ')
1维线性插值 1维3次样条插值 1维3次插值 2维线性插值 2维3次插值 2维最近邻插值