Re: About natural cubic spline

[Sergei Steshenko]

________________________________
From: José Luis García Pallero <address@hidden>
To: Help GNU Octave <address@hidden>; Octave Help <address@hidden> 
Sent: Sunday, November 26, 2017 7:04 PM
Subject: About natural cubic spline



Hello:


Reading the help text of spline() I can see that the computed spline

is the not-a-knot versión (continuous third derivatives at the second

and the penultimate points), but there is not an option to determine

the natural spline (second derivatives equal to zero at the first and

end points). This is the case in Octave and Matlab, but the last one

have the function csape() in the Curve Fitting toolbox which provides

the option to compute the natural spline.


Why spline() has not such an option in Octave? Is only for Matlab

compatibility? Has Octave other function to compute the natural

spline?


Thanks


-- 

*****************************************

José Luis García Pallero

address@hidden

(o<

/ / \

V_/_

Use Debian GNU/Linux and enjoy!

*****************************************


_______________________________________________



My reply will be somewhat off topic, but maybe useful.

At the moment I'm playing with interpolation polynomials - I need to 
interpolate solutions of differential equations.

Performing web search I stumbled upon "Polynomial Interpolators for 
High-Quality Resampling of Oversampled Audio" - 
http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf . The article has 
formulas for a number of "classical" polynomials, including splines, Lagrange, 
Hermit, as well as for polynomials created by the author specifically for the 
needs audio reasmpling AFTER a sync interpolator. It also has "C" code snippets 
which are easy to convert into real "C" code, but be careful - negative indexes 
are used.


I've tested quality of interpolation using 2nd order LTI as source of "golden" 
data (per 
http://academic.csuohio.edu/richter_h/courses/mce441/mce441_7n.pdf ). I used 
natural omega equal to 1 and

___zetas[] = {0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10};
.

X (i.e. time) range was 0 .. 4, and the step for "golden" data was 1/128. 
Interpolation was performed at 1/8 of the golden data step.

The best data:

p_interp_p4_o3_lagrange_x_form rms_delta=3.98165e-07 max_abs_delta=3.60401e-05
p_interp_p4_o3_lagrange_z_form rms_delta=3.98165e-07 max_abs_delta=3.60401e-05
p_interp_p6_o5_lagrange_x_form rms_delta=7.96166e-08 max_abs_delta=7.721e-06
p_interp_p6_o5_lagrange_z_form rms_delta=7.96166e-08 max_abs_delta=7.721e-06
p_interp_p4_o3_hermite_x_form rms_delta=7.49479e-07 max_abs_delta=5.31025e-05
p_interp_p4_o3_hermite_z_form rms_delta=7.49479e-07 max_abs_delta=5.31025e-05
p_interp_p6_o3_hermite_x_form rms_delta=6.83208e-08 max_abs_delta=8.20603e-06
p_interp_p6_o3_hermite_z_form rms_delta=6.83208e-08 max_abs_delta=8.20603e-06
p_interp_p6_o5_hermite_x_form rms_delta=1.017e-07 max_abs_delta=1.15619e-05
p_interp_p6_o5_hermite_z_form rms_delta=1.017e-07 max_abs_delta=1.15619e-05
p_interp_p4_o5_osculating_x_form rms_delta=1.03598e-06 max_abs_delta=6.74535e-05
p_interp_p4_o5_osculating_z_form rms_delta=1.03598e-06 max_abs_delta=6.74535e-05
p_interp_p6_o5_osculating_x_form rms_delta=1.18253e-07 max_abs_delta=1.40449e-05
p_interp_p6_o5_osculating_z_form rms_delta=1.18253e-07 max_abs_delta=1.40449e-05
.

In the above in pN (e.g. p6) N means number of points used to construct the 
interpolation polynomial (for p6 it means 6 points were used) and oM (e.g. o5) 
means order of interpolation polynomial (for o5 order is 5).

rms_delta means root mean square delta, where delta is the difference between 
"golden" value and interpolated value.

max_abs_delta means max absolute delta.

So, of the above delta I quite like

p_interp_p6_o3_hermite_x_form rms_delta=6.83208e-08 max_abs_delta=8.20603e-06

because of the smallest of all rms_delta.

Splines did not work so well for me:


p_interp_p4_o3_bspline_x_form rms_delta=0.000152163 max_abs_delta=0.000434817
p_interp_p4_o3_bspline_z_form rms_delta=0.000152163 max_abs_delta=0.000434817
p_interp_p6_o5_bspline_x_form rms_delta=0.000228235 max_abs_delta=0.000674847
.

I didn't change anything in the formulas from the article - just converted the 
code into full fledged "C" and tested it. Maybe I've introduced some bugs, but 
unlikely - because the code is very simple and I checked results at some 
intermediate points.

My point is that maybe splines is not the best solutions. Another point is that 
thanks to presence of "C" code snippets it's easy to transform the code into 
Octave one.

--Sergei.