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