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