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