1 ## Copyright (C) 2001-2012 Kai Habel
3 ## This file is part of Octave.
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.
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.
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/>.
20 ## @deftypefn {Function File} {@var{pp} =} pchip (@var{x}, @var{y})
21 ## @deftypefnx {Function File} {@var{yi} =} pchip (@var{x}, @var{y}, @var{xi})
22 ## Return the Piecewise Cubic Hermite Interpolating Polynomial (pchip) of
23 ## points @var{x} and @var{y}.
25 ## If called with two arguments, return the piecewise polynomial @var{pp}
26 ## that may be used with @code{ppval} to evaluate the polynomial at specific
27 ## points. When called with a third input argument, @code{pchip} evaluates
28 ## the pchip polynomial at the points @var{xi}. The third calling form is
29 ## equivalent to @code{ppval (pchip (@var{x}, @var{y}), @var{xi})}.
31 ## The variable @var{x} must be a strictly monotonic vector (either
32 ## increasing or decreasing) of length @var{n}. @var{y} can be either a
33 ## vector or array. If @var{y} is a vector then it must be the same length
34 ## @var{n} as @var{x}. If @var{y} is an array then the size of @var{y} must
37 ## $$[s_1, s_2, \cdots, s_k, n]$$
40 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
42 ## The array is reshaped internally to a matrix where the leading
43 ## dimension is given by
45 ## $$s_1 s_2 \cdots s_k$$
48 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
50 ## and each row of this matrix is then treated separately. Note that this
51 ## is exactly opposite to @code{interp1} but is done for @sc{matlab}
54 ## @seealso{spline, ppval, mkpp, unmkpp}
57 ## Author: Kai Habel <kai.habel@gmx.de>
60 ## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom)
64 ## S_k(x_k+1) = y_k+1;
66 ## S_k'(x_k+1) = y_k+1';
68 function ret = pchip (x, y, xi)
70 if (nargin < 2 || nargin > 3)
78 ## Check the size and shape of y
80 y = y(:).'; ##row vector
82 if !(size_equal (x, y))
83 error ("pchip: length of X and Y must match")
88 error ("pchip: length of X and last dimension of Y must match")
90 y = reshape (y, [prod(szy(1:end-1)), szy(end)]);
99 error("pchip: X must be strictly monotonic");
104 ## Compute derivatives.
105 d = __pchip_deriv__ (x, y, 2);
109 ## This is taken from SLATEC.
112 delta = diff (y, 1, 2) / h;
113 del1 = (d1 - delta) / h;
114 del2 = (d2 - delta) / h;
118 coeffs = cat (3, c3, c2, d1, f1);
120 ret = mkpp (x, coeffs, szy(1:end-1));
123 ret = ppval (ret, xi);
130 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
132 %! yspline = spline(x,y,xi);
133 %! ypchip = pchip(x,y,xi);
134 %! title("pchip and spline fit to discontinuous function");
135 %! plot(xi,yspline,xi,ypchip,"-",x,y,"+");
136 %! legend ("spline","pchip","data");
137 %! %-------------------------------------------------------------------
138 %! % confirm that pchip agreed better to discontinuous data than spline
140 %!shared x,y,y2,pp,yi1,yi2,yi3
142 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
143 %!assert (pchip(x,y,x), y);
144 %!assert (pchip(x,y,x'), y');
145 %!assert (pchip(x',y',x'), y');
146 %!assert (pchip(x',y',x), y);
147 %!assert (isempty(pchip(x',y',[])));
148 %!assert (isempty(pchip(x,y,[])));
149 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
150 %!assert (pchip(x,[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
151 %!assert (pchip(x',[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
152 %!assert (pchip(x',[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
154 %! x=(0:8)*pi/4;y=[sin(x);cos(x)];
155 %! y2(:,:,1)=y;y2(:,:,2)=y+1;y2(:,:,3)=y-1;
156 %! pp=pchip(x,shiftdim(y2,2));
157 %! yi1=ppval(pp,(1:4)*pi/4);
158 %! yi2=ppval(pp,repmat((1:4)*pi/4,[5,1]));
159 %! yi3=ppval(pp,[pi/2,pi]);
160 %!assert(size(pp.coefs),[48,4]);
161 %!assert(pp.pieces,8);
162 %!assert(pp.order,4);
163 %!assert(pp.dim,[3,2]);
164 %!assert(ppval(pp,pi),[0,-1;1,0;-1,-2],1e-14);
165 %!assert(yi3(:,:,2),ppval(pp,pi),1e-14);
166 %!assert(yi3(:,:,1),[1,0;2,1;0,-1],1e-14);
167 %!assert(squeeze(yi1(1,2,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
168 %!assert(size(yi2),[3,2,5,4]);
169 %!assert(squeeze(yi2(1,2,3,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
171 %!error (pchip (1,2));
172 %!error (pchip (1,2,3));