1########################################################################
2##
3## Copyright (C) 2000-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {@var{yi} =} interp1 (@var{x}, @var{y}, @var{xi})
28## @deftypefnx {} {@var{yi} =} interp1 (@var{y}, @var{xi})
29## @deftypefnx {} {@var{yi} =} interp1 (@dots{}, @var{method})
30## @deftypefnx {} {@var{yi} =} interp1 (@dots{}, @var{extrap})
31## @deftypefnx {} {@var{yi} =} interp1 (@dots{}, "left")
32## @deftypefnx {} {@var{yi} =} interp1 (@dots{}, "right")
33## @deftypefnx {} {@var{pp} =} interp1 (@dots{}, "pp")
34##
35## One-dimensional interpolation.
36##
37## Interpolate input data to determine the value of @var{yi} at the points
38## @var{xi}.  If not specified, @var{x} is taken to be the indices of @var{y}
39## (@code{1:length (@var{y})}).  If @var{y} is a matrix or an N-dimensional
40## array, the interpolation is performed on each column of @var{y}.
41##
42## The interpolation @var{method} is one of:
43##
44## @table @asis
45## @item @qcode{"nearest"}
46## Return the nearest neighbor.
47##
48## @item @qcode{"previous"}
49## Return the previous neighbor.
50##
51## @item @qcode{"next"}
52## Return the next neighbor.
53##
54## @item @qcode{"linear"} (default)
55## Linear interpolation from nearest neighbors.
56##
57## @item @qcode{"pchip"}
58## Piecewise cubic Hermite interpolating polynomial---shape-preserving
59## interpolation with smooth first derivative.
60##
61## @item @qcode{"cubic"}
62## Cubic interpolation (same as @qcode{"pchip"}).
63##
64## @item @qcode{"spline"}
65## Cubic spline interpolation---smooth first and second derivatives
66## throughout the curve.
67## @end table
68##
69## Adding '*' to the start of any method above forces @code{interp1}
70## to assume that @var{x} is uniformly spaced, and only @code{@var{x}(1)}
71## and @code{@var{x}(2)} are referenced.  This is usually faster,
72## and is never slower.  The default method is @qcode{"linear"}.
73##
74## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values
75## beyond the endpoints using the current @var{method}.  If @var{extrap} is a
76## number, then replace values beyond the endpoints with that number.  When
77## unspecified, @var{extrap} defaults to @code{NA}.
78##
79## If the string argument @qcode{"pp"} is specified, then @var{xi} should not
80## be supplied and @code{interp1} returns a piecewise polynomial object.  This
81## object can later be used with @code{ppval} to evaluate the interpolation.
82## There is an equivalence, such that @code{ppval (interp1 (@var{x},
83## @var{y}, @var{method}, @qcode{"pp"}), @var{xi}) == interp1 (@var{x},
84## @var{y}, @var{xi}, @var{method}, @qcode{"extrap"})}.
85##
86## Duplicate points in @var{x} specify a discontinuous interpolant.  There
87## may be at most 2 consecutive points with the same value.
88## If @var{x} is increasing, the default discontinuous interpolant is
89## right-continuous.  If @var{x} is decreasing, the default discontinuous
90## interpolant is left-continuous.
91## The continuity condition of the interpolant may be specified by using
92## the options @qcode{"left"} or @qcode{"right"} to select a left-continuous
93## or right-continuous interpolant, respectively.
94## Discontinuous interpolation is only allowed for @qcode{"nearest"} and
95## @qcode{"linear"} methods; in all other cases, the @var{x}-values must be
96## unique.
97##
98## An example of the use of @code{interp1} is
99##
100## @example
101## @group
102## xf = [0:0.05:10];
103## yf = sin (2*pi*xf/5);
104## xp = [0:10];
105## yp = sin (2*pi*xp/5);
106## lin = interp1 (xp, yp, xf);
107## near = interp1 (xp, yp, xf, "nearest");
108## pch = interp1 (xp, yp, xf, "pchip");
109## spl = interp1 (xp, yp, xf, "spline");
110## plot (xf,yf,"r", xf,near,"g", xf,lin,"b", xf,pch,"c", xf,spl,"m",
111##       xp,yp,"r*");
112## legend ("original", "nearest", "linear", "pchip", "spline");
113## @end group
114## @end example
115##
116## @seealso{pchip, spline, interpft, interp2, interp3, interpn}
117## @end deftypefn
118
119function yi = interp1 (x, y, varargin)
120
121  if (nargin < 2 || nargin > 6)
122    print_usage ();
123  endif
124
125  method = "linear";
126  extrap = NA;
127  xi = [];
128  ispp = false;
129  have_xi = false;
130  rightcontinuous = NaN;
131
132  if (nargin > 2)
133    for i_arg = 1:length (varargin)
134      arg = varargin{i_arg};
135      if (ischar (arg))
136        arg = tolower (arg);
137        switch (arg)
138          case "extrap"
139            extrap = "extrap";
140          case "pp"
141            ispp = true;
142          case {"right", "-right"}
143            rightcontinuous = true;
144          case {"left", "-left"}
145            rightcontinuous = false;
146          otherwise
147            method = arg;
148        endswitch
149      else
150        if (i_arg == 1)
151          xi = arg;
152          have_xi = true;
153        else
154          extrap = arg;
155        endif
156      endif
157    endfor
158  endif
159
160  if (! have_xi && ! ispp)
161    xi = y;
162    y = x;
163    if (isvector (y))
164      x = 1:numel (y);
165    else
166      x = 1:rows (y);
167    endif
168  endif
169
170  ## reshape matrices for convenience
171  x = x(:);
172  nx = rows (x);
173  szx = size (xi);
174  if (isvector (y))
175    y = y(:);
176  endif
177
178  szy = size (y);
179  y = y(:,:);
180  [ny, nc] = size (y);
181  xi = xi(:);
182
183  ## determine sizes
184  if (nx < 2 || ny < 2)
185    error ("interp1: minimum of 2 points required in each dimension");
186  endif
187
188  ## check whether x is sorted; sort if not.
189  if (! issorted (x, "either"))
190    [x, p] = sort (x);
191    y = y(p,:);
192  endif
193
194  if (any (strcmp (method, {"previous", "*previous", "next", "*next"})))
195    rightcontinuous = NaN; # needed for these methods to work
196  endif
197
198  if (isnan (rightcontinuous))
199    ## If not specified, set the continuity condition
200    if (x(end) < x(1))
201      rightcontinuous = false;
202    else
203      rightcontinuous = true;
204    endif
205  elseif ((rightcontinuous && (x(end) < x(1)))
206          || (! rightcontinuous && (x(end) > x(1))))
207    ## Switch between left-continuous and right-continuous
208    x = flipud (x);
209    y = flipud (y);
210  endif
211
212  ## Because of the way mkpp works, it's easiest to implement "next"
213  ## by running "previous" with vectors flipped.
214  if (strcmp (method, "next"))
215    x = flipud (x);
216    y = flipud (y);
217    method = "previous";
218  elseif (strcmp (method, "*next"))
219    x = flipud (x);
220    y = flipud (y);
221    method = "*previous";
222  endif
223
224  starmethod = method(1) == "*";
225
226  if (starmethod)
227    dx = x(2) - x(1);
228  else
229    jumps = x(1:end-1) == x(2:end);
230    have_jumps = any (jumps);
231    if (have_jumps)
232      if (strcmp (method, "linear") || strcmp (method, ("nearest")))
233        if (any (jumps(1:nx-2) & jumps(2:nx-1)))
234          warning ("interp1: multiple discontinuities at the same X value");
235        endif
236      else
237        error ("interp1: discontinuities not supported for METHOD '%s'",
238                                                                   method);
239      endif
240    endif
241  endif
242
243  ## Proceed with interpolating by all methods.
244  switch (method)
245
246    case "nearest"
247      pp = mkpp ([x(1); (x(1:nx-1)+x(2:nx))/2; x(nx)],
248                 shiftdim (y, 1), szy(2:end));
249      pp.orient = "first";
250
251      if (ispp)
252        yi = pp;
253      else
254        yi = ppval (pp, reshape (xi, szx));
255      endif
256
257    case "*nearest"
258      pp = mkpp ([x(1), x(1)+[0.5:(nx-1)]*dx, x(nx)],
259                 shiftdim (y, 1), szy(2:end));
260      pp.orient = "first";
261
262      if (ispp)
263        yi = pp;
264      else
265        yi = ppval (pp, reshape (xi, szx));
266      endif
267
268    case "previous"
269      pp = mkpp ([x(1:nx); 2*x(nx)-x(nx-1)],
270                 shiftdim (y, 1), szy(2:end));
271      pp.orient = "first";
272
273      if (ispp)
274        yi = pp;
275      else
276        yi = ppval (pp, reshape (xi, szx));
277      endif
278
279    case "*previous"
280      pp = mkpp (x(1)+[0:nx]*dx,
281                 shiftdim (y, 1), szy(2:end));
282      pp.orient = "first";
283
284      if (ispp)
285        yi = pp;
286      else
287        yi = ppval (pp, reshape (xi, szx));
288      endif
289
290    case "linear"
291
292      xx = x;
293      nxx = nx;
294      yy = y;
295      dy = diff (yy);
296      if (have_jumps)
297        ## Omit zero-size intervals.
298        xx(jumps) = [];
299        nxx = rows (xx);
300        yy(jumps, :) = [];
301        dy(jumps, :) = [];
302      endif
303
304      dx = diff (xx);
305      dx = repmat (dx, [1 size(dy)(2:end)]);
306
307      coefs = [(dy./dx).', yy(1:nxx-1, :).'];
308
309      pp = mkpp (xx, coefs, szy(2:end));
310      pp.orient = "first";
311
312      if (ispp)
313        yi = pp;
314      else
315        yi = ppval (pp, reshape (xi, szx));
316      endif
317
318    case "*linear"
319      dy = diff (y);
320      coefs = [(dy/dx).'(:), y(1:nx-1, :).'(:)];
321      pp = mkpp (x, coefs, szy(2:end));
322      pp.orient = "first";
323
324      if (ispp)
325        yi = pp;
326      else
327        yi = ppval (pp, reshape (xi, szx));
328      endif
329
330    case {"pchip", "*pchip", "cubic", "*cubic"}
331      if (nx == 2 || starmethod)
332        x = linspace (x(1), x(nx), ny);
333      endif
334
335      if (ispp)
336        y = shiftdim (reshape (y, szy), 1);
337        yi = pchip (x, y);
338        yi.orient = "first";
339      else
340        y = shiftdim (y, 1);
341        yi = pchip (x, y, reshape (xi, szx));
342        if (! isvector (y))
343          yi = shiftdim (yi, 1);
344        endif
345      endif
346
347    case {"spline", "*spline"}
348      if (nx == 2 || starmethod)
349        x = linspace (x(1), x(nx), ny);
350      endif
351
352      if (ispp)
353        y = shiftdim (reshape (y, szy), 1);
354        yi = spline (x, y);
355        yi.orient = "first";
356      else
357        y = shiftdim (y, 1);
358        yi = spline (x, y, reshape (xi, szx));
359        if (! isvector (y))
360          yi = shiftdim (yi, 1);
361        endif
362      endif
363
364    otherwise
365      error ("interp1: invalid METHOD '%s'", method);
366
367  endswitch
368
369  if (! ispp && isnumeric (extrap))
370    ## determine which values are out of range and set them to extrap,
371    ## unless extrap == "extrap".
372    minx = min (x(1), x(nx));
373    maxx = max (x(1), x(nx));
374
375    xi = reshape (xi, szx);
376    outliers = (xi < minx) | ! (xi <= maxx);  # this even catches NaNs
377    if (size_equal (outliers, yi))
378      yi(outliers) = extrap;
379      yi = reshape (yi, szx);
380    elseif (! isscalar (yi))
381      yi(outliers, :) = extrap;
382    else
383      warning ("interp1: Unreachable state.  Please submit data that produced this warning to bugs.octave.org");
384      yi(outliers.') = extrap;
385    endif
386
387  endif
388
389endfunction
390
391
392%!demo
393%! clf;
394%! xf = 0:0.05:10;  yf = sin (2*pi*xf/5);
395%! xp = 0:10;       yp = sin (2*pi*xp/5);
396%! lin = interp1 (xp,yp,xf, 'linear');
397%! spl = interp1 (xp,yp,xf, 'spline');
398%! pch = interp1 (xp,yp,xf, 'pchip');
399%! near= interp1 (xp,yp,xf, 'nearest');
400%! plot (xf,yf,'r',xf,near,'g',xf,lin,'b',xf,pch,'c',xf,spl,'m',xp,yp,'r*');
401%! legend ('original', 'nearest', 'linear', 'pchip', 'spline');
402%! title ('Interpolation of continuous function sin (x) w/various methods');
403%! %--------------------------------------------------------
404%! % confirm that interpolated function matches the original
405
406%!demo
407%! clf;
408%! xf = 0:0.05:10;  yf = sin (2*pi*xf/5);
409%! xp = 0:10;       yp = sin (2*pi*xp/5);
410%! lin = interp1 (xp,yp,xf, '*linear');
411%! spl = interp1 (xp,yp,xf, '*spline');
412%! pch = interp1 (xp,yp,xf, '*pchip');
413%! near= interp1 (xp,yp,xf, '*nearest');
414%! plot (xf,yf,'r',xf,near,'g',xf,lin,'b',xf,pch,'c',xf,spl,'m',xp,yp,'r*');
415%! legend ('*original', '*nearest', '*linear', '*pchip', '*spline');
416%! title ('Interpolation of continuous function sin (x) w/various *methods');
417%! %--------------------------------------------------------
418%! % confirm that interpolated function matches the original
419
420%!demo
421%! clf;
422%! fstep = @(x) x > 1;
423%! xf = 0:0.05:2;  yf = fstep (xf);
424%! xp = linspace (0,2,10);  yp = fstep (xp);
425%! pch = interp1 (xp,yp,xf, 'pchip');
426%! spl = interp1 (xp,yp,xf, 'spline');
427%! plot (xf,yf,'r',xf,pch,'b',xf,spl,'m',xp,yp,'r*');
428%! title ({'Interpolation of step function with discontinuity at x==1', ...
429%!         'Note: "pchip" is shape-preserving, "spline" (continuous 1st, 2nd derivatives) is not'});
430%! legend ('original', 'pchip', 'spline');
431
432%!demo
433%! clf;
434%! t = 0 : 0.3 : pi; dt = t(2)-t(1);
435%! n = length (t); k = 100; dti = dt*n/k;
436%! ti = t(1) + [0 : k-1]*dti;
437%! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
438%! ddys = diff (diff (interp1 (t,y,ti, 'spline'))./dti)./dti;
439%! ddyp = diff (diff (interp1 (t,y,ti, 'pchip')) ./dti)./dti;
440%! ddyc = diff (diff (interp1 (t,y,ti, 'cubic')) ./dti)./dti;
441%! plot (ti(2:end-1),ddys,'b*', ti(2:end-1),ddyp,'c^', ti(2:end-1),ddyc,'g+');
442%! title ({'Second derivative of interpolated "sin (4*t + 0.3) .* cos (3*t - 0.1)"', ...
443%!         'Note: "spline" has continuous 2nd derivative, others do not'});
444%! legend ('spline', 'pchip', 'cubic');
445
446%!demo
447%! clf;
448%! xf = 0:0.05:10;                yf = sin (2*pi*xf/5) - (xf >= 5);
449%! xp = [0:.5:4.5,4.99,5:.5:10];  yp = sin (2*pi*xp/5) - (xp >= 5);
450%! lin = interp1 (xp,yp,xf, 'linear');
451%! near= interp1 (xp,yp,xf, 'nearest');
452%! plot (xf,yf,'r', xf,near,'g', xf,lin,'b', xp,yp,'r*');
453%! legend ('original', 'nearest', 'linear');
454%! %--------------------------------------------------------
455%! % confirm that interpolated function matches the original
456
457%!demo
458%! clf;
459%! x = 0:0.5:3;
460%! x1 = [3 2 2 1];
461%! x2 = [1 2 2 3];
462%! y1 = [1 1 0 0];
463%! y2 = [0 0 1 1];
464%! h = plot (x, interp1 (x1, y1, x), 'b', x1, y1, 'sb');
465%! hold on
466%! g = plot (x, interp1 (x2, y2, x), 'r', x2, y2, '*r');
467%! axis ([0.5 3.5 -0.5 1.5]);
468%! legend ([h(1), g(1)], {'left-continuous', 'right-continuous'}, ...
469%!         'location', 'northwest')
470%! legend boxoff
471%! %--------------------------------------------------------
472%! % red curve is left-continuous and blue is right-continuous at x = 2
473
474##FIXME: add test for N-d arguments here
475
476## For each type of interpolated test, confirm that the interpolated
477## value at the knots match the values at the knots.  Points away
478## from the knots are requested, but only "nearest" and "linear"
479## confirm they are the correct values.
480
481%!shared xp, yp, xi, style
482%! xp = 0:2:10;
483%! yp = sin (2*pi*xp/5);
484%! xi = [-1, 0, 2.2, 4, 6.6, 10, 11];
485
486## The following BLOCK/ENDBLOCK section is repeated for each style
487##    nearest, previous, next, linear, cubic, spline, pchip
488## The test for ppval of cubic has looser tolerance, but otherwise
489## the tests are identical.
490## Note that the block checks style and *style; if you add more tests
491## be sure to add them to both sections of each block.  One test,
492## style vs. *style, occurs only in the first section.
493## There is an ENDBLOCKTEST after the final block
494
495%!test style = "nearest";
496## BLOCK
497%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
498%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
499%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
500%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
501%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
502%!assert (isempty (interp1 (xp',yp',[],style)))
503%!assert (isempty (interp1 (xp,yp,[],style)))
504%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
505%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
506%!assert (interp1 (xp,yp,xi,style),...
507%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
508%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
509%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
510%!error interp1 (1,1,1, style)
511%!assert (interp1 (xp,[yp',yp'],xi,style),
512%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
513%!test style = ["*",style];
514%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
515%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
516%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
517%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
518%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
519%!assert (isempty (interp1 (xp',yp',[],style)))
520%!assert (isempty (interp1 (xp,yp,[],style)))
521%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
522%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
523%!assert (interp1 (xp,yp,xi,style),...
524%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
525%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
526%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
527%!assert (interp1 (yp, xi, style, 0), ...
528%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
529%!error interp1 (1,1,1, style)
530## ENDBLOCK
531
532%!test style = "previous";
533## BLOCK
534%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
535%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
536%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
537%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
538%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
539%!assert (isempty (interp1 (xp',yp',[],style)))
540%!assert (isempty (interp1 (xp,yp,[],style)))
541%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
542%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
543## This test is expected to fail, so commented out.
544## "previous" and "next" options are not symmetric w.r.t to flipping xp,yp
545#%!assert (interp1 (xp,yp,xi,style),...
546#%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
547%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
548%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
549%!error interp1 (1,1,1, style)
550%!assert (interp1 (xp,[yp',yp'],xi,style),
551%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
552%!test style = ["*",style];
553%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
554%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
555%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
556%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
557%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
558%!assert (isempty (interp1 (xp',yp',[],style)))
559%!assert (isempty (interp1 (xp,yp,[],style)))
560%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
561%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
562#%!assert (interp1 (xp,yp,xi,style),...
563#%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
564%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
565%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
566%!assert (interp1 (yp, xi, style, 0), ...
567%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
568%!error interp1 (1,1,1, style)
569## ENDBLOCK
570
571%!test style = "next";
572## BLOCK
573%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
574%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
575%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
576%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
577%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
578%!assert (isempty (interp1 (xp',yp',[],style)))
579%!assert (isempty (interp1 (xp,yp,[],style)))
580%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
581%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
582#%!assert (interp1 (xp,yp,xi,style),...
583#%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
584%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
585%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
586%!error interp1 (1,1,1, style)
587%!assert (interp1 (xp,[yp',yp'],xi,style),
588%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
589%!test style = ["*",style];
590%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
591%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
592%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
593%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
594%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
595%!assert (isempty (interp1 (xp',yp',[],style)))
596%!assert (isempty (interp1 (xp,yp,[],style)))
597%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
598%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
599#%!assert (interp1 (xp,yp,xi,style),...
600#%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
601%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
602%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
603%!assert (interp1 (yp, xi, style, 0), ...
604%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
605%!error interp1 (1,1,1, style)
606## ENDBLOCK
607
608%!test style = "linear";
609## BLOCK
610%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
611%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
612%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
613%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
614%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
615%!assert (isempty (interp1 (xp',yp',[],style)))
616%!assert (isempty (interp1 (xp,yp,[],style)))
617%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
618%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
619%!assert (interp1 (xp,yp,xi,style),...
620%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
621%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
622%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
623%!error interp1 (1,1,1, style)
624%!assert (interp1 (xp,[yp',yp'],xi,style),
625%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
626%!test style = ['*',style];
627%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
628%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
629%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
630%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
631%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
632%!assert (isempty (interp1 (xp',yp',[],style)))
633%!assert (isempty (interp1 (xp,yp,[],style)))
634%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
635%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
636%!assert (interp1 (xp,yp,xi,style),...
637%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
638%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
639%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
640%!assert (interp1 (yp, xi, style, 0), ...
641%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
642%!assert (interp1 ([1 2 2 3], [1 2 3 4], 2), 3)
643%!assert (interp1 ([3 2 2 1], [4 3 2 1], 2), 2)
644%!error interp1 (1,1,1, style)
645## ENDBLOCK
646
647%!test style = "cubic";
648## BLOCK
649%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
650%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
651%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
652%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
653%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
654%!assert (isempty (interp1 (xp',yp',[],style)))
655%!assert (isempty (interp1 (xp,yp,[],style)))
656%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
657%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
658%!assert (interp1 (xp,yp,xi,style),...
659%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
660%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
661%!        interp1 (xp,yp,xi,style,"extrap"),100*eps)
662%!error interp1 (1,1,1, style)
663%!assert (interp1 (xp,[yp',yp'],xi,style),
664%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
665%!test style = ["*",style];
666%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
667%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
668%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
669%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
670%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
671%!assert (isempty (interp1 (xp',yp',[],style)))
672%!assert (isempty (interp1 (xp,yp,[],style)))
673%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
674%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
675%!assert (interp1 (xp,yp,xi,style),...
676%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
677%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
678%!        interp1 (xp,yp,xi,style,"extrap"),100*eps)
679%!assert (interp1 (yp, xi, style, 0), ...
680%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
681%!error interp1 (1,1,1, style)
682## ENDBLOCK
683
684%!test style = "pchip";
685## BLOCK
686%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
687%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
688%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
689%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
690%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
691%!assert (isempty (interp1 (xp',yp',[],style)))
692%!assert (isempty (interp1 (xp,yp,[],style)))
693%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
694%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
695%!assert (interp1 (xp,yp,xi,style),...
696%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
697%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
698%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
699%!error interp1 (1,1,1, style)
700%!assert (interp1 (xp,[yp',yp'],xi,style),
701%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
702%!test style = ["*",style];
703%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
704%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
705%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
706%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
707%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
708%!assert (isempty (interp1 (xp',yp',[],style)))
709%!assert (isempty (interp1 (xp,yp,[],style)))
710%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
711%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
712%!assert (interp1 (xp,yp,xi,style),...
713%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
714%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
715%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
716%!assert (interp1 (yp, xi, style, 0), ...
717%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
718%!error interp1 (1,1,1, style)
719## ENDBLOCK
720
721%!test style = "spline";
722## BLOCK
723%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
724%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
725%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
726%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
727%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
728%!assert (isempty (interp1 (xp',yp',[],style)))
729%!assert (isempty (interp1 (xp,yp,[],style)))
730%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
731%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
732%!assert (interp1 (xp,yp,xi,style),...
733%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
734%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
735%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
736%!error interp1 (1,1,1, style)
737%!assert (interp1 (xp,[yp',yp'],xi,style),
738%!        interp1 (xp,[yp',yp'],xi,["*",style]),100*eps)
739%!test style = ["*",style];
740%!assert (interp1 (xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA])
741%!assert (interp1 (xp,yp,xp,style), yp, 100*eps)
742%!assert (interp1 (xp,yp,xp',style), yp', 100*eps)
743%!assert (interp1 (xp',yp',xp',style), yp', 100*eps)
744%!assert (interp1 (xp',yp',xp,style), yp, 100*eps)
745%!assert (isempty (interp1 (xp',yp',[],style)))
746%!assert (isempty (interp1 (xp,yp,[],style)))
747%!assert (interp1 (xp,[yp',yp'],xi(:),style),...
748%!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
749%!assert (interp1 (xp,yp,xi,style),...
750%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
751%!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
752%!        interp1 (xp,yp,xi,style,"extrap"),10*eps)
753%!assert (interp1 (yp, xi, style, 0), ...
754%!        interp1 (1:numel (yp), yp, xi, style, 0), 10*eps)
755%!error interp1 (1,1,1, style)
756## ENDBLOCK
757## ENDBLOCKTEST
758
759## test extrapolation
760%!assert (interp1 ([1:5],[3:2:11],[0,6],"linear","extrap"), [1, 13], eps)
761%!assert (interp1 ([1:5],[3:2:11],[0,6],"nearest","extrap"), [3, 11], eps)
762%!assert (interp1 ([1:5],[3:2:11],[0,6],"previous","extrap"), [3, 11], eps)
763%!assert (interp1 ([1:5],[3:2:11],[0,6],"next","extrap"), [3, 11], eps)
764%!assert (interp1 (xp, yp, [-1, max(xp)+1],"linear",5), [5, 5])
765%!assert (interp1 ([0,1],[1,0],[0.1,0.9;0.2,1.1]), [0.9 0.1; 0.8 NA], eps)
766%!assert (interp1 ([0,1],[1,0],[0.1,0.9;0.2,1]), [0.9 0.1; 0.8 0], eps)
767
768## Basic sanity checks
769%!assert (interp1 (1:2,1:2,1.4,"nearest"), 1)
770%!assert (interp1 (1:2,1:2,1.6,"previous"), 1)
771%!assert (interp1 (1:2,1:2,1.4,"next"), 2)
772%!assert (interp1 (1:2,1:2,1.4,"linear"), 1.4)
773%!assert (interp1 (1:4,1:4,1.4,"cubic"), 1.4)
774%!assert (interp1 (1:2,1:2,1.1,"spline"), 1.1)
775%!assert (interp1 (1:3,1:3,1.4,"spline"), 1.4)
776
777%!assert (interp1 (1:2:4,1:2:4,1.4,"*nearest"), 1)
778%!assert (interp1 (1:2:4,1:2:4,2.2,"*previous"), 1)
779%!assert (interp1 (1:2:4,1:2:4,1.4,"*next"), 3)
780%!assert (interp1 (1:2:4,1:2:4,[0,1,1.4,3,4],"*linear"), [NA,1,1.4,3,NA])
781%!assert (interp1 (1:2:8,1:2:8,1.4,"*cubic"), 1.4)
782%!assert (interp1 (1:2,1:2,1.3, "*spline"), 1.3)
783%!assert (interp1 (1:2:6,1:2:6,1.4,"*spline"), 1.4)
784
785%!assert (interp1 ([3,2,1],[3,2,2],2.5), 2.5)
786
787%!assert (interp1 ([4,4,3,2,0],[0,1,4,2,1],[1.5,4,4.5], "linear"), [1.75,1,NA])
788%!assert (interp1 (0:4, 2.5), 1.5)
789
790## Left and Right discontinuities
791%!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap", "right"), [-2,0.5,4,3,1.5])
792%!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap", "left"), [-2,0.5,1,3,1.5])
793
794## Test input validation
795%!error interp1 ()
796%!error interp1 (1,2,3,4,5,6,7)
797%!error <minimum of 2 points required> interp1 (1,1,1, "linear")
798%!error <minimum of 2 points required> interp1 (1,1,1, "*nearest")
799%!error <minimum of 2 points required> interp1 (1,1,1, "*linear")
800%!error <minimum of 2 points required> interp1 (1,1,1, "previous")
801%!error <minimum of 2 points required> interp1 (1,1,1, "*previous")
802%!warning <multiple discontinuities> interp1 ([1 1 1 2], [1 2 3 4], 1);
803%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "next")
804%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "pchip")
805%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "cubic")
806%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "spline")
807%!error <invalid METHOD 'invalid'> interp1 (1:2,1:2,1, "invalid")
808