1@c Copyright (C) 2007-2019 John W. Eaton 2@c 3@c This file is part of Octave. 4@c 5@c Octave is free software: you can redistribute it and/or modify it 6@c under the terms of the GNU General Public License as published by 7@c the Free Software Foundation, either version 3 of the License, or 8@c (at your option) any later version. 9@c 10@c Octave is distributed in the hope that it will be useful, but 11@c WITHOUT ANY WARRANTY; without even the implied warranty of 12@c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13@c GNU General Public License for more details. 14@c 15@c You should have received a copy of the GNU General Public License 16@c along with Octave; see the file COPYING. If not, see 17@c <https://www.gnu.org/licenses/>. 18 19@node Interpolation 20@chapter Interpolation 21 22@menu 23* One-dimensional Interpolation:: 24* Multi-dimensional Interpolation:: 25@end menu 26 27@node One-dimensional Interpolation 28@section One-dimensional Interpolation 29 30Octave supports several methods for one-dimensional interpolation, most 31of which are described in this section. @ref{Polynomial Interpolation} 32and @ref{Interpolation on Scattered Data} describe additional methods. 33 34@DOCSTRING(interp1) 35 36There are some important differences between the various interpolation 37methods. The @qcode{"spline"} method enforces that both the first and second 38derivatives of the interpolated values have a continuous derivative, 39whereas the other methods do not. This means that the results of the 40@qcode{"spline"} method are generally smoother. If the function to be 41interpolated is in fact smooth, then @qcode{"spline"} will give excellent 42results. However, if the function to be evaluated is in some manner 43discontinuous, then @qcode{"pchip"} interpolation might give better results. 44 45This can be demonstrated by the code 46 47@example 48@group 49t = -2:2; 50dt = 1; 51ti =-2:0.025:2; 52dti = 0.025; 53y = sign (t); 54ys = interp1 (t,y,ti,"spline"); 55yp = interp1 (t,y,ti,"pchip"); 56ddys = diff (diff (ys)./dti) ./ dti; 57ddyp = diff (diff (yp)./dti) ./ dti; 58figure (1); 59plot (ti,ys,"r-", ti,yp,"g-"); 60legend ("spline", "pchip", 4); 61figure (2); 62plot (ti,ddys,"r+", ti,ddyp,"g*"); 63legend ("spline", "pchip"); 64@end group 65@end example 66 67@ifnotinfo 68@noindent 69The result of which can be seen in @ref{fig:interpderiv1} and 70@ref{fig:interpderiv2}. 71 72@float Figure,fig:interpderiv1 73@center @image{interpderiv1,4in} 74@caption{Comparison of @qcode{"pchip"} and @qcode{"spline"} interpolation methods for a 75step function} 76@end float 77 78@float Figure,fig:interpderiv2 79@center @image{interpderiv2,4in} 80@caption{Comparison of the second derivative of the @qcode{"pchip"} and @qcode{"spline"} 81interpolation methods for a step function} 82@end float 83@end ifnotinfo 84 85Fourier interpolation, is a resampling technique where a signal is 86converted to the frequency domain, padded with zeros and then 87reconverted to the time domain. 88 89@DOCSTRING(interpft) 90 91There are two significant limitations on Fourier interpolation. First, 92the function signal is assumed to be periodic, and so non-periodic 93signals will be poorly represented at the edges. Second, both the 94signal and its interpolation are required to be sampled at equispaced 95points. An example of the use of @code{interpft} is 96 97@example 98@group 99t = 0 : 0.3 : pi; dt = t(2)-t(1); 100n = length (t); k = 100; 101ti = t(1) + [0 : k-1]*dt*n/k; 102y = sin (4*t + 0.3) .* cos (3*t - 0.1); 103yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1); 104plot (ti, yp, "g", ti, interp1 (t, y, ti, "spline"), "b", ... 105 ti, interpft (y, k), "c", t, y, "r+"); 106legend ("sin(4t+0.3)cos(3t-0.1)", "spline", "interpft", "data"); 107@end group 108@end example 109 110@noindent 111@ifinfo 112which demonstrates the poor behavior of Fourier interpolation for non-periodic 113functions. 114@end ifinfo 115@ifnotinfo 116which demonstrates the poor behavior of Fourier interpolation for non-periodic 117functions, as can be seen in @ref{fig:interpft}. 118 119@float Figure,fig:interpft 120@center @image{interpft,4in} 121@caption{Comparison of @code{interp1} and @code{interpft} for non-periodic data} 122@end float 123@end ifnotinfo 124 125In addition, the support functions @code{spline} and @code{lookup} that 126underlie the @code{interp1} function can be called directly. 127 128@DOCSTRING(spline) 129 130@node Multi-dimensional Interpolation 131@section Multi-dimensional Interpolation 132 133There are three multi-dimensional interpolation functions in Octave, with 134similar capabilities. Methods using Delaunay tessellation are described 135in @ref{Interpolation on Scattered Data}. 136 137@DOCSTRING(interp2) 138 139@DOCSTRING(interp3) 140 141@DOCSTRING(interpn) 142 143A significant difference between @code{interpn} and the other two 144multi-dimensional interpolation functions is the fashion in which the 145dimensions are treated. For @code{interp2} and @code{interp3}, the y-axis is 146considered to be the columns of the matrix, whereas the x-axis corresponds to 147the rows of the array. As Octave indexes arrays in column major order, the 148first dimension of any array is the columns, and so @code{interpn} effectively 149reverses the 'x' and 'y' dimensions. Consider the example, 150 151@example 152@group 153x = y = z = -1:1; 154f = @@(x,y,z) x.^2 - y - z.^2; 155[xx, yy, zz] = meshgrid (x, y, z); 156v = f (xx,yy,zz); 157xi = yi = zi = -1:0.1:1; 158[xxi, yyi, zzi] = meshgrid (xi, yi, zi); 159vi = interp3 (x, y, z, v, xxi, yyi, zzi, "spline"); 160[xxi, yyi, zzi] = ndgrid (xi, yi, zi); 161vi2 = interpn (x, y, z, v, xxi, yyi, zzi, "spline"); 162mesh (zi, yi, squeeze (vi2(1,:,:))); 163@end group 164@end example 165 166@noindent 167where @code{vi} and @code{vi2} are identical. The reversal of the 168dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions 169respectively. 170@ifnotinfo 171The result of this code can be seen in @ref{fig:interpn}. 172 173@float Figure,fig:interpn 174@center @image{interpn,4in} 175@caption{Demonstration of the use of @code{interpn}} 176@end float 177@end ifnotinfo 178