Mercurial > hg > octave-lyh

comparison scripts/general/interp1.m @ 5838:376e02b2ce70

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-06-01 20:23:53 by jwe]
author jwe
date Thu, 01 Jun 2006 20:23:54 +0000
parents 55404f3b0da1
children 06f26e174fc9
comparison
equal deleted inserted replaced
5837:55404f3b0da1 5838:376e02b2ce70
76 ## @end example 76 ## @end example
77 ## 77 ##
78 ## @seealso{interpft} 78 ## @seealso{interpft}
79 ## @end deftypefn 79 ## @end deftypefn
80 80
81 ## 2000-03-25 Paul Kienzle 81 ## Author: Paul Kienzle
82 ## Date: 2000-03-25
82 ## added 'nearest' as suggested by Kai Habel 83 ## added 'nearest' as suggested by Kai Habel
83 ## 2000-07-17 Paul Kienzle 84 ## 2000-07-17 Paul Kienzle
84 ## added '*' methods and matrix y 85 ## added '*' methods and matrix y
85 ## check for proper table lengths 86 ## check for proper table lengths
86 ## 2002-01-23 Paul Kienzle 87 ## 2002-01-23 Paul Kienzle
87 ## fixed extrapolation 88 ## fixed extrapolation
88 89
89 function yi = interp1(x, y, varargin) 90 function yi = interp1 (x, y, varargin)
90 91
91 if ( nargin < 3 || nargin > 6) 92 if (nargin < 3 || nargin > 6)
92 print_usage (); 93 print_usage ();
93 endif 94 endif
94 95
95 method = "linear"; 96 method = "linear";
96 extrap = NaN; 97 extrap = NaN;
97 xi = []; 98 xi = [];
98 pp = false; 99 pp = false;
99 firstnumeric = true; 100 firstnumeric = true;
100 101
101 if (nargin > 2) 102 if (nargin > 2)
102 for i = 1:length(varargin) 103 for i = 1:length (varargin)
103 arg = varargin{i}; 104 arg = varargin{i};
104 if (ischar(arg)) 105 if (ischar (arg))
105 arg = tolower (arg); 106 arg = tolower (arg);
106 if (strcmp("extrap",arg)) 107 if (strcmp ("extrap", arg))
107 extrap = "extrap"; 108 extrap = "extrap";
108 elseif (strcmp("pp",arg)) 109 elseif (strcmp ("pp", arg))
109 pp = true; 110 pp = true;
110 else 111 else
111 method = arg; 112 method = arg;
112 endif 113 endif
113 else 114 else
121 endfor 122 endfor
122 endif 123 endif
123 124
124 ## reshape matrices for convenience 125 ## reshape matrices for convenience
125 x = x(:); 126 x = x(:);
126 nx = size(x,1); 127 nx = size (x, 1);
127 if (isvector(y) && size (y, 1) == 1) 128 if (isvector(y) && size (y, 1) == 1)
128 y = y (:); 129 y = y(:);
129 endif 130 endif
130 ndy = ndims (y); 131 ndy = ndims (y);
131 szy = size(y); 132 szy = size (y);
132 ny = szy(1); 133 ny = szy(1);
133 nc = prod (szy(2:end)); 134 nc = prod (szy(2:end));
134 y = reshape (y, ny, nc); 135 y = reshape (y, ny, nc);
135 szx = size(xi); 136 szx = size (xi);
136 xi = xi(:); 137 xi = xi(:);
137 138
138 ## determine sizes 139 ## determine sizes
139 if (nx < 2 || ny < 2) 140 if (nx < 2 || ny < 2)
140 error ("interp1: table too short"); 141 error ("interp1: table too short");
141 endif 142 endif
142 143
143 ## determine which values are out of range and set them to extrap, 144 ## determine which values are out of range and set them to extrap,
144 ## unless extrap=="extrap" in which case, extrapolate them like we 145 ## unless extrap=="extrap" in which case, extrapolate them like we
145 ## should be doing in the first place. 146 ## should be doing in the first place.
148 dx = x(2) - x(1); 149 dx = x(2) - x(1);
149 maxx = minx + (ny-1)*dx; 150 maxx = minx + (ny-1)*dx;
150 else 151 else
151 maxx = x(nx); 152 maxx = x(nx);
152 endif 153 endif
153 if (!pp) 154
154 if ischar(extrap) && strcmp(extrap,"extrap") 155 if (! pp)
155 range=1:size(xi,1); 156 if (ischar (extrap) && strcmp (extrap, "extrap"))
156 yi = zeros(size(xi,1), size(y,2)); 157 range = 1:size (xi, 1);
157 else 158 yi = zeros (size (xi, 1), size (y, 2));
158 range = find(xi >= minx & xi <= maxx); 159 else
159 yi = extrap*ones(size(xi,1), size(y,2)); 160 range = find (xi >= minx & xi <= maxx);
160 if isempty(range), 161 yi = extrap * ones (size (xi, 1), size (y, 2));
161 if (!isvector(y) && length(szx) == 2 && (szx(1) == 1 || szx(2) == 1)) 162 if (isempty (range))
163 if (! isvector (y) && length (szx) == 2
164 && (szx(1) == 1 || szx(2) == 1))
162 if (szx(1) == 1) 165 if (szx(1) == 1)
163 yi = reshape (yi, [szx(2), szy(2:end)]); 166 yi = reshape (yi, [szx(2), szy(2:end)]);
164 else 167 else
165 yi = reshape (yi, [szx(1), szy(2:end)]); 168 yi = reshape (yi, [szx(1), szy(2:end)]);
166 endif 169 endif
171 endif 174 endif
172 xi = xi(range); 175 xi = xi(range);
173 endif 176 endif
174 endif 177 endif
175 178
176 if strcmp(method, "nearest") 179 if (strcmp (method, "nearest"))
177 if (pp) 180 if (pp)
178 yi = mkpp ([x(1);(x(1:end-1)+x(2:end))/2;x(end)],y,szy(2:end)); 181 yi = mkpp ([x(1); (x(1:end-1)+x(2:end))/2; x(end)], y, szy(2:end));
179 else 182 else
180 idx = lookup(0.5*(x(1:nx-1)+x(2:nx)), xi)+1; 183 idx = lookup (0.5*(x(1:nx-1)+x(2:nx)), xi) + 1;
181 yi(range,:) = y(idx,:); 184 yi(range,:) = y(idx,:);
182 endif 185 endif
183 elseif strcmp(method, "*nearest") 186 elseif (strcmp (method, "*nearest"))
184 if (pp) 187 if (pp)
185 yi = mkpp ([minx; minx + [0.5:(ny-1)]'*dx; maxx],y,szy(2:end)); 188 yi = mkpp ([minx; minx+[0.5:(ny-1)]'*dx; maxx], y, szy(2:end));
186 else 189 else
187 idx = max(1,min(ny,floor((xi-minx)/dx+1.5))); 190 idx = max (1, min (ny, floor((xi-minx)/dx+1.5)));
188 yi(range,:) = y(idx,:); 191 yi(range,:) = y(idx,:);
189 endif 192 endif
190 elseif strcmp(method, "linear") 193 elseif (strcmp (method, "linear"))
191 dy = y(2:ny,:) - y(1:ny-1,:); 194 dy = y(2:ny,:) - y(1:ny-1,:);
192 dx = x(2:nx) - x(1:nx-1); 195 dx = x(2:nx) - x(1:nx-1);
193 if (pp) 196 if (pp)
194 yi = mkpp(x, [dy./dx, y(1:end-1)],szy(2:end)); 197 yi = mkpp (x, [dy./dx, y(1:end-1)], szy(2:end));
195 else 198 else
196 ## find the interval containing the test point 199 ## find the interval containing the test point
197 idx = lookup (x(2:nx-1), xi)+1; 200 idx = lookup (x(2:nx-1), xi)+1;
198 # 2:n-1 so that anything beyond the ends 201 # 2:n-1 so that anything beyond the ends
199 # gets dumped into an interval 202 # gets dumped into an interval
200 ## use the endpoints of the interval to define a line 203 ## use the endpoints of the interval to define a line
201 s = (xi - x(idx))./dx(idx); 204 s = (xi - x(idx))./dx(idx);
202 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); 205 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:);
203 endif 206 endif
204 elseif strcmp(method, "*linear") 207 elseif (strcmp (method, "*linear"))
205 if (pp) 208 if (pp)
206 dy = [y(2:ny,:) - y(1:ny-1,:)]; 209 dy = [y(2:ny,:) - y(1:ny-1,:)];
207 yi = mkpp (minx + [0:ny-1]*dx, [dy ./ dx, y(1:end-1)], szy(2:end)); 210 yi = mkpp (minx + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end));
208 else 211 else
209 ## find the interval containing the test point 212 ## find the interval containing the test point
210 t = (xi - minx)/dx + 1; 213 t = (xi - minx)/dx + 1;
211 idx = max(1,min(ny,floor(t))); 214 idx = max(1,min(ny,floor(t)));
212 215
213 ## use the endpoints of the interval to define a line 216 ## use the endpoints of the interval to define a line
214 dy = [y(2:ny,:) - y(1:ny-1,:); y(ny,:) - y(ny-1,:)]; 217 dy = [y(2:ny,:) - y(1:ny-1,:); y(ny,:) - y(ny-1,:)];
215 s = t - idx; 218 s = t - idx;
216 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); 219 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:);
217 endif 220 endif
218 elseif strcmp(method, "pchip") || strcmp(method, "*pchip") 221 elseif (strcmp (method, "pchip") || strcmp (method, "*pchip"))
219 if (nx == 2 || method(1) == "*") 222 if (nx == 2 || method(1) == "*")
220 x = linspace(minx, maxx, ny); 223 x = linspace (minx, maxx, ny);
221 endif 224 endif
222 ## Note that pchip's arguments are transposed relative to interp1 225 ## Note that pchip's arguments are transposed relative to interp1
223 if (pp) 226 if (pp)
224 yi = pchip(x.', y.'); 227 yi = pchip (x.', y.');
225 yi.d = szy(2:end); 228 yi.d = szy(2:end);
226 else 229 else
227 yi(range,:) = pchip(x.', y.', xi.').'; 230 yi(range,:) = pchip (x.', y.', xi.').';
228 endif 231 endif
229 232
230 elseif strcmp(method, "cubic") || (strcmp(method, "*cubic") && pp) 233 elseif (strcmp (method, "cubic") || (strcmp (method, "*cubic") && pp))
231 ## FIXME Is there a better way to treat pp return return and *cubic 234 ## FIXME Is there a better way to treat pp return return and *cubic
232 if (method(1) == "*") 235 if (method(1) == "*")
233 x = linspace(minx, maxx, ny).'; 236 x = linspace (minx, maxx, ny).';
234 nx = ny; 237 nx = ny;
235 endif 238 endif
236 239
237 if (nx < 4 || ny < 4) 240 if (nx < 4 || ny < 4)
238 error ("interp1: table too short"); 241 error ("interp1: table too short");
239 endif 242 endif
240 idx = lookup(x(3:nx-2), xi) + 1; 243 idx = lookup (x(3:nx-2), xi) + 1;
241 244
242 ## Construct cubic equations for each interval using divided 245 ## Construct cubic equations for each interval using divided
243 ## differences (computation of c and d don't use divided differences 246 ## differences (computation of c and d don't use divided differences
244 ## but instead solve 2 equations for 2 unknowns). Perhaps 247 ## but instead solve 2 equations for 2 unknowns). Perhaps
245 ## reformulating this as a lagrange polynomial would be more efficient. 248 ## reformulating this as a lagrange polynomial would be more efficient.
246 i=1:nx-3; 249 i = 1:nx-3;
247 J = ones(1,nc); 250 J = ones (1, nc);
248 dx = diff(x); 251 dx = diff (x);
249 dx2 = x(i+1).^2 - x(i).^2; 252 dx2 = x(i+1).^2 - x(i).^2;
250 dx3 = x(i+1).^3 - x(i).^3; 253 dx3 = x(i+1).^3 - x(i).^3;
251 a=diff(y,3)./dx(i,J).^3/6; 254 a = diff (y, 3)./dx(i,J).^3/6;
252 b=(diff(y(1:nx-1,:),2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2; 255 b = (diff (y(1:nx-1,:), 2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
253 c=(diff(y(1:nx-2,:),1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J); 256 c = (diff (y(1:nx-2,:), 1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
254 d=y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J); 257 d = y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
255 258
256 if (pp) 259 if (pp)
257 xs = [x(1);x(3:nx-2)]; 260 xs = [x(1);x(3:nx-2)];
258 yi = mkpp ([x(1);x(3:nx-2);x(nx)], 261 yi = mkpp ([x(1);x(3:nx-2);x(nx)],
259 [a(:), (b(:) + 3.*xs(:,J).*a(:)), ... 262 [a(:), (b(:) + 3.*xs(:,J).*a(:)), ...
262 xs(:,J).^3.*a(:))], szy(2:end)); 265 xs(:,J).^3.*a(:))], szy(2:end));
263 else 266 else
264 yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ... 267 yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ...
265 + c(idx,:)).*xi(:,J) + d(idx,:); 268 + c(idx,:)).*xi(:,J) + d(idx,:);
266 endif 269 endif
267 elseif strcmp(method, "*cubic") 270 elseif (strcmp (method, "*cubic"))
268 if (nx < 4 || ny < 4) 271 if (nx < 4 || ny < 4)
269 error ("interp1: table too short"); 272 error ("interp1: table too short");
270 endif 273 endif
271 274
272 ## From: Miloje Makivic 275 ## From: Miloje Makivic
273 ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html 276 ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
274 t = (xi - minx)/dx + 1; 277 t = (xi - minx)/dx + 1;
275 idx = max(min(floor(t), ny-2), 2); 278 idx = max (min (floor (t), ny-2), 2);
276 t = t - idx; 279 t = t - idx;
277 t2 = t.*t; 280 t2 = t.*t;
278 tp = 1 - 0.5*t; 281 tp = 1 - 0.5*t;
279 a = (1 - t2).*tp; 282 a = (1 - t2).*tp;
280 b = (t2 + t).*tp; 283 b = (t2 + t).*tp;
281 c = (t2 - t).*tp/3; 284 c = (t2 - t).*tp/3;
282 d = (t2 - 1).*t/6; 285 d = (t2 - 1).*t/6;
283 J = ones(1,nc); 286 J = ones (1, nc);
284 287
285 yi(range,:) = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ... 288 yi(range,:) = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ...
286 + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:); 289 + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:);
287 290
288 elseif strcmp(method, "spline") || strcmp(method, "*spline") 291 elseif (strcmp (method, "spline") || strcmp (method, "*spline"))
289 if (nx == 2 || method(1) == "*") 292 if (nx == 2 || method(1) == "*")
290 x = linspace(minx, maxx, ny); 293 x = linspace(minx, maxx, ny);
291 endif 294 endif
292 ## Note that spline's arguments are transposed relative to interp1 295 ## Note that spline's arguments are transposed relative to interp1
293 if (pp) 296 if (pp)
294 yi = spline(x.', y.'); 297 yi = spline (x.', y.');
295 yi.d = szy(2:end); 298 yi.d = szy(2:end);
296 else 299 else
297 yi(range,:) = spline(x.', y.', xi.').'; 300 yi(range,:) = spline (x.', y.', xi.').';
298 endif 301 endif
299 else 302 else
300 error(["interp1 doesn't understand method '", method, "'"]); 303 error ("interp1: invalid method '%s'", method);
301 endif 304 endif
302 305
303 if (!pp) 306 if (! pp)
304 if (!isvector(y) && length(szx) == 2 && (szx(1) == 1 || szx(2) == 1)) 307 if (! isvector (y) && length (szx) == 2 && (szx(1) == 1 || szx(2) == 1))
305 if (szx(1) == 1) 308 if (szx(1) == 1)
306 yi = reshape (yi, [szx(2), szy(2:end)]); 309 yi = reshape (yi, [szx(2), szy(2:end)]);
307 else 310 else
308 yi = reshape (yi, [szx(1), szy(2:end)]); 311 yi = reshape (yi, [szx(1), szy(2:end)]);
309 endif 312 endif