[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
More information about the SciPy-user
mailing list