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