octave_packages/m/general/interp1q.m

   1 ## Copyright (C) 2008-2012 David Bateman
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{yi} =} interp1q (@var{x}, @var{y}, @var{xi})
  21 ## One-dimensional linear interpolation without error checking.
  22 ## Interpolates @var{y}, defined at the points @var{x}, at the points
  23 ## @var{xi}.  The sample points @var{x} must be a strictly monotonically
  24 ## increasing column vector.  If @var{y} is an array, treat the columns
  25 ## of @var{y} separately.  If @var{y} is a vector, it must be a column
  26 ## vector of the same length as @var{x}.
  27 ##
  28 ## Values of @var{xi} beyond the endpoints of the interpolation result
  29 ## in NA being returned.
  30 ##
  31 ## Note that the error checking is only a significant portion of the
  32 ## execution time of this @code{interp1} if the size of the input arguments
  33 ## is relatively small.  Therefore, the benefit of using @code{interp1q}
  34 ## is relatively small.
  35 ## @seealso{interp1}
  36 ## @end deftypefn
  37
  38 function yi = interp1q (x, y, xi)
  39   x = x(:);
  40   nx = size (x, 1);
  41   szy = size (y);
  42   y = y(:,:);
  43   [ny, nc] = size (y);
  44   szx = size (xi);
  45   xi = xi (:);
  46   dy = diff (y);
  47   dx = diff (x);
  48   idx = lookup (x, xi, "lr");
  49   s = (xi - x (idx)) ./ dx (idx);
  50   yi = bsxfun (@times, s, dy(idx,:)) + y(idx,:);
  51   range = xi < x(1) | !(xi <= x(nx));
  52   yi(range,:) = NA;
  53   if (length (szx) == 2 && any (szx == 1))
  54     yi = reshape (yi, [max(szx), szy(2:end)]);
  55   else
  56     yi = reshape (yi, [szx, szy(2:end)]);
  57   endif
  58 endfunction
  59
  60 %!shared xp, yp, xi, yi
  61 %! xp=[0:2:10].';      yp = sin(2*pi*xp/5);
  62 %! xi = [-1; 0; 2.2; 4; 6.6; 10; 11];
  63 %! yi = interp1 (xp,yp,xi);
  64 %!assert (interp1q(xp, yp, [min(xp)-1; max(xp)+1]), [NA; NA]);
  65 %!assert (interp1q(xp,yp,xp), yp, 100*eps);
  66 %!assert (isempty(interp1q(xp,yp,[])));
  67 %!assert (interp1q(xp,yp,xi), yi);
  68 %!assert (interp1q(xp,[yp,yp],xi), [yi, yi]);
  69 %!assert (interp1q(xp,yp,[xi,xi]), [yi, yi]);
  70 %!assert (interp1q(xp,[yp,yp],[xi,xi]), cat (3, [yi, yi], [yi, yi]));