Implementation of Akima's univariate interpolation method (Journal of the ACM, Vol. 17, No. 4, October 1970, pages 589-602).
N. Shamsundar, University of Houston
This is a great implementation of the Akima 1970 interpolation method (Akima-70). This gives less ringing and overshooting than the FFT interpolations, or natural, cubic, and not-a-knot spline algorithms, while also not introducing the broadening of apodized FFT interpolations or other convolution based interpolations.
Do take caution when data rises or falls quickly to a baseline, this algorithm will include some ringing at the top and bottom of steep segments. There is an improved Akima algorithm (see 1986 US Dept. of Commerce NTIA report 86-208 by Hiroshi Akima). The improved Akima algorithm (Akima-86) is much less sensitive to those types of conditions, while also having the positive features of the Akima 1970 algorithm. It's also guaranteed to fit with a cubic polynomial accuracy, whereas the Akima-70 algorithm only guarantees fits to a second order polynomial (despite the cubic polynomial interpolation of the first derivatives).
I modified it for input with descending order. Now it can accept all orders
% Usage: yi=akima_nosequence(x,y,xi)
% Given vectors x and y (of the same length)
% and the array xi at which to interpolate,
% fits piecewise cubic polynomials and returns
% the interpolated values yi at xi.
% Ref. : Hiroshi Akima, Journal of the ACM, Vol. 17, No. 4, October 1970,
% pages 589-602.
% Programmer: N. Shamsundar, University of Houston, 6/2002
% Correction to lines 32-33, 9/2004, motivated by Gilford Ward,
% to make routine work correctly for linear data.
% Notes: Use only for precise data, as the fitted curve passes through the
% given points exactly. This routine is useful for plotting a pleasingly
% smooth curve through a few given points for purposes of plotting.
x=x(:); y=y(:); xi=xi(:); n=length(x);% change row to column
if n~=length(y), error('input x and y arrays must be of same length'), end
if any(dx <= 0)
% error('input x-array must be in strictly ascending order'),
% if the sequence of x array is descending, then fliplr(x) fliplr(y)
% it is different from akima.m, it needs x array to be ascending
if any(xi<x(1)) | any(xi > x(n))
warning('All interpolation points xi must lie between x(1) and x(n)')
mm=2*m(1)-m(2); mmm=2*mm-m(1); % augment at left
mp=2*m(n-1)-m(n-2); mpp=2*mp-m(n-1); % augment at right
m1=[mmm; mm; m; mp; mpp]; % slopes
dm=abs(diff(m1)); f1=dm(3:n+2); f2=dm(1:n); f12=f1+f2;
id=find(f12 > 1e-8*max(f12)); b=m1(2:n+1);
% Loop replaced by vector ops following tip from firstname.lastname@example.org
if rowcolumn(2)>rowcolumn(1)% if xi is a row, returns a row
Hi, has anyone modified the Akima spline for 2D interpolation? My purpose is to use it for spatial interpolation. Thanks
Works well, It easily be made compatible with the ppval function too.
Excellent! Excellent! Excellent! I was looking for this exact capability, and the routine performed beautifully. Thank you!
With the Korsawe comment it's a 5 star function. More useful than Akima's orginal fortran algorithm
Great work! Besides, the routine can be accelerated a lot by using the lines
instead of the loop at the end. Speed-Up of 500 for 50.000 interpolated points.
It works great for interpolation of ocean temperature and salinity profiles. Thanks!
I'm trying to get in touch with Jesper Dietz. The Jesper I'm looking for was crazy about sailing but sold his boat in order to go to Switzerland for a year, learn French, ski and have a fabulous time with new friends. Unfortunately, we've lost touch over the years and I'm now trying to reconnect. If you happen to be this Jesper or know him, please write back....
jeo it works nice
One line corrected, to correctly handle special case where input data fall on straight line (9/2004)