Code covered by the BSD License

# Interpolation Utilities

### Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

cubiconv(x, y, xi)
```function [yi, ypi, yppi] = cubiconv(x, y, xi)

% CUBICONV Cubic convolution interpolation
%    CUBICONV(X,Y,XI) interpolates to find YI, the value of
%    the underlying function Y at the points in the uniform array
%    XI, using cubic convolution interpolation.  X and Y must be
%    vectors of length N.
%
%    [YI,YPI,YPPI] = CUBICONV() also returns the interpolated
%    quartic derivative and cubic second derivative of the underlying
%    function Y at points XI.

% Joe Henning - Fall 2011

% Cubic Convolution Interpolation for Digital Image Processing
% Robert G. Keys
% IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-29, No. 6
% December, 1981

n = length(x);
%h = (x(n) - x(1))/n;

for i = 1:length(xi)
if (xi(i) < x(1) || xi(i) > x(n))
fprintf('??? Bad x input to cubiconv ==> x(1) <= x <= x(n)\n');
yi(i) = NaN;
continue;
end

% Find the right place in the table by means of a bisection.
klo = 1;
khi = n;
while (khi-klo > 1)
k = fix((khi+klo)/2.0);
if (x(k) > xi(i))
khi = k;
else
klo = k;
end
end

h = x(khi) - x(klo);
if (h == 0.0)
fprintf('??? Bad x input to cubiconv ==> x values must be distinct\n');
yi(i) = NaN;
ypi(i) = NaN;
yppi(i) = NaN;
continue;
end

km1 = klo - 1;
kp1 = khi + 1;
if (km1 < 1)
a = 3*y(klo) - 3*y(khi) + y(khi+1);
else
a = y(km1);
end
if (kp1 > n)
d = 3*y(khi) - 3*y(klo) + y(klo-1);
else
d = y(kp1);
end
b = y(klo);
c = y(khi);

% Evaluate cubic polynomial
t = (xi(i) - x(klo))/h;
t2 = t*t;
t3 = t2*t;
h2 = h*h;
c00 = (-t3 + 2*t2 - t)/2.0;
c10 = (3*t3 - 5*t2 + 2)/2.0;
c20 = (-3*t3 + 4*t2 + t)/2.0;
c30 = (t3 - t2)/2.0;

yi(i) = a*c00 + b*c10 + c*c20 + d*c30;

% Differentiate to find the second-order interpolant
c00 = (-3*t2 + 4*t - 1)/2.0;
c10 = (9*t2 - 10*t)/2.0;
c20 = (-9*t2 + 8*t + 1)/2.0;
c30 = (3*t2 - 2*t)/2.0;

ypi(i) = a*c00/h + b*c10/h + c*c20/h + d*c30/h;

% Differentiate to find the first-order interpolant
c00 = (-6*t + 4)/2.0;
c10 = (18*t - 10)/2.0;
c20 = (-18*t + 8)/2.0;
c30 = (6*t - 2)/2.0;

yppi(i) = a*c00/h2 + b*c10/h2 + c*c20/h2 + d*c30/h2;
end
```