About natural cubic spline

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

About natural cubic spline

José Luis García Pallero
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
[hidden email]
(o<
/ / \
V_/_
Use Debian GNU/Linux and enjoy!
*****************************************

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: About natural cubic spline

Nir Krakauer-3
Octave has csape in the Splines package.
https://octave.sourceforge.io/splines/function/csape.html

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: About natural cubic spline

Sergei Steshenko
In reply to this post by José Luis García Pallero





________________________________
From: José Luis García Pallero <[hidden email]>
To: Help GNU Octave <[hidden email]>; Octave Help <[hidden email]>
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

[hidden email]

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

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave