changeset 7785:9a656cd34f34

Add the interp1q function
author David Bateman <dbateman@free.fr>
date Wed, 14 May 2008 22:51:19 +0200
parents d3a7882fa0b3
children 37ff0c21c17d
files scripts/ChangeLog scripts/general/Makefile.in scripts/general/interp1q.m
diffstat 3 files changed, 76 insertions(+), 1 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Tue May 20 16:24:50 2008 -0400
+++ b/scripts/ChangeLog	Wed May 14 22:51:19 2008 +0200
@@ -1,3 +1,8 @@
+2008-05-20  David Bateman  <dbateman@free.fr>
+
+	* general/interp1q.m: New function.
+	* general/Makefile.in (SOURCES): Add it here.
+
 2008-05-20  Bill Denney  <bill@denney.ws>
 
 	* time/datenum.m: Allow mixed scalar and vector/matrix input.
--- a/scripts/general/Makefile.in	Tue May 20 16:24:50 2008 -0400
+++ b/scripts/general/Makefile.in	Wed May 14 22:51:19 2008 +0200
@@ -38,7 +38,7 @@
   cart2sph.m cell2mat.m celldisp.m circshift.m common_size.m \
   cplxpair.m cumtrapz.m dblquad.m deal.m del2.m diff.m flipdim.m fliplr.m \
   flipud.m genvarname.m gradient.m idivide.m ind2sub.m int2str.m \
-  interp1.m interp2.m interp3.m interpn.m interpft.m \
+  interp1.m interp1q.m interp2.m interp3.m interpn.m interpft.m \
   is_duplicate_entry.m isa.m isdefinite.m isdir.m isequal.m \
   isequalwithequalnans.m isscalar.m issquare.m issymmetric.m \
   isvector.m logical.m logspace.m mod.m nargchk.m \
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/interp1q.m	Wed May 14 22:51:19 2008 +0200
@@ -0,0 +1,70 @@
+## Copyright (C) 2008 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{yi} =} interp1q (@var{x}, @var{y}, @var{xi})
+## One-dimensional linear interpolation without error checking.
+## Interpolates @var{y}, defined at the points @var{x}, at the points
+## @var{xi}. The sample points @var{x} must be a strictly monotonically
+## increasing column vector. If @var{y} is an array, treat the columns
+## of @var{y} separately. If @var{y} is a vector, it must be a column
+## vector of the same length as @var{x}.
+##
+## Values of @var{xi} beyond the endpoints of the interpolation result
+## in NA being returned.
+##
+## Note that the error checking is only a significant portion of the
+## execution time of this @code{interp1} if the size of the input arguments
+## is relatively small. Therefore, the benefit of using @code{interp1q}
+## is relatively small.
+## @seealso{interp1}
+## @end deftypefn
+
+function yi = interp1q (x, y, xi)
+  x = x(:);
+  nx = size (x, 1);
+  szy = size (y);
+  ny = szy (1);
+  nc = prod (szy (2 : end));
+  y = reshape (y, ny, nc);
+  szx = size (xi);
+  xi = xi (:);
+  range = find (xi >= x (1) & xi <= x (nx));
+  yi = NA (size(xi, 1), size (y, 2));
+  xi = xi (range);
+  dy = y (2 : ny, :) - y (1 : ny - 1, :);
+  dx = x (2 : nx) - x (1 : nx - 1);
+  idx = lookup (x, xi, "lr");
+  s = (xi - x (idx)) ./ dx (idx);
+  yi (range, :) = s (:, ones (1, nc)) .* dy (idx, :) + y (idx, :);
+  if (length (szx) == 2 && any (szx == 1))
+    yi = reshape (yi, [max (szx), szy (2 : end)]);
+  else
+    yi = reshape (yi, [szx, szy(2:end)]);
+  endif
+endfunction
+
+%!shared xp, yp, xi
+%! xp=[0:2:10].';      yp = sin(2*pi*xp/5);
+%! xi = [-1; 0; 2.2; 4; 6.6; 10; 11];
+%!assert (interp1(xp, yp, [min(xp)-1; max(xp)+1]), [NA; NA]);
+%!assert (interp1(xp,yp,xp), yp, 100*eps);
+%!assert (isempty(interp1(xp,yp,[])));
+%!assert (interp1(xp,[yp,yp],xi), [interp1(xp,yp,xi),interp1(xp,yp,xi)]);
+%!assert (interp1(xp,yp,[xi,xi]), [interp1(xp,yp,xi),interp1(xp,yp,xi)]);
+%!assert (interp1(xp,[yp,yp],[xi,xi]), [interp1(xp,yp,xi),interp1(xp,yp,xi)]);