[SciPy-user] cubic splines in scipy

Travis Oliphant oliphant.travis at ieee.org
Sat Jan 28 18:42:03 CST 2006

```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

```