www.pudn.com > 精通Matlab综合辅导与指南-源程序.zip > mmtable.m

function y=mmtable(a,b,c)
%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 as 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 matches the usage of INTERP1(X,Y,XI)
%
% This routine is 10X faster than TABLE1 which is called by INTERP1.
% For linearly spaced data, ILINEAR is 10X faster yet.

% D.C. Hanselman, University of Maine, Orono ME,  04469
% 8/21/95
% Copyright (c) 1996 by Prentice-Hall, Inc. 


if nargin==2          % MMTABLE(TAB,VALS) syntax
	tab=a;vals=b(:);col=1;cflag=1;
elseif length(b)>1    % MMTABLE(X,Y,XI) syntax
	[rb,cb]=size(b);
	a=a(:);
	la=length(a);
	if rb==la,tab=[a b];
	elseif cb==la,tab=[a b.'];
	end
	col=1;vals=c(:);cflag=1;
else                  % MMTABLE(TAB,COL,VALS) syntax
	tab=a;col=b(1);vals=c(:);cflag=0;
end

[rt,ct]=size(tab);
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
ict=1:ct;      % index for table columns
nv=length(vals);

dx=diff(tab);
dx=[dx;dx(rt-1,:)];  % differences within each column.
id=sign(dx(1,col));  % sign of differences in col

if any(sign(dx(:,col))~=id)
	error(sprintf('Column %.0f of table is not monotonic.',col))
end

inan=valsmax(tab(1,col),tab(rt,col));
xo=vals(~inan);  % pick values within range of tab(:,col)
nxo=length(xo);
if nxo==0  % all values are outside range of tab(:,col)!
	y=NaN*ones(nv,ct);
	y(:,col)=vals;
else  % interpolate using algorithm posted by Hans Olsson of Lund University

	
	[svals,sidx]=sort(id*[tab(2:rt-1,col);xo]);
	ttab=sidx<=rt-2;      % True at tab(:,col) values in svals
	cidx=1+cumsum(ttab);  % index of xo in svals
	sidx(ttab)=[];        % throw out indices of tab(:,col) values
	yidx=zeros(nxo,1);    % storage for xo index in tab(:,col)
	yidx(sidx-(rt-2))=cidx(~ttab);  % xo index in tab(:,col)
	ndx=(xo-tab(yidx,col))./dx(yidx,col); % normalized dx increment
	y=tab(yidx,ict)+ndx(:,ones(1,ct)).*dx(yidx,ict);
	if ~isempty(inan)  % must poke in NaNs where they belong
		yt=zeros(nv,ct);
		yt(~inan,:)=y;
		yt(inan,:)=NaN*ones(nv-nxo,ct);
		yt(:,col)=vals;
		y=yt;
	end
end
if cflag,y(:,col)=[];end