Mercurial > hg > octave-lyh
changeset 6721:01036667884a
[project @ 2007-06-14 06:56:41 by dbateman]
author | dbateman |
---|---|
date | Thu, 14 Jun 2007 06:56:42 +0000 |
parents | fa2f5d4e55db |
children | 5b09d433171c |
files | doc/ChangeLog doc/interpreter/Makefile.in doc/interpreter/interp.txi doc/interpreter/interpimages.m scripts/ChangeLog scripts/general/__splinen__.m scripts/general/interp1.m scripts/general/interpn.m scripts/polynomial/mkpp.m |
diffstat | 9 files changed, 171 insertions(+), 18 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,12 @@ +2007-06-14 David Bateman <dbateman@free.fr> + + * 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-12 David Bateman <dbateman@free.fr> * interpreter/numbers.txi: Document that 64-bit arithmetic is
--- a/doc/interpreter/Makefile.in +++ b/doc/interpreter/Makefile.in @@ -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)) +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)
--- a/doc/interpreter/interp.txi +++ b/doc/interpreter/interp.txi @@ -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 + +@example +@group +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'); +@end group +@end example + +@ifnotinfo +@noindent +The result of which can be seen in @ref{fig:interpderiv}. + +@float Figure,fig:interpderiv +@image{interpderiv,8cm} +@caption{Comparison of second derivative of interpolated values for +various interpolation methods} +@end float +@end 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 +@ifinfo which demonstrates the poor behavior of Fourier interpolation for non periodic functions. +@end ifinfo +@ifnotinfo +which demonstrates the poor behavior of Fourier interpolation for non +periodic functions, as can be seen in @ref{fig:interpft}. + +@float Figure,fig:interpft +@image{interpft,8cm} +@caption{Comparison of @code{interp1} and @code{interpft} for non +periodic data} +@end float +@end 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. +@ifnotinfo +The result of this code can be seen in @ref{fig:interpn}. + +@float Figure,fig:interpn +@image{interpn,8cm} +@caption{Demonstration of the use of @code{interpn}} +@end float +@end ifnotinfo In additional the support function @code{bicubic} that underlies the cubic interpolation of @code{interp2} function can be called directly.
new file mode 100644 --- /dev/null +++ b/doc/interpreter/interpimages.m @@ -0,0 +1,45 @@ +function interpimages (nm, typ) + bury_output (); + if (strcmp (nm, "interpft")) + t = 0 : 0.3 : pi; dt = t(2)-t(1); + n = length (t); k = 100; + ti = t(1) + [0 : k-1]*dt*n/k; + y = sin (4*t + 0.3) .* cos (3*t - 0.1); + yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1); + plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ... + ti, interpft (y, k), 'c', t, y, 'r+'); + legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data'); + print (strcat (nm, ".", typ), strcat ("-d", typ)) + elseif (strcmp (nm, "interpn")) + 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,:,:))); + print (strcat (nm, ".", typ), strcat ("-d", typ)) + elseif (strcmp (nm, "interpderiv")) + 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'); + print (strcat (nm, ".", typ), strcat ("-d", typ)) + endif + bury_output (); +endfunction + +## Use this function before plotting commands and after every call to +## print since print() resets output to stdout (unfortunately, gnpulot +## can't pop output as it can the terminal type). +function bury_output () + f = figure (1); + set (f, "visible", "off"); +endfunction
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,13 @@ +2007-06-14 David Bateman <dbateman@free.fr> + + * 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. + 2007-06-13 John W. Eaton <jwe@octave.org> * miscellaneous/mkoctfile.m: Quote args too.
--- a/scripts/general/__splinen__.m +++ b/scripts/general/__splinen__.m @@ -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{:});
--- a/scripts/general/interp1.m +++ b/scripts/general/interp1.m @@ -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'
--- a/scripts/general/interpn.m +++ b/scripts/general/interpn.m @@ -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,:,:))); +