[SciPy-user] cubic splines in scipy

[Travis Oliphant]

Travis Oliphant wrote:

>Instead, use the fact that the acutal interpolating function is:
>
>y(x) = sum(c_j*beta_3(x/deltax - j))  
>
>where beta_3 is a cubic-spline interpolator function (a symmetric 
>kernel)  which scipy.signal.bspline(x/deltax-j,3)  will evaluate for 
>you, and cj are the (spline coefficients) returned from cspline1d(y).
>
Incidentally, the scipy.signal.bspline function is probably not the most 
intelligently written code (and it's generic). It is a "vectorized" 
function so that it works like a ufunc.  But, you may find something 
that works better. 

The "closed-form" expression for the bspline interpolator is

beta_n(x) = Delta^(n+1){x**n u(x)} / n! 

where u(x) is the step-function that is 0 when u(x)<0  and

Delta^(n+1){f(x)} is equivalent to Delta^n{Delta{f(x)}}

and

Delta{f(x)} == f(x+1/2)-f(x-1/2).

Thus:

beta_3(x) = [(x+2)**3 u(x+2) - 4(x+1)**3 u(x+1) +6 x**3 u(x) - 4(x-1)**3 
u(x-1) + (x-2)**3 u(x-2)]/6

  Breaking this up into intervals we can show:

beta_3(x):

 |x| > 2                  0
1 < |x| <= 2         [(2-|x|)**3]  / 6
0 <= |x| <= 1      [(2-|x|)**3 - 4(1-|x|)**3] / 6

Proving this to yourself from the general formula takes a bit of 
calculation...

So, for any given point we are trying to find the interpolated value for 
we need at most use the spline coefficients strictly within two 
knot-distances away from the current point. 

We should definitely write up an interpolation class using these splines 
and include it with the bspline material, so this kind of calculation 
doesn't have to be refigured all the time.

-Travis