%% Interpolation with the interp Family of Functions
%% The interp1 Function
% This example shows how the interp1 function can be used to interpolate a set of sample values using the 'pchip' method.
% The function interp1 performs one-dimensional interpolation. Its most general form is:
% Vq = interp1(X,V,Xq,method)
% where X is a vector of coordinates and V is a vector containing the values at those coordinates. Xq is a vector containing the query points at which to interpolate, and method optionally specifies any of four interpolation methods: 'nearest', 'linear', 'pchip', or 'spline'.
% Create a set of 1-D grid points X and corresponding sample values V.

X = [1 2 3 4 5];
V = [12 16 31 10 6];

%
% Interpolate over finer intervals with 0.1 spacing.

Xq = (1:0.1:5);
Vq = interp1(X,V,Xq,'pchip');

%
% Plot the samples and interpolated values.

plot(X,V,'o');
hold on
plot(Xq,Vq,'-');
legend('samples','pchip');
hold off

%
%% 1-D Extrapolation With interp1
% This example shows how to use the 'extrap' option to interpolate beyond the domain of your sample points.
% Define the sample points and values.

X = [1 2 3 4 5];
V = [12 16 31 10 6];

%
% Specify the query points, Xq, that extend beyond the domain of X.

Xq = (0:0.1:6);
Vq = interp1(X,V,Xq,'pchip','extrap');

%
% Plot the results.

figure
plot(X,V,'o');
hold on
plot(Xq,Vq,'-');
legend('samples','pchip');
hold off

%
% In an alternative approach, you can introduce additional points to gain more control over the behavior in the extrapolated regions. For example, you can constrain the curve to remain flat in the extrapolated region by extending the domain with repeated values.

X = [0 1 2 3 4 5 6];
V = [12 12 16 31 10 6 6];

%
% Specify the query points, Xq, that extend further beyond the domain of X.

Xq = (-1:0.1:7);

%
% Interpolate using 'pchip'. You can omit the 'extrap' option because it is the default with the 'pchip', 'makima', and 'spline' methods.

Vq = interp1(X,V,Xq,'pchip');

%
% Plot the results.

figure
plot(X,V,'o');
hold on
plot(Xq,Vq,'-');
legend('samples','pchip');
hold off

%
%% The interp2 Function
% This example shows how the interp2 function can be used to interpolate the coarsely sampled peaks function over a finer grid.
% The interp2 and interp3 functions perform two and three-dimensional interpolation respectively, and they interpolate grids in the meshgrid format. The calling syntax for interp2 has the general form:
% Vq = interp2(X,Y,V,Xq,Yq,method)
% where X and Y are arrays of coordinates that define a grid in meshgrid format, and V is an array containing the values at the grid points. Xq and Yq are arrays containing the coordinates of the query points at which to interpolate. method optionally specifies any of four interpolation methods: 'nearest', 'linear', 'cubic', or 'spline'.
% The grid points that comprise X and Y must be monotonically increasing and should conform to the meshgrid format.
% Create a coarse grid and corresponding sample values.

[X,Y] = meshgrid(-3:1:3);
V = peaks(X,Y);

%
% Plot the sample values.

surf(X,Y,V)
title('Sample Grid');

%
% Generate a finer grid for interpolation.

[Xq,Yq] = meshgrid(-3:0.25:3);

%
% Use interp2 to interpolate at the query points.

Vq = interp2(X,Y,V,Xq,Yq,'linear');

%
% Plot the results.

surf(Xq,Yq,Vq);
title('Refined Grid');

%
%% The interp3 Function
% This example shows how to use interp3 to interpolate a 3-D function at a single query point and compare it to the value generated by an analytic expression.
% interp3 works in the same way as interp2 except that it takes two additional arguments: one for the third dimension in the sample grid and the other for the third dimension in the query points,
% Vq = interp3(X,Y,Z,V,Xq,Yq,Zq,method).
% As is the case with interp2, the grid points you supply to interp3 must be monotonically increasing and should conform to the meshgrid format.
% Define a function that generates values for X, Y, and Z input.

generatedvalues = @(X,Y,Z)(X.^2 + Y.^3 + Z.^4);

%
% Create the sample data.

[X,Y,Z] = meshgrid((-5:.25:5));
V = generatedvalues(X,Y,Z);

%
% Interpolate at a specific query point.

Vq = interp3(X,Y,Z,V,2.35,1.76,0.23,'cubic')

%
% Compare Vq to the value generated by the analytic expression.

V_actual = generatedvalues(2.35,1.76,0.23)

%
%% The interpn Function
% This example shows how the interpn function can be used to interpolate a coarsely sampled function over a finer grid using the 'cubic' method.
% The function interpn performs n-dimensional interpolation on grids that are in the ndgrid format. Its most general form is:
% Vq = interpn(X1,X2,X3,...Xn,V,Y1,Y2,Y3,...,Yn,method)
% where X1,X2,X3,...,Xn are arrays of coordinates that define a grid in ndgrid format, and V is an array containing the values at the grid points. Y1,Y2,Y3,...,Yn are arrays containing the coordinates of the query points at which to interpolate. method optionally specifies any of four interpolation methods: 'nearest', 'linear', 'cubic', or 'spline'.
% The grid points that comprise X1,X2,X3,...Xn must be monotonically increasing and should conform to the ndgrid format.
% Create a set of grid points and corresponding sample values.

[X1,X2] = ndgrid((-5:1:5));
R = sqrt(X1.^2 + X2.^2)+ eps;
V = sin(R)./(R);

%
% Plot the sample values.

mesh(X1,X2,V)
title('Sample Grid');

%
% Create a finer grid for interpolation.

[Y1,Y2] = ndgrid((-5:.5:5));

%
% Interpolate over the finer grid and plot the results.

Vq = interpn(X1,X2,V,Y1,Y2,'cubic');
mesh(Y1,Y2,Vq)
title('Refined Grid');

% interpn has an alternate syntax: Vq = interpn(V,ntimes,method) that allows you to interpolate over a grid that is an integer number of times finer than the sample grid. In the previous code, Y1 and Y2 queried the interpolant over a grid that contained one extra point between each of the samples. The following code demonstrates how you can achieve the same result with ntimes=1.
%
% Interpolate over a finer grid using ntimes=1.

Vq = interpn(V,1,'cubic');

%
% Plot the result.

mesh(Vq)
title('Refined Grid with NTIMES');

% Notice that the plot is scaled differently than in the previous example. This is because we called mesh and passed the values only. The mesh function used a default grid based on the dimensions of Vq. The output values are the same in both cases.
% You can also supply a nonuniform grid of query points. This can be useful if you are interested in querying the interpolant at a higher resolution in one region of the grid. The following code shows how this can be done.
%
% Interpolate over a biased grid.

[Y1, Y2] = ndgrid([-5 -4 -3 -2.5 -2 -1.5 -1.25 -1 -0.75 -0.5 -0.20 0]);
Vq = interpn(X1,X2,V,Y1,Y2,'cubic');

%
% Plot the result.

mesh(Y1,Y2,Vq);
title('Biased Grid');

% Copyright 2012 The MathWorks, Inc.