[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
More text and figures for interpolation chapter
From: |
David Bateman |
Subject: |
More text and figures for interpolation chapter |
Date: |
Thu, 14 Jun 2007 03:51:43 +0200 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
Here is a patch that adds a little more text to the interpolation
chapter. It also adds additional figures to this chapter. It uses a
script to create these in the same manner as the sparse figures.
D.
Index: doc/interpreter/Makefile.in
===================================================================
RCS file: /usr/local/cvsroot/octave/doc/interpreter/Makefile.in,v
retrieving revision 1.147
diff -u -r1.147 Makefile.in
--- doc/interpreter/Makefile.in 13 Jun 2007 05:57:02 -0000 1.147
+++ doc/interpreter/Makefile.in 14 Jun 2007 01:39:36 -0000
@@ -18,7 +18,7 @@
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_DATA = @INSTALL_DATA@
-SCRIPT_SOURCES = sparseimages.m
+SCRIPT_SOURCES = sparseimages.m interpimages.m
EXAMPLE_FILES_NODIR = \
addtwomatrices.cc \
@@ -43,16 +43,20 @@
EXAMPLE_FILES = $(addprefix $(top_srcdir)/examples/, $(EXAMPLE_FILES_NODIR))
-SPARSEIMAGES_1 = gplot grid spmatrix spchol spcholperm
+INTERPIMAGES = interpft interpn interpderiv
+INTERPIMAGES_EPS = $(addsuffix .eps, $(INTERPIMAGES))
+INTERPIMAGES_PDF = $(addsuffix .pdf, $(INTERPIMAGES))
+INTERPIMAGES_PNG = $(addsuffix .png, $(INTERPIMAGES))
+SPARSEIMAGES_1 = gplot grid spmatrix spchol spcholperm
SPARSEIMAGES_EPS = $(addsuffix .eps, $(SPARSEIMAGES_1))
SPARSEIMAGES_PDF = $(addsuffix .pdf, $(SPARSEIMAGES_1))
SPARSEIMAGES_PNG = $(addsuffix .png, $(SPARSEIMAGES_1))
SPARSEIMAGES_TXT = $(addsuffix .txt, $(SPARSEIMAGES_1))
-IMAGES_EPS = $(SPARSEIMAGES_EPS)
-IMAGES_PDF = $(SPARSEIMAGES_PDF)
-IMAGES_PNG = $(SPARSEIMAGES_PNG)
+IMAGES_EPS = $(SPARSEIMAGES_EPS) $(INTERPIMAGES_EPS)
+IMAGES_PDF = $(SPARSEIMAGES_PDF) $(INTERPIMAGES_PDF)
+IMAGES_PNG = $(SPARSEIMAGES_PNG) $(INTERPIMAGES_PNG)
IMAGES_TXT = $(SPARSEIMAGES_TXT)
HTML_IMAGES_PNG = $(addprefix HTML/, $(IMAGES_PNG))
@@ -207,6 +211,9 @@
--eval "$(notdir $(basename $<)) ('$(notdir $(basename $@))', '$(patsubst
.%,%, $(suffix $@))'); sleep (1);"
endef
+$(INTERPIMAGES_EPS) $(INTERPIMAGES_PNG) $(INTERPIMAGES_TXT): interpimages.m
+ $(run-octave)
+
$(SPARSEIMAGES_EPS) $(SPARSEIMAGES_PNG) $(SPARSEIMAGES_TXT): sparseimages.m
$(run-octave)
Index: doc/interpreter/interp.txi
===================================================================
RCS file: /usr/local/cvsroot/octave/doc/interpreter/interp.txi,v
retrieving revision 1.3
diff -u -r1.3 interp.txi
--- doc/interpreter/interp.txi 12 Jun 2007 21:39:27 -0000 1.3
+++ doc/interpreter/interp.txi 14 Jun 2007 01:39:36 -0000
@@ -15,6 +15,43 @@
@DOCSTRING(interp1)
+There are some important differences between the various interpolation
+methods. The 'spline' method enforces that both the first and second
+derivatives of the interpolated values have a continuous derivative,
+whereas the other methods do not. This can be demonstrated by the code
+
address@hidden
address@hidden
+t = 0 : 0.3 : pi; dt = t(2)-t(1);
+n = length (t); k = 100; dti = dt*n/k;
+ti = t(1) + [0 : k-1]*dti;
+y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+ ti(2:end-1),ddyp,'c^');
+legend('cubic','spline','pchip');
address@hidden group
address@hidden example
+
address@hidden
address@hidden
+The result of which can be seen in @ref{fig:interpderiv}.
+
address@hidden Figure,fig:interpderiv
address@hidden,8cm}
address@hidden of second derivative of interpolated values for
+various interpolation methods}
address@hidden float
address@hidden ifnotinfo
+
+This means that in general the 'spline' method results in smooth
+results. If the function to be interpolated is in fact smooth, then
+'spline' will give excellent results. However, if the function to be
+evaluated is in some manner discontinuous, then 'cubic' or 'pchip'
+interpolation might give better results.
+
Fourier interpolation, is a resampling technique where a signal is
converted to the frequency domain, padded with zeros and then
reconverted to the time domain.
@@ -40,8 +77,20 @@
@end group
@end example
address@hidden
which demonstrates the poor behavior of Fourier interpolation for non
periodic functions.
address@hidden ifinfo
address@hidden
+which demonstrates the poor behavior of Fourier interpolation for non
+periodic functions, as can be seen in @ref{fig:interpft}.
+
address@hidden Figure,fig:interpft
address@hidden,8cm}
address@hidden of @code{interp1} and @code{interpft} for non
+periodic data}
address@hidden float
address@hidden ifnotinfo
In additional the support function @code{spline} and @code{lookup} that
underlie the @code{interp1} function can be called directly.
@@ -81,11 +130,12 @@
f = @@(x,y,z) x.^2 - y - z.^2;
[xx, yy, zz] = meshgrid (x, y, z);
v = f (xx,yy,zz);
-xi = yi = zi = -1:0.5:1;
+xi = yi = zi = -1:0.1:1;
[xxi, yyi, zzi] = meshgrid (xi, yi, zi);
-vi = interp3(x, y, z, v, xxi, yyi, zzi);
+vi = interp3(x, y, z, v, xxi, yyi, zzi, 'spline');
[xxi, yyi, zzi] = ndgrid (xi, yi, zi);
-vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
+vi2 = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+mesh (yi, zi, squeeze (vi2(1,:,:)));
@end group
@end example
@@ -93,6 +143,14 @@
where @code{vi} and @code{vi2} are identical. The reversal of the
dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions
respectively.
address@hidden
+The result of this code can be seen in @ref{fig:interpn}.
+
address@hidden Figure,fig:interpn
address@hidden,8cm}
address@hidden of the use of @code{interpn}}
address@hidden float
address@hidden ifnotinfo
In additional the support function @code{bicubic} that underlies the
cubic interpolation of @code{interp2} function can be called directly.
Index: scripts/general/__splinen__.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/__splinen__.m,v
retrieving revision 1.1
diff -u -r1.1 __splinen__.m
--- scripts/general/__splinen__.m 12 Jun 2007 21:39:27 -0000 1.1
+++ scripts/general/__splinen__.m 14 Jun 2007 01:39:37 -0000
@@ -28,14 +28,14 @@
if (nargin != 5)
error ("Incorrect number of arguments");
endif
-
- if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (@isvector, x)) ||
- !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (@isvector, xi)))
+ isvec = @(x) prod(size(x)) == length(x); # ND isvector function
+ if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (isvec, x)) ||
+ !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (isvec, xi)))
error ("%s: non gridded data or dimensions inconsistent", f);
endif
yi = y;
for i = length(x):-1:1
- yi = spline (x{i}, yi, xi{i}).';
+ yi = permute (spline (x{i}, yi, xi{i}), [length(x),1:length(x)-1]);
endfor
[xi{:}] = ndgrid (xi{:});
Index: scripts/general/interp1.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/interp1.m,v
retrieving revision 1.5
diff -u -r1.5 interp1.m
--- scripts/general/interp1.m 12 Jun 2007 21:39:27 -0000 1.5
+++ scripts/general/interp1.m 14 Jun 2007 01:39:37 -0000
@@ -345,6 +345,19 @@
%! %--------------------------------------------------------
%! % confirm that interpolated function matches the original
+%!demo
+%! t = 0 : 0.3 : pi; dt = t(2)-t(1);
+%! n = length (t); k = 100; dti = dt*n/k;
+%! ti = t(1) + [0 : k-1]*dti;
+%! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+%! ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+%! ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+%! ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+%! plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+%! ti(2:end-1),ddyp,'c^');
+%! legend('cubic','spline','pchip');
+%! title("Second derivative of interpolated 'sin (4*t + 0.3) .* cos (3*t -
0.1)'");
+
## For each type of interpolated test, confirm that the interpolated
## value at the knots match the values at the knots. Points away
## from the knots are requested, but only 'nearest' and 'linear'
Index: scripts/general/interpn.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/interpn.m,v
retrieving revision 1.1
diff -u -r1.1 interpn.m
--- scripts/general/interpn.m 12 Jun 2007 21:39:27 -0000 1.1
+++ scripts/general/interpn.m 14 Jun 2007 01:39:37 -0000
@@ -128,22 +128,22 @@
endif
idx (1 : nd) = {1};
idx (i) = ":";
- x{i} = x{i}(idx{:});
+ x{i} = x{i}(idx{:})(:);
endfor
idx (1 : nd) = {1};
idx (1) = ":";
- x{1} = x{1}(idx{:});
+ x{1} = x{1}(idx{:})(:);
endif
if (strcmp (method, "linear") || strcmp (method, "nearest"))
if (all (cellfun (@isvector, y)))
[y{:}] = ndgrid (y{:});
endif
- elseif (any (! cellfun (@isvector, x)))
+ elseif (any (! cellfun (@isvector, y)))
for i = 1 : nd
idx (1 : nd) = {1};
idx (i) = ":";
- y{i} = y{i}(idx{:});
+ y{i} = y{i}(idx{:})(:).';
endfor
endif
@@ -170,7 +170,7 @@
endfor
vi(idx) = extrapval;
vi = reshape (vi, yshape);
- elseif (strcmp (method, "spline"))
+ elseif (strcmp (method, "spline"))
vi = __splinen__ (x, v, y, extrapval, "interpn");
elseif (strcmp (method, "cubic"))
error ("cubic interpolation not yet implemented");
@@ -216,3 +216,14 @@
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!demo
+%! x = y = z = -1:1;
+%! f = @(x,y,z) x.^2 - y - z.^2;
+%! [xx, yy, zz] = meshgrid (x, y, z);
+%! v = f (xx,yy,zz);
+%! xi = yi = zi = -1:0.1:1;
+%! [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
+%! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+%! mesh (yi, zi, squeeze (vi(1,:,:)));
+
Index: scripts/polynomial/mkpp.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/polynomial/mkpp.m,v
retrieving revision 1.3
diff -u -r1.3 mkpp.m
--- scripts/polynomial/mkpp.m 22 Jan 2007 17:29:53 -0000 1.3
+++ scripts/polynomial/mkpp.m 14 Jun 2007 01:39:37 -0000
@@ -52,7 +52,7 @@
d = round (rows (P) / pp.n);
endif
pp.d = d;
- if (pp.n*d != rows (P))
+ if (pp.n * prod (d) != rows (P))
error ("mkpp: num intervals in x doesn't match num polynomials in P");
endif
endfunction
2007-06-14 David Bateman <address@hidden>
* interpreter/Makefile.in (SCRIPT_SORCES): add interimages.m
(INTERPIMAGES*): New variables. Add targets for them
(IMAGES_EPS, IMAGES_PDF, IMAGES_PNG): Add the INTERPIMAGES.
* interpreter/interpimages.m: New function
* interpreter/interp.txi: Add text about second derivation of
splines and add figures.
2007-06-14 David Bateman <address@hidden>
* general/__splinen__.m: Check also for ND vectors. Fix for N > 2,
as permutation of results was incorrect.
* general/interp1.m: Add demo on second derivative
* general/interpn.m: Convert "y" to vectors for __splinen__
call. Add 3D demo.
* polynomial/mkpp.m: Correction for matrices of 3 or more dimensions.
- More text and figures for interpolation chapter,
David Bateman <=