View License

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

» Watch video

Highlights from
Akima Interpolation

Join the 15-year community celebration.

Play games and win prizes!

» Learn more

4.66667
4.7 | 9 ratings Rate this file 18 Downloads (last 30 days) File Size: 2.29 KB File ID: #1814 Version: 1.0
image thumbnail

Akima Interpolation

by

N. Shamsundar (view profile)

 

10 Jun 2002 (Updated )

Interpolate smooth curve through given points on a plane.

| Watch this File

File Information
Description

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

MATLAB release MATLAB 6.0 (R12)
Other requirements Works with Matlab Versions 3 and upwards.
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (12)
07 Sep 2016 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

17 Mar 2016 Wai Heng

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

Comment only
16 Oct 2009 Chris

Chris (view profile)

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

25 Dec 2006 Kosit Tea

Very Excellent

10 Nov 2006 Jeff Barton

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

05 Nov 2006 mehdi ravanbakhsh

very good

Comment only
14 Oct 2006 Zeke Cummings

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

19 Jan 2006 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

02 Dec 2005 ericka solano

good

Comment only
08 Sep 2004 Heyjin Kim

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

23 Aug 2004 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

07 Oct 2003 jesper dietz

jeo it works nice

Updates
09 Sep 2004

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

Contact us