octave_packages/m/general/interp1.m

   1 ## Copyright (C) 2000-2012 Paul Kienzle
   2 ## Copyright (C) 2009 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {@var{yi} =} interp1 (@var{x}, @var{y}, @var{xi})
  22 ## @deftypefnx {Function File} {@var{yi} =} interp1 (@var{y}, @var{xi})
  23 ## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{method})
  24 ## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{extrap})
  25 ## @deftypefnx {Function File} {@var{pp} =} interp1 (@dots{}, 'pp')
  26 ##
  27 ## One-dimensional interpolation.  Interpolate @var{y}, defined at the
  28 ## points @var{x}, at the points @var{xi}.  The sample points @var{x}
  29 ## must be monotonic.  If not specified, @var{x} is taken to be the
  30 ## indices of @var{y}.  If @var{y} is an array, treat the columns
  31 ## of @var{y} separately.
  32 ##
  33 ## Method is one of:
  34 ##
  35 ## @table @asis
  36 ## @item 'nearest'
  37 ## Return the nearest neighbor.
  38 ##
  39 ## @item 'linear'
  40 ## Linear interpolation from nearest neighbors
  41 ##
  42 ## @item 'pchip'
  43 ## Piecewise cubic Hermite interpolating polynomial
  44 ##
  45 ## @item 'cubic'
  46 ## Cubic interpolation (same as @code{pchip})
  47 ##
  48 ## @item 'spline'
  49 ## Cubic spline interpolation---smooth first and second derivatives
  50 ## throughout the curve
  51 ## @end table
  52 ##
  53 ## Appending '*' to the start of the above method forces @code{interp1}
  54 ## to assume that @var{x} is uniformly spaced, and only @code{@var{x}
  55 ## (1)} and @code{@var{x} (2)} are referenced.  This is usually faster,
  56 ## and is never slower.  The default method is 'linear'.
  57 ##
  58 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond
  59 ## the endpoints.  If @var{extrap} is a number, replace values beyond the
  60 ## endpoints with that number.  If @var{extrap} is missing, assume NA.
  61 ##
  62 ## If the string argument 'pp' is specified, then @var{xi} should not be
  63 ## supplied and @code{interp1} returns the piecewise polynomial that
  64 ## can later be used with @code{ppval} to evaluate the interpolation.
  65 ## There is an equivalence, such that @code{ppval (interp1 (@var{x},
  66 ## @var{y}, @var{method}, 'pp'), @var{xi}) == interp1 (@var{x}, @var{y},
  67 ## @var{xi}, @var{method}, 'extrap')}.
  68 ##
  69 ## Duplicate points in @var{x} specify a discontinuous interpolant.  There
  70 ## should be at most 2 consecutive points with the same value.
  71 ## The discontinuous interpolant is right-continuous if @var{x} is increasing,
  72 ## left-continuous if it is decreasing.
  73 ## Discontinuities are (currently) only allowed for "nearest" and "linear"
  74 ## methods; in all other cases, @var{x} must be strictly monotonic.
  75 ##
  76 ## An example of the use of @code{interp1} is
  77 ##
  78 ## @example
  79 ## @group
  80 ## xf = [0:0.05:10];
  81 ## yf = sin (2*pi*xf/5);
  82 ## xp = [0:10];
  83 ## yp = sin (2*pi*xp/5);
  84 ## lin = interp1 (xp, yp, xf);
  85 ## spl = interp1 (xp, yp, xf, "spline");
  86 ## cub = interp1 (xp, yp, xf, "cubic");
  87 ## near = interp1 (xp, yp, xf, "nearest");
  88 ## plot (xf, yf, "r", xf, lin, "g", xf, spl, "b",
  89 ##       xf, cub, "c", xf, near, "m", xp, yp, "r*");
  90 ## legend ("original", "linear", "spline", "cubic", "nearest");
  91 ## @end group
  92 ## @end example
  93 ##
  94 ## @seealso{interpft}
  95 ## @end deftypefn
  96
  97 ## Author: Paul Kienzle
  98 ## Date: 2000-03-25
  99 ##    added 'nearest' as suggested by Kai Habel
 100 ## 2000-07-17 Paul Kienzle
 101 ##    added '*' methods and matrix y
 102 ##    check for proper table lengths
 103 ## 2002-01-23 Paul Kienzle
 104 ##    fixed extrapolation
 105
 106 function yi = interp1 (x, y, varargin)
 107
 108   if (nargin < 2 || nargin > 6)
 109     print_usage ();
 110   endif
 111
 112   method = "linear";
 113   extrap = NA;
 114   xi = [];
 115   ispp = false;
 116   firstnumeric = true;
 117
 118   if (nargin > 2)
 119     for i = 1:length (varargin)
 120       arg = varargin{i};
 121       if (ischar (arg))
 122         arg = tolower (arg);
 123         if (strcmp ("extrap", arg))
 124           extrap = "extrap";
 125         elseif (strcmp ("pp", arg))
 126           ispp = true;
 127         else
 128           method = arg;
 129         endif
 130       else
 131         if (firstnumeric)
 132           xi = arg;
 133           firstnumeric = false;
 134         else
 135           extrap = arg;
 136         endif
 137       endif
 138     endfor
 139   endif
 140
 141   if (isempty (xi) && firstnumeric && ! ispp)
 142     xi = y;
 143     y = x;
 144     x = 1:numel(y);
 145   endif
 146
 147   ## reshape matrices for convenience
 148   x = x(:);
 149   nx = rows (x);
 150   szx = size (xi);
 151   if (isvector (y))
 152     y = y(:);
 153   endif
 154
 155   szy = size (y);
 156   y = y(:,:);
 157   [ny, nc] = size (y);
 158   xi = xi(:);
 159
 160   ## determine sizes
 161   if (nx < 2 || ny < 2)
 162     error ("interp1: table too short");
 163   endif
 164
 165   ## check whether x is sorted; sort if not.
 166   if (! issorted (x, "either"))
 167     [x, p] = sort (x);
 168     y = y(p,:);
 169   endif
 170
 171   starmethod = method(1) == "*";
 172
 173   if (starmethod)
 174     dx = x(2) - x(1);
 175   else
 176     jumps = x(1:nx-1) == x(2:nx);
 177     have_jumps = any (jumps);
 178     if (have_jumps)
 179       if (any (strcmp (method, {"nearest", "linear"})))
 180         if (any (jumps(1:nx-2) & jumps(2:nx-1)))
 181           warning ("interp1: extra points in discontinuities");
 182         endif
 183       else
 184         error ("interp1: discontinuities not supported for method %s", method);
 185       endif
 186     endif
 187   endif
 188
 189   ## Proceed with interpolating by all methods.
 190
 191   switch (method)
 192   case "nearest"
 193     pp = mkpp ([x(1); (x(1:nx-1)+x(2:nx))/2; x(nx)], shiftdim (y, 1), szy(2:end));
 194     pp.orient = "first";
 195
 196     if (ispp)
 197       yi = pp;
 198     else
 199       yi = ppval (pp, reshape (xi, szx));
 200     endif
 201   case "*nearest"
 202     pp = mkpp ([x(1), x(1)+[0.5:(nx-1)]*dx, x(nx)], shiftdim (y, 1), szy(2:end));
 203     pp.orient = "first";
 204     if (ispp)
 205       yi = pp;
 206     else
 207       yi = ppval(pp, reshape (xi, szx));
 208     endif
 209   case "linear"
 210     dy = diff (y);
 211     dx = diff (x);
 212     dx = repmat (dx, [1 size(dy)(2:end)]);
 213     coefs = [(dy./dx).'(:), y(1:nx-1, :).'(:)];
 214     xx = x;
 215
 216     if (have_jumps)
 217       ## Omit zero-size intervals.
 218       coefs(jumps, :) = [];
 219       xx(jumps) = [];
 220     endif
 221
 222     pp = mkpp (xx, coefs, szy(2:end));
 223     pp.orient = "first";
 224
 225     if (ispp)
 226       yi = pp;
 227     else
 228       yi = ppval(pp, reshape (xi, szx));
 229     endif
 230
 231   case "*linear"
 232     dy = diff (y);
 233     coefs = [(dy/dx).'(:), y(1:nx-1, :).'(:)];
 234     pp = mkpp (x, coefs, szy(2:end));
 235     pp.orient = "first";
 236
 237     if (ispp)
 238       yi = pp;
 239     else
 240       yi = ppval(pp, reshape (xi, szx));
 241     endif
 242
 243   case {"pchip", "*pchip", "cubic", "*cubic"}
 244     if (nx == 2 || starmethod)
 245       x = linspace (x(1), x(nx), ny);
 246     endif
 247
 248     if (ispp)
 249       y = shiftdim (reshape (y, szy), 1);
 250       yi = pchip (x, y);
 251     else
 252       y = shiftdim (y, 1);
 253       yi = pchip (x, y, reshape (xi, szx));
 254     endif
 255   case {"spline", "*spline"}
 256     if (nx == 2 || starmethod)
 257       x = linspace(x(1), x(nx), ny);
 258     endif
 259
 260     if (ispp)
 261       y = shiftdim (reshape (y, szy), 1);
 262       yi = spline (x, y);
 263     else
 264       y = shiftdim (y, 1);
 265       yi = spline (x, y, reshape (xi, szx));
 266     endif
 267   otherwise
 268     error ("interp1: invalid method '%s'", method);
 269   endswitch
 270
 271   if (! ispp)
 272     if (! ischar (extrap))
 273       ## determine which values are out of range and set them to extrap,
 274       ## unless extrap == "extrap".
 275       minx = min (x(1), x(nx));
 276       maxx = max (x(1), x(nx));
 277
 278       outliers = xi < minx | ! (xi <= maxx); # this catches even NaNs
 279       if (size_equal (outliers, yi))
 280         yi(outliers) = extrap;
 281         yi = reshape (yi, szx);
 282       elseif (!isvector (yi))
 283         if (strcmp (method, "pchip") || strcmp (method, "*pchip")
 284           ||strcmp (method, "cubic") || strcmp (method, "*cubic")
 285           ||strcmp (method, "spline") || strcmp (method, "*spline"))
 286           yi(:, outliers) = extrap;
 287           yi = shiftdim(yi, 1);
 288         else
 289           yi(outliers, :) = extrap;
 290         endif
 291       else
 292         yi(outliers.') = extrap;
 293       endif
 294     endif
 295   else
 296     yi.orient = "first";
 297   endif
 298
 299 endfunction
 300
 301 %!demo
 302 %! xf=0:0.05:10; yf = sin(2*pi*xf/5);
 303 %! xp=0:10;      yp = sin(2*pi*xp/5);
 304 %! lin=interp1(xp,yp,xf,"linear");
 305 %! spl=interp1(xp,yp,xf,"spline");
 306 %! cub=interp1(xp,yp,xf,"pchip");
 307 %! near=interp1(xp,yp,xf,"nearest");
 308 %! plot(xf,yf,"r",xf,near,"g",xf,lin,"b",xf,cub,"c",xf,spl,"m",xp,yp,"r*");
 309 %! legend ("original","nearest","linear","pchip","spline")
 310 %! %--------------------------------------------------------
 311 %! % confirm that interpolated function matches the original
 312
 313 %!demo
 314 %! xf=0:0.05:10; yf = sin(2*pi*xf/5);
 315 %! xp=0:10;      yp = sin(2*pi*xp/5);
 316 %! lin=interp1(xp,yp,xf,"*linear");
 317 %! spl=interp1(xp,yp,xf,"*spline");
 318 %! cub=interp1(xp,yp,xf,"*cubic");
 319 %! near=interp1(xp,yp,xf,"*nearest");
 320 %! plot(xf,yf,"r",xf,near,"g",xf,lin,"b",xf,cub,"c",xf,spl,"m",xp,yp,"r*");
 321 %! legend ("*original","*nearest","*linear","*cubic","*spline")
 322 %! %--------------------------------------------------------
 323 %! % confirm that interpolated function matches the original
 324
 325 %!demo
 326 %! t = 0 : 0.3 : pi; dt = t(2)-t(1);
 327 %! n = length (t); k = 100; dti = dt*n/k;
 328 %! ti = t(1) + [0 : k-1]*dti;
 329 %! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
 330 %! ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
 331 %! ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
 332 %! ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
 333 %! plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
 334 %!       ti(2:end-1),ddyp,'c^');
 335 %! legend('cubic','spline','pchip');
 336 %! title("Second derivative of interpolated 'sin (4*t + 0.3) .* cos (3*t - 0.1)'");
 337
 338 %!demo
 339 %! xf=0:0.05:10; yf = sin(2*pi*xf/5) - (xf >= 5);
 340 %! xp=[0:.5:4.5,4.99,5:.5:10];      yp = sin(2*pi*xp/5) - (xp >= 5);
 341 %! lin=interp1(xp,yp,xf,"linear");
 342 %! near=interp1(xp,yp,xf,"nearest");
 343 %! plot(xf,yf,"r",xf,near,"g",xf,lin,"b",xp,yp,"r*");
 344 %! legend ("original","nearest","linear")
 345 %! %--------------------------------------------------------
 346 %! % confirm that interpolated function matches the original
 347
 348 ##FIXME: add test for n-d arguments here
 349
 350 ## For each type of interpolated test, confirm that the interpolated
 351 ## value at the knots match the values at the knots.  Points away
 352 ## from the knots are requested, but only 'nearest' and 'linear'
 353 ## confirm they are the correct values.
 354
 355 %!shared xp, yp, xi, style
 356 %! xp=0:2:10;      yp = sin(2*pi*xp/5);
 357 %! xi = [-1, 0, 2.2, 4, 6.6, 10, 11];
 358
 359
 360 ## The following BLOCK/ENDBLOCK section is repeated for each style
 361 ##    nearest, linear, cubic, spline, pchip
 362 ## The test for ppval of cubic has looser tolerance, but otherwise
 363 ## the tests are identical.
 364 ## Note that the block checks style and *style; if you add more tests
 365 ## before to add them to both sections of each block.  One test,
 366 ## style vs. *style, occurs only in the first section.
 367 ## There is an ENDBLOCKTEST after the final block
 368 %!test style = "nearest";
 369 ## BLOCK
 370 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 371 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 372 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 373 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 374 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 375 %!assert (isempty(interp1(xp',yp',[],style)));
 376 %!assert (isempty(interp1(xp,yp,[],style)));
 377 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 378 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 379 %!assert (interp1(xp,yp,xi,style),...
 380 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 381 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 382 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 383 %!error interp1(1,1,1, style);
 384 %!assert (interp1(xp,[yp',yp'],xi,style),
 385 %!        interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
 386 %!test style=['*',style];
 387 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 388 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 389 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 390 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 391 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 392 %!assert (isempty(interp1(xp',yp',[],style)));
 393 %!assert (isempty(interp1(xp,yp,[],style)));
 394 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 395 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 396 %!assert (interp1(xp,yp,xi,style),...
 397 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 398 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 399 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 400 %!error interp1(1,1,1, style);
 401 ## ENDBLOCK
 402 %!test style='linear';
 403 ## BLOCK
 404 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 405 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 406 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 407 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 408 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 409 %!assert (isempty(interp1(xp',yp',[],style)));
 410 %!assert (isempty(interp1(xp,yp,[],style)));
 411 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 412 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 413 %!assert (interp1(xp,yp,xi,style),...
 414 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 415 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 416 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 417 %!error interp1(1,1,1, style);
 418 %!assert (interp1(xp,[yp',yp'],xi,style),
 419 %!        interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
 420 %!test style=['*',style];
 421 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 422 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 423 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 424 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 425 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 426 %!assert (isempty(interp1(xp',yp',[],style)));
 427 %!assert (isempty(interp1(xp,yp,[],style)));
 428 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 429 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 430 %!assert (interp1(xp,yp,xi,style),...
 431 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 432 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 433 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 434 %!error interp1(1,1,1, style);
 435 ## ENDBLOCK
 436 %!test style='cubic';
 437 ## BLOCK
 438 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 439 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 440 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 441 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 442 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 443 %!assert (isempty(interp1(xp',yp',[],style)));
 444 %!assert (isempty(interp1(xp,yp,[],style)));
 445 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 446 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 447 %!assert (interp1(xp,yp,xi,style),...
 448 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 449 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 450 %!        interp1(xp,yp,xi,style,"extrap"),100*eps);
 451 %!error interp1(1,1,1, style);
 452 %!assert (interp1(xp,[yp',yp'],xi,style),
 453 %!        interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
 454 %!test style=['*',style];
 455 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 456 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 457 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 458 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 459 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 460 %!assert (isempty(interp1(xp',yp',[],style)));
 461 %!assert (isempty(interp1(xp,yp,[],style)));
 462 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 463 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 464 %!assert (interp1(xp,yp,xi,style),...
 465 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 466 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 467 %!        interp1(xp,yp,xi,style,"extrap"),100*eps);
 468 %!error interp1(1,1,1, style);
 469 ## ENDBLOCK
 470 %!test style='pchip';
 471 ## BLOCK
 472 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 473 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 474 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 475 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 476 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 477 %!assert (isempty(interp1(xp',yp',[],style)));
 478 %!assert (isempty(interp1(xp,yp,[],style)));
 479 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 480 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 481 %!assert (interp1(xp,yp,xi,style),...
 482 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 483 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 484 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 485 %!error interp1(1,1,1, style);
 486 %!assert (interp1(xp,[yp',yp'],xi,style),
 487 %!        interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
 488 %!test style=['*',style];
 489 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 490 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 491 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 492 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 493 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 494 %!assert (isempty(interp1(xp',yp',[],style)));
 495 %!assert (isempty(interp1(xp,yp,[],style)));
 496 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 497 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 498 %!assert (interp1(xp,yp,xi,style),...
 499 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 500 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 501 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 502 %!error interp1(1,1,1, style);
 503 ## ENDBLOCK
 504 %!test style='spline';
 505 ## BLOCK
 506 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 507 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 508 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 509 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 510 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 511 %!assert (isempty(interp1(xp',yp',[],style)));
 512 %!assert (isempty(interp1(xp,yp,[],style)));
 513 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 514 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 515 %!assert (interp1(xp,yp,xi,style),...
 516 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 517 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 518 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 519 %!error interp1(1,1,1, style);
 520 %!assert (interp1(xp,[yp',yp'],xi,style),
 521 %!        interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
 522 %!test style=['*',style];
 523 %!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 524 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 525 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 526 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
 527 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
 528 %!assert (isempty(interp1(xp',yp',[],style)));
 529 %!assert (isempty(interp1(xp,yp,[],style)));
 530 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 531 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
 532 %!assert (interp1(xp,yp,xi,style),...
 533 %!        interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
 534 %!assert (ppval(interp1(xp,yp,style,"pp"),xi),
 535 %!        interp1(xp,yp,xi,style,"extrap"),10*eps);
 536 %!error interp1(1,1,1, style);
 537 ## ENDBLOCK
 538 ## ENDBLOCKTEST
 539
 540 %!# test linear extrapolation
 541 %!assert (interp1([1:5],[3:2:11],[0,6],"linear","extrap"), [1, 13], eps);
 542 %!assert (interp1(xp, yp, [-1, max(xp)+1],"linear",5), [5, 5]);
 543
 544 %!error interp1
 545 %!error interp1(1:2,1:2,1,"bogus")
 546
 547 %!assert (interp1(1:2,1:2,1.4,"nearest"),1);
 548 %!error interp1(1,1,1, "linear");
 549 %!assert (interp1(1:2,1:2,1.4,"linear"),1.4);
 550 %!assert (interp1(1:4,1:4,1.4,"cubic"),1.4);
 551 %!assert (interp1(1:2,1:2,1.1, "spline"), 1.1);
 552 %!assert (interp1(1:3,1:3,1.4,"spline"),1.4);
 553
 554 %!error interp1(1,1,1, "*nearest");
 555 %!assert (interp1(1:2:4,1:2:4,1.4,"*nearest"),1);
 556 %!error interp1(1,1,1, "*linear");
 557 %!assert (interp1(1:2:4,1:2:4,[0,1,1.4,3,4],"*linear"),[NA,1,1.4,3,NA]);
 558 %!assert (interp1(1:2:8,1:2:8,1.4,"*cubic"),1.4);
 559 %!assert (interp1(1:2,1:2,1.3, "*spline"), 1.3);
 560 %!assert (interp1(1:2:6,1:2:6,1.4,"*spline"),1.4);
 561
 562 %!assert (interp1([3,2,1],[3,2,2],2.5),2.5)
 563
 564 %!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap"), [-2,0.5,4,3,1.5])
 565 %!assert (interp1 ([4,4,3,2,0],[0,1,4,2,1],[1.5,4,4.5], "linear"), [1.75,1,NA])
 566 %!assert (interp1 (0:4, 2.5), 1.5)