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