File Exchange

image thumbnail

Interpolation Utilities

version 1.34.0.0 (44.9 KB) by Joe Henning
A variety of interpolation utilities

16 Downloads

Updated 12 Jan 2020

View License

Editor's Note: This file was selected as MATLAB Central Pick of the Week

This zip file contains 40 functions related to interpolation. The functions are:
1) analyticint.m performs piecewise analytic interpolation
2) baryinv.m performs barycentric interpolation with inverse distance weighting
3) barylag.m performs barycentric lagrange interpolation
4) bulirschstoer.m performs performs rational Bulirsch-Stoer interpolation
5) cakima.m performs piecewise cubic Akima interpolation
6) cbezier.m performs piecewise cubic Bezier spline interpolation
7) chermite.m performs piecewise cubic Hermite spline interpolation
8) constint performs piecewise constant interpolation
9) cosint.m performs piecewise cosine interpolation
10) cspline.m performs piecewise cubic spline interpolation
11) cubicfarrow.m performs Cubic Farrow interpolation
12) cubiconv.m performs cubic convolution interpolation
13) divdiff.m calculates divided differences
14) expint.m calculates piecewise exponential interpolation
15) floaterhormann.m performs rational interpolation using the Floater-Hormann Method
16) fractint.m performs self-affine fractal interpolation for specific interpolation points
17) hermint.m performs piecewise Hermite interpolation
18) hermitediv.m calculates Hermite divided differences
19) hilbint.m performs piecewise Hilbert transform interpolation
20) interpdct.m performs interpolation using the DCT method
21) lagint.m performs piecewise Lagrange interpolation
22) linint.m performs piecewise linear interpolation
23) minmaxpoly.m calculates the optimal polynomial via the Remez algorithm
24) mqspline.m performs piecewise monotone quadratic spline interpolation
25) neville.m performs interpolation using Neville's Method
26) newtint.m performs interpolation of equally-spaced points
27) polycof.m performs non-Vandermonde polynomial interpolation coefficients
28) qhermite.m performs piecewise quintic Hermite interpolation
29) qspline.m performs piecewise quadratic spline interpolation
30) ratint.m performs simple rational polynomial interpolation
31) rchermite.m performs piecewise rational cubic Hermite spline interpolation
32) said.m performs piecewise Said interpolation
33) schwerner.m performs rational interpolation using the Schneider-Werner Method
34) shermite.m performs piecewise septic Hermite interpolation
35) sincdint.m performs piecewise discrete sinc interpolation
36) sincint.m performs piecewise sinc interpolation
37) steffen.m performs monotonic Steffen interpolation
38) stineman.m performs monotonic Stineman interpolation
39) thiele.m performs rational Thiele interpolation
40) trigint.m performs piecewise trigonometric interpolation

Cite As

Joe Henning (2020). Interpolation Utilities (https://www.mathworks.com/matlabcentral/fileexchange/36800-interpolation-utilities), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (13)

Once more in addition to stineman.m issue reported below:

Both, dylo and dyhi, can become 0.

The implementation will then do this:
prod = dylo*dyhi;
if (abs(prod) < eps)
yi(i) = 0;

The proposed fix was:
if (prod >= 0)
yi(i) = y0 + prod/(dylo + dyhi);
else
...

But this will result in a division by 0 if dylo == 0 and dyhi == 0.

I suggest to interpolate linear in this case:

if dylo == 0 && dyhi == 0
yi(i)=interp1(x,y,xi(i));
else
% test product
prod = dylo*dyhi;
if (prod > 0)
yi(i) = y0 + prod/(dylo + dyhi);
else % prod < 0
yi(i) = y0 + prod*(xi(i) - x(klo) + xi(i) - x(khi))/((dylo - dyhi)*(x(khi) - x(klo)));
end
end

I detected an issue with stineman.m.

% test product
prod = dylo*dyhi;
if (abs(prod) < eps)
yi(i) = 0;
elseif (prod > 0)
yi(i) = y0 + prod/(dylo + dyhi);
else % prod < 0
yi(i) = y0 + prod*(xi(i) - x(klo) + xi(i) - x(khi))/((dylo - dyhi)*(x(khi) - x(klo)));
end

'prod' can become 0 for a lot of samples in 'xi' depending on the slope of the nodes and the resulting values for 'dylo' or 'dyhi'.

In this case, yi(i) becomes 0 until the next node is being evaluated, hence the interpolation fails in these cases.

Try with

x = [0;10;20;30;40;50;60;70;80;90;100];
y = [0;2;5;8;13;13;12;10;5;1;0];
xi = [0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100]
yi = stineman(x,y,xi);
This results in
yi = [0;0.16;0.33;0.51;0.7;0.89;1.09;1.3;1.53;1.76;2;0;0;0;0;0;0;0;0;4.70;5;0;0;0;0;0;0;0;0;0;8;8.42;8.88;9.39;9.93;10.5;11.08;11.66;12.2;12.67;13;13.15;13.21;13.23;13.22;13.2;13.17;13.14;13.09;13.05;13;12.95;12.88;12.8;12.72;12.62;12.52;12.4;12.28;12.14;12;11.85;11.69;11.53;11.36;11.18;11;10.79;10.56;10.3;10;9.6;9.11;8.58;8.04;7.50;6.97;6.45;5.94;5.46;5;4.55;4.11;3.67;3.24;2.81;2.4;2;1.63;1.28;1;0.79;0.63;0.49;0.380;0.29;0.22;0.150;0.09;0.04;0]

Note the row of zeros for samples 12 to 19 and 22 to 29.

The following change seems to work (at least for the samples I evaluated);

if (prod >= 0)
yi(i) = y0 + prod/(dylo + dyhi);
else % prod < 0
yi(i) = y0 + prod*(xi(i) - x(klo) + xi(i) - x(khi))/((dylo - dyhi)*(x(khi) - x(klo)));
end

This results in
yi = [;0.16;0.33;0.51;0.7;0.89;1.09;1.3;1.53;1.76;2;2.3;2.6;2.9;3.2;3.5;3.8;4.1;4.4;4.7;5;5.3;5.6;5.9;6.2;6.5;6.8;7.1;7.4;7.7;8;8.42;8.88;9.39;9.93;10.5;11.08;11.66;12.2;12.67;13;13.15;13.21;13.23;13.22;13.2;13.17;13.14;13.09;13.05;13;12.95;12.88;12.8;12.72;12.62;12.52;12.4;12.28;12.14;12;11.85;11.69;11.53;11.36;11.18;11;10.79;10.56;10.3;1;9.6;9.11;8.58;8.04;7.5;6.97;6.45;5.94;5.46;5;4.55;4.11;3.67;3.24;2.81;2.4;2;1.63;1.28;1;0.79;0.63;0.49;0.38;0.29;0.22;0.15;0.09;0.04;0]

Looks OK to me.

I detected an issue with cakima.m today.

If I interpolate my test data (see below) then I get a wrong cakima result at the second value of the interpolated data:
yi(2) = 1.42108547152020e-14
I know that it is wrong because it is almost 0.0, which should not be the case.
The interpolateted testdata with two other interpolation methods (Matlab spline.m and another akima implementation) is significantly greater than 0 and similar to each other:
"correct" spline result: yi(2) = 0.00271325527955031
"correct" akima result: yi(2) = 0.00271252847004375;
Both "correct" results match each other very good.

Testdata:
% x = [0,0.383865456343585,0.875278890624145,1.36427252920419,1.85081835187610,2.33488834388682,2.81645449615647,3.30875890213734,3.79835845825298,4.28522277619577,4.76932147508651,5.25062418178334,5.74199281280865,6.23034736093939,6.71565495089055,7.19788271754666,7.68956376377395,8.17793064649527,8.66294787356091,9.14457996625234,9.63500799909241,10.1217995214661,10.6049163032000,11.0962573240730,11.5836562937140,12.0670721037405,12.5581055251765,13.0448718631453,13.5273270056942,14.0167567534549,14.5015740863952,14.9928424018129,15.4791812466965,15.9605400950058,16.4478769606337,16.9316225458693,17.4226260008344,17.9100601265062,18.3938857474978,18.8844615704438,19.3711877214937,19.8540225326813,20.3430707442944,20.8279725627409,21.3186536879540,21.8049207561054,22.2867268308059,22.7737185539466,23.2559668189102,23.7429218533379,24.2248375555182,24.7109596658190,25.2009256290159,25.6852762204520,26.1729296222031,26.6546312329393,27.1390713996938,27.6258431737713,28.1145263715792,28.5962973164040,29.0799101072968,29.5682134217497,30.0532118115966,30.5428786780211,31.0287993707382,31.5108086957148,31.9966717684461,32.4860039798981,32.9706233800920,33.4580453058658,33.9402246355093,34.4244948230778,34.9103906038430,35.3974242185693,35.8850848300977,36.3728379566002,36.8601249269250,37.3463623646877,37.8309417093219,38.3198612775001,38.8056842031725,39.2919888135210,39.7803723764687,40.2640775113193,40.7491256903732,41.2350618585002,41.7214094438180,42.2076699872172,42.6933228135684,43.1778247542674,43.6606099329746,44.1467591268708,44.6297817364699,45.1144803480354,45.5998309222389,46.0847504576474,46.5683936251959,47.0551384312226,47.5412722989731,48.0262287069314,48.5142636923878,48.9997706685131,49.4867899964211,49.9743934278086,50.4615990620718,50.9473717752685,51.4349240110933,51.9185927438624,52.4053029508257,52.8921013292307,53.3760150810979,53.8603654349716,54.3479873565715,54.8339068283263,55.3206875437889,55.8069397540581,56.2912005379835,56.7753013193231,57.2605826862390,57.7456811721952,58.2322124132459,58.7189567020615,59.2045520712042,59.6905098438515,60.1751399480906,60.6594670572142,61.1441423530091,61.6285900861004,62.1150404667240,62.6017175524324,63.0866803791180,63.5727020736038,64.0573974194261,64.5432355441556,65.0292287444836,65.5149361691336,66.0009613156384,66.4856679831040,66.9714192369191,67.4578484318308,67.9441154814398,68.4298100507477,68.9155633505963,69.4018524878555,69.8879903378167,70.3745555031874,70.8605295359052,71.3459796590441,71.8324501349799,72.3176294031312,72.8031576200045,73.2891547482618,73.7742443474892,74.2597120932594,74.7454616011734,75.2311682218484,75.7164490076203,76.2015216512112,76.6875222596191,77.1734718945258,77.6595457975884,78.1456639825593,78.6314657186151,79.1167854553498,79.6025120178138,80.0880722626421,80.5739079410733,81.0598064095473,81.5458964826904,82.0316497431547,82.5172967286252,83.0027807026799,83.4888527758970,83.9744830846248,84.4601077835561,84.9456785380408,85.4308846729043,85.9162000130151,86.4016236731773,86.8872279239453,87.3726208872652,87.8580816871142,88.3436974798867,88.8297349391148,89.3155639946326,89.8014454598648,90.2873243185875,90.7731678201164,91.2589159039694,91.7447477870474,92.2301295031583,92.7155594666181,93.2009292992969,93.6864134812118,94.1719035880400,94.6573933152480,95.1431768116826,95.6286463995482,96.1145543195491,96.6000447049951,97.0858297556004,97.5714143233177,98.0571418994533,98.5429323444118,99.0286507496608,99.5143536715751,100];
% y = [0,0.0427252613107498,0.0976577401388568,0.152590218966964,0.207522697795071,0.262455176623178,0.317387655451284,0.373846036469061,0.430304417486838,0.486762798504614,0.543221179522391,0.599679560540167,0.657663843747613,0.715648126955060,0.773632410162506,0.831616693369952,0.891126878767068,0.950637064164184,1.01014724956130,1.06965743495842,1.13069352254520,1.19172961013199,1.25276569771877,1.31532768749523,1.37788967727168,1.44045166704814,1.50453955901426,1.56862745098039,1.63271534294651,1.69832913710231,1.76394293125810,1.83108262760356,1.89822232394903,1.96536202029449,2.03402761882964,2.10269321736477,2.17288471808958,2.24307621881438,2.31326771953918,2.38498512245366,2.45670252536813,2.52841992828260,2.60166323338675,2.67490653849089,2.74967574578470,2.82444495307851,2.89921416037232,2.97550926985581,3.05180437933929,3.12962539101244,3.20744640268559,3.28679331654841,3.36766613260090,3.44853894865339,3.53093766689555,3.61333638513771,3.69726100556954,3.78271152819104,3.86968795300221,3.95666437781338,4.04516670481422,4.13519493400473,4.22522316319524,4.31677729457542,4.40833142595560,4.49988555733577,4.59296559090562,4.68757152666514,4.78217746242466,4.87830930037384,4.97444113832303,5.07209887846189,5.17128252079041,5.27199206530861,5.37422751201648,5.47798886091401,5.58327611200122,5.69008926527809,5.79842832074463,5.90981918059052,6.02273594262607,6.13717860685131,6.25314717326620,6.36911573968109,6.48661020828565,6.60563057907989,6.72617685206379,6.84824902723736,6.97184710460060,7.09697108415351,7.22362096589609,7.35332265201801,7.48455024032960,7.61882963302053,7.75616083009079,7.89654383154040,8.03997863736934,8.18646524757763,8.33447775997558,8.48401617456321,8.63660639353017,8.79072251468681,8.94789044022278,9.10811017013809,9.27138170443274,9.43770504310673,9.60860608834973,9.78255893797208,9.96261539635310,10.1472495613031,10.3334096284428,10.5226214999619,10.7164110780499,10.9132524605173,11.1146715495537,11.3206683451591,11.5312428473335,11.7479209582666,11.9722285801480,12.2026398107881,12.4376287479973,12.6771953917754,12.9213397421225,13.1715877012283,13.4279392690928,13.6919203479057,13.9650568398566,14.2458228427558,14.5326924544137,14.8256656748302,15.1247425040055,15.4329747463188,15.7503624017700,16.0799572747387,16.4171816586557,16.7605096513313,17.1114671549554,17.4700541695277,17.8393224994278,18.2192721446555,18.6053253986420,18.9974822613871,19.3972686350805,19.8062104219120,20.2243076218814,20.6500343327993,21.0833905546654,21.5259021896696,21.9806210421912,22.4444953078508,22.9175249866484,23.4012359807736,23.8956282902266,24.4007019150072,24.9118791485466,25.4291599908446,25.9525444419013,26.4805065995270,27.0130464637217,27.5516899366751,28.0994888227665,28.6564431219959,29.2225528343633,29.7993438620584,30.3898680094606,30.9956511787594,31.6182192721447,32.2575722896162,32.9106584267948,33.5759517814908,34.2534523537041,34.9416342412451,35.6435492484932,36.3607232776379,37.0962081330587,37.8515297169451,38.6266880292973,39.4216830701152,40.2380407415885,41.0727092393378,41.9256885633631,42.7969787136645,43.6896314946212,44.6112764171817,45.5497062638285,46.5003433279927,47.4631876096742,48.4428168154421,49.4514381628138,50.4814221408408,51.5419241626612,52.6436255436027,53.8078889143206,55.0301365682460,56.3073167009995,57.6287479972534,59.0218966964218,60.5142290379187,62.1957732509346,64.1901274128328,66.6971847104601,69.9229419394217,73.1578545815213,76.0234988937209,80.0259403372244,93.2097352559701,100];
% xi = 0.0:100.0/4095:100.0;

JOSE ROMAY

HIT

In thiele interpolation equation what is mean by X,Y,Xi,
i want to remove noise from image by using newton thiele interpolation model. please suggest me the ideas... thanks in advance.....

Andreas

Hi,

One more issue with mqspline. In code line 84, it should be
yp(1) = (3*d(1) - yp(2))/2;
instead of
yp(1) = -(3*d(1) - yp(2))/2;

The additional minus is a typo in Judd (1998), "Numerical Methods in Economis". I asked him and he confirmed. In the original paper (Schumaker (1983), "On Shape Preserving Quadratic Spline Interpolation") Kenneth Judd is referring to, there is no minus.

Obviously, this matters for the interpolation at the left end.

Kind regards,
Andreas

Andreas

Hi,

Thanks a lot for the codes. I'm using mqspline.m.

I think that there is an error in both lines 126 and 129. It should be

eps = x(klo) + h*(b-del)/(b-a);
eps = x(khi) + h*(a-del)/(b-a);

instead of

eps = x(klo) + 2*h*(b-del)/(b-a);
eps = x(khi) + 2*h*(a-del)/(b-a);

See Schumaker (1983), "On Shape Preserving Quadratic Spline Interpolation", or the exposition in Judd (1998), "Numerical Methods in Economis", for the formulae.

Kind regards,
Andreas

Grzesiek

Sorry for my last comment. Everything is alright. My mistake. Data were not equally spaced.

Very good utilities. Thanks

Grzesiek

in newtint.m file there is an error in Everett's method (line 163) and Stirling's method (line 136).

Error msg: newtint: A(I,J): row index out of bounds; value 5 out of bound 4

I tried to execute the following command:

y = [308.6 362.6 423.3 491.4];
x = [0.055389 0.047485 0.040914 0.035413];
c = 4; %c=3
newtint(x,y,[0.04],c)

Updates

1.34.0.0

removed unnecessary variable printout in schwerner.m

1.33.0.0

added analyticint.m, hilbint.m, and references to help sections
fixed baryinv.m, bulirschstoer.m, thiele.m

1.32.0.0

bug fix for stineman.m

1.31.0.0

added cubicfarrow.m

1.30.0.0

Added barylag.m, constint.m, cspline.m, linint.m, minmaxpoly.m, polycof.m
Removed safif.m

1.29.0.0

Updated bulirschstoer.m

1.28.0.0

Corrected mqsline.m and added bulirschstoer.m, hermitediv.m, qspline.m

1.27.0.0

Added ratint.m and thiele.m for rational polynomial interplation

1.26.0.0

Added steffen.m and stineman.m for monotonic interpolation

1.25.0.0

Changed kernel tolerance in schwerner.m

1.24.0.0

Added boundary conditions to helper function in newtint.m

1.22.0.0

Added fractint.m and safif.m and added a fractional delay filter to interpdct.m

1.21.0.0

Modified window condition in sincint.m

1.20.0.0

Added expint.m; added second derivative calculations to lagint.m and hermint.m

1.19.0.0

added expint.m

1.18.0.0

added mqspline.m and harmonic mean monotone options

1.16.0.0

Updated output arguments when error conditions occur

1.15.0.0

updated monotone condition in chermite.m

1.14.0.0

updated trigint.m to handle an even number of datapoints

1.13.0.0

updated trigint.m

1.12.0.0

added comment to said.m

1.11.0.0

Added said.m for piecewise Said interpolation

1.9.0.0

updated trigint.m

1.8.0.0

updated epsilon computations

1.7.0.0

Added finite difference function to rchermite.m

1.6.0.0

Updated sincint.m

1.2.0.0

Modified sincdint.m and sincint.m and added trigint.m

1.1.0.0

Changed the finite difference calculation in chermite.m
Added a septic hermite interpolator, shermite.m

MATLAB Release Compatibility
Created with R14SP3
Compatible with any release
Platform Compatibility
Windows macOS Linux
Categories