File Exchange

image thumbnail

Akima Interpolation

version 1.0 (2.29 KB) by

Interpolate smooth curve through given points on a plane.

22 Downloads

Updated

View License

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

Comments and Ratings (14)

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).

chunlin fu

geiliable

Tsinghua

I modified it for input with descending order. Now it can accept all orders

function yi=akima_nosequence(x,y,xi)
%
% 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.
%
rowcolumn=size(xi);
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

dx=diff(x);
Sign=0;

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
Sign=1;
x=flipud(x);
y=flipud(y);
xi=flipud(xi);
end

if any(xi<x(1)) | any(xi > x(n))
  warning('All interpolation points xi must lie between x(1) and x(n)')
end

if Sign==1
dx=diff(x);
end

m=diff(y)./dx;
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);
b(id)=(f1(id).*m1(id+1)+f2(id).*m1(id+2))./f12(id);
c=(3*m-2*b(1:n-1)-b(2:n))./dx;
d=(b(1:n-1)+b(2:n)-2*m)./dx.^2;

% Loop replaced by vector ops following tip from johannes.korsawe@volkswagen.de
% 1/19/2006
%
[ncnt,bin]=histc(xi,x);
bin=min(bin,n-1);
bb=bin(1:length(xi));
wj=xi-x(bb);
yi=((wj.*d(bb) +c(bb)).*wj+b(bb)).*wj+y(bb);

if Sign==1
    yi=flipud(yi);
end
if rowcolumn(2)>rowcolumn(1)% if xi is a row, returns a row
yi=yi';
end

Wai Heng

Hi, has anyone modified the Akima spline for 2D interpolation? My purpose is to use it for spatial interpolation. Thanks

Chris

Chris (view profile)

Works well, It easily be made compatible with the ppval function too.

Kosit Tea

Very Excellent

Jeff Barton

Excellent! Excellent! Excellent! I was looking for this exact capability, and the routine performed beautifully. Thank you!

mehdi ravanbakhsh

very good

Zeke Cummings

With the Korsawe comment it's a 5 star function. More useful than Akima's orginal fortran algorithm

Johannes Korsawe

Great work! Besides, the routine can be accelerated a lot by using the lines

[nerd,bin]=histc(xi,x);bin=min(bin,n-1);
bb=bin(1:length(xi))';wj=xi-x(bb);yi=((wj.*d(bb) +c(bb)).*wj+b(bb)).*wj+y(bb);

instead of the loop at the end. Speed-Up of 500 for 50.000 interpolated points.

Best,
Johannes Korsawe

ericka solano

good

Heyjin Kim

It works great for interpolation of ocean temperature and salinity profiles. Thanks!

Catherine Martel

Hello,

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....

Cheers,

Catherine Martel

jesper dietz

jeo it works nice

Updates

One line corrected, to correctly handle special case where input data fall on straight line (9/2004)

MATLAB Release
MATLAB 6.0 (R12)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video