Mercurial > forge

comparison main/general/interp1.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children e6b6904a399e
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 ## Copyright (C) 2000 Paul Kienzle
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, write to the Free Software
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 ## usage: yi = interp1(x, y, xi [, 'method' [, 'extrap']])
18 ##
19 ## Interpolate the function y=f(x) at the points xi. The sample
20 ## points x must be strictly monotonic. If y is a matrix with
21 ## length(x) rows, yi will be a matrix of size rows(xi) by columns(y),
22 ## or its transpose if xi is a row vector.
23 ##
24 ## Method is one of:
25 ## 'nearest': return nearest neighbour.
26 ## 'linear': linear interpolation from nearest neighbours
27 ## 'pchip': piece-wise cubic hermite interpolating polynomial
28 ## 'cubic': cubic interpolation from four nearest neighbours
29 ## 'spline': cubic spline interpolation--smooth first and second
30 ## derivatives throughout the curve
31 ## ['*' method]: same as method, but assumes x is uniformly spaced
32 ## only uses x(1) and x(2); usually faster, never slower
33 ##
34 ## Method defaults to 'linear'.
35 ##
36 ## If extrap is the string 'extrap', then extrapolate values beyond
37 ## the endpoints. If extrap is a number, replace values beyond the
38 ## endpoints with that number. If extrap is missing, assume NaN.
39 ##
40 ## Example:
41 ## xf=[0:0.05:10]; yf = sin(2*pi*xf/5);
42 ## xp=[0:10]; yp = sin(2*pi*xp/5);
43 ## lin=interp1(xp,yp,xf);
44 ## spl=interp1(xp,yp,xf,'spline');
45 ## cub=interp1(xp,yp,xf,'cubic');
46 ## near=interp1(xp,yp,xf,'nearest');
47 ## plot(xf,yf,';original;',xf,lin,';linear;',xf,spl,';spline;',...
48 ## xf,cub,';cubic;',xf,near,';nearest;',xp,yp,'*;;');
49 ##
50 ## See also: interp
51
52 ## 2000-03-25 Paul Kienzle
53 ## added 'nearest' as suggested by Kai Habel
54 ## 2000-07-17 Paul Kienzle
55 ## added '*' methods and matrix y
56 ## check for proper table lengths
57
58 function yi = interp1(x, y, xi, method, extrap)
59
60 if nargin<3 || nargin>5
61 usage("yi = interp1(x, y, xi [, 'method' [, 'extrap']])");
62 endif
63
64 if nargin < 4,
65 method = 'linear';
66 else
67 method = tolower(method);
68 endif
69
70 if nargin < 5
71 extrap = NaN;
72 endif
73
74 ## reshape matrices for convenience
75 x = x(:);
76 if size(y,1)==1, y=y(:); endif
77 transposed = (size(xi,1)==1);
78 xi = xi(:);
79
80 ## determine sizes
81 nx = size(x,1);
82 [ny, nc] = size(y);
83 if (nx < 2 || ny < 2)
84 error ("interp1: table too short");
85 endif
86
87 ## determine which values are out of range and set them to extrap,
88 ## unless extrap=='extrap' in which case, extrapolate them like we
89 ## should have done in the first place.
90 minx = x(1);
91 if (method(1) == '*')
92 dx = x(2) - x(1);
93 maxx = minx + (ny-1)*dx;
94 else
95 maxx = x(nx);
96 endif
97 if strcmp(extrap,"extrap")
98 range=1:nx;
99 else
100 range = find(xi >= minx & xi <= maxx);
101 yi = extrap*ones(size(xi,1), size(y,2));
102 if isempty(range),
103 if transposed, yi = yi.'; endif
104 return;
105 endif
106 xi = xi(range);
107 endif
108
109 if strcmp(method, 'nearest')
110 idx = lookup(0.5*(x(1:nx-1)+x(2:nx)), xi)+1;
111 yi(range,:) = y(idx,:);
112
113 elseif strcmp(method, '*nearest')
114 idx = floor((xi-minx)/dx+1.5);
115 yi(range,:) = y(idx,:);
116
117 elseif strcmp(method, 'linear')
118 ## find the interval containing the test point
119 idx = lookup (x(2:nx-1), xi)+1;
120 # 2:n-1 so that anything beyond the ends
121 # gets dumped into an interval
122
123 ## use the endpoints of the interval to define a line
124 dy = y(2:ny,:) - y(1:ny-1,:);
125 dx = x(2:nx) - x(1:nx-1);
126 s = (xi - x(idx))./dx(idx);
127 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:);
128
129 elseif strcmp(method, '*linear')
130 ## find the interval containing the test point
131 t = (xi - minx)/dx + 1;
132 idx = floor(t);
133
134 ## use the endpoints of the interval to define a line
135 dy = [y(2:ny,:) - y(1:ny-1,:); zeros(1,nc)];
136 s = (t - idx)./dx;
137 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:);
138
139 elseif strcmp(method, 'pchip')
140 if (nx == 2) x = linspace(minx, maxx, ny); endif
141 yi(range,:) = pchip(x, y, xi);
142
143 elseif strcmp(method, 'cubic')
144 if (nx < 4 || ny < 4)
145 error ("interp1: table too short");
146 endif
147 idx = lookup(x(3:nx-2), xi) + 1;
148
149 ## Construct cubic equations for each interval using divided
150 ## differences (computation of c and d don't use divided differences
151 ## but instead solve 2 equations for 2 unknowns). Perhaps
152 ## reformulating this as a lagrange polynomial would be more efficient.
153 i=1:nx-3;
154 J = ones(1,nc);
155 dx = diff(x);
156 dx2 = x(i+1).^2 - x(i).^2;
157 dx3 = x(i+1).^3 - x(i).^3;
158 a=diff(y,3)./dx(i,J).^3/6;
159 b=(diff(y(1:nx-1,:),2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
160 c=(diff(y(1:nx-2,:),1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
161 d=y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
162 yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ...
163 + c(idx,:)).*xi(:,J) + d(idx,:);
164
165 elseif strcmp(method, '*cubic')
166 if (nx < 4 || ny < 4)
167 error ("interp1: table too short");
168 endif
169
170 ## From: Miloje Makivic
171 ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
172 t = (xi - minx)/dx + 1;
173 idx = max(min(floor(t), ny-2), 2);
174 t = t - idx;
175 t2 = t.*t;
176 tp = 1 - 0.5*t;
177 a = (1 - t2).*tp;
178 b = (t2 + t).*tp;
179 c = (t2 - t).*tp/3;
180 d = (t2 - 1).*t/6;
181 J = ones(1,nc);
182 yi(range,:) = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ...
183 + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:);
184
185 elseif strcmp(method, 'spline') || strcmp(method, '*spline')
186 if (nx == 2) x = linspace(minx, maxx, ny); endif
187 yi(range,:) = spline(x, y, xi);
188
189 else
190 error(["interp1 doesn't understand method '", method, "'"]);
191 endif
192 if transposed, yi=yi.'; endif
193
194 endfunction
195
196 %!demo
197 %! xf=0:0.05:10; yf = sin(2*pi*xf/5);
198 %! xp=0:10; yp = sin(2*pi*xp/5);
199 %! lin=interp1(xp,yp,xf,'linear');
200 %! spl=interp1(xp,yp,xf,'spline');
201 %! cub=interp1(xp,yp,xf,'pchip');
202 %! near=interp1(xp,yp,xf,'nearest');
203 %! plot(xf,yf,';original;',xf,near,';nearest;',xf,lin,';linear;',...
204 %! xf,cub,';pchip;',xf,spl,';spline;',xp,yp,'*;;');
205 %! %--------------------------------------------------------
206 %! % confirm that interpolated function matches the original
207
208 %!demo
209 %! xf=0:0.05:10; yf = sin(2*pi*xf/5);
210 %! xp=0:10; yp = sin(2*pi*xp/5);
211 %! lin=interp1(xp,yp,xf,'*linear');
212 %! spl=interp1(xp,yp,xf,'*spline');
213 %! cub=interp1(xp,yp,xf,'*cubic');
214 %! near=interp1(xp,yp,xf,'*nearest');
215 %! plot(xf,yf,';*original;',xf,near,';*nearest;',xf,lin,';*linear;',...
216 %! xf,cub,';*cubic;',xf,spl,';*spline;',xp,yp,'*;;');
217 %! %--------------------------------------------------------
218 %! % confirm that interpolated function matches the original
219
220 %!shared xp, yp, xi, style
221 %! xp=0:5; yp = sin(2*pi*xp/5);
222 %! xi = sort([-1, max(xp)*rand(1,6), max(xp)+1]);
223
224 %!test style = 'nearest';
225 %!assert (interp1(xp, yp, [-1, max(xp)]), [NaN, NaN]);
226 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
227 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
228 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
229 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
230 %!assert (isempty(interp1(xp',yp',[],style)));
231 %!assert (isempty(interp1(xp,yp,[],style)));
232 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
233 %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
234 %!assert (interp1(xp,[yp',yp'],xi,style),
235 %! interp1(xp,[yp',yp'],xi,["*",style]));
236
237 %!test style = 'linear';
238 %!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]);
239 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
240 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
241 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
242 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
243 %!assert (isempty(interp1(xp',yp',[],style)));
244 %!assert (isempty(interp1(xp,yp,[],style)));
245 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
246 %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
247 %!assert (interp1(xp,[yp',yp'],xi,style),
248 %! interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
249
250 %!test style = 'cubic';
251 %!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]);
252 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
253 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
254 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
255 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
256 %!assert (isempty(interp1(xp',yp',[],style)));
257 %!assert (isempty(interp1(xp,yp,[],style)));
258 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
259 %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
260 %!assert (interp1(xp,[yp',yp'],xi,style),
261 %! interp1(xp,[yp',yp'],xi,["*",style]),1000*eps);
262
263 %!test style = 'spline';
264 %!assert (interp1(xp, yp, [-1, max(xp) + 1]), [NaN, NaN]);
265 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
266 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
267 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
268 %!assert (interp1(xp',yp',xp,style), yp, 100*eps);
269 %!assert (isempty(interp1(xp',yp',[],style)));
270 %!assert (isempty(interp1(xp,yp,[],style)));
271 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
272 %! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
273 %!assert (interp1(xp,[yp',yp'],xi,style),
274 %! interp1(xp,[yp',yp'],xi,["*",style]),10*eps);
275
276 %!error interp1
277 %!error interp1(1:2,1:2,1,'bogus')
278
279 %!error interp1(1,1,1, 'nearest');
280 %!assert (interp1(1:2,1:2,1.4,'nearest'),1);
281 %!error interp1(1,1,1, 'linear');
282 %!assert (interp1(1:2,1:2,1.4,'linear'),1.4);
283 %!error interp1(1:3,1:3,1, 'cubic');
284 %!assert (interp1(1:4,1:4,1.4,'cubic'),1.4);
285 %!error interp1(1:2,1:2,1, 'spline');
286 %!assert (interp1(1:3,1:3,1.4,'spline'),1.4);
287
288 %!error interp1(1,1,1, '*nearest');
289 %!assert (interp1(1:2,1:2,1.4,'*nearest'),1);
290 %!error interp1(1,1,1, '*linear');
291 %!assert (interp1(1:2,1:2,1.4,'*linear'),1.4);
292 %!error interp1(1:3,1:3,1, '*cubic');
293 %!assert (interp1(1:4,1:4,1.4,'*cubic'),1.4);
294 %!error interp1(1:2,1:2,1, '*spline');
295 %!assert (interp1(1:3,1:3,1.4,'*spline'),1.4);