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