[SciPy-User] Splines in scipy.signal vs scipy.interpolation

[Tony S Yu]

On Jan 28, 2010, at 9:20 AM, denis wrote:

> 
> 
> On Jan 20, 11:56 pm, Tony S Yu <ton... at MIT.EDU> wrote:
>> I'm having trouble making splines from scipy.signal work with those in scipy.interpolation.
>> 
>> Both packages have functions for creating (`signal.cspline1d`/`interpolate.splrep`) and evaluating (`signal.cspline1d_eval`/`interpolate.splev`) splines. There are, of course, huge differences between these functions, which is why I'm trying to get them to talk to each other.
>> 
>> In particular, I'd like to create a smoothing spline using `cspline1d` (which allows easier smoothing) and evaluate using `splev` (which allows me to get derivatives of the spline).
> 
> Tony,
>  bouncing between two murky packages doesn't sound as though it'll
> converge ...

Agreed. This was more of a naive attempt to try and get the results that I wanted.

> interpolate though has both smoothing and derivs --
>        interpolator = interpolate.UnivariateSpline( x, y, k=3, s=s )
> # s=0 interpolates
>        yy = interpolator( xx )
>        y1 = interpolator( xx, 1 )  # deriv

You're right. When I originally read the docs for splrep, I had it in my head that the splines in scipy.interpolation didn't provide the "right" type of smoothing (don't ask me what "right" means---I have no idea). After taking some time to understand the interpolation module, I realize it does what I want. Thanks, Denis!

> Just curious, are your real knots uniform, how many ?

I'm actually converting some matlab code which tries to find the optimal smoothed-spline, so the number of knots really depends on the data (the data is uniformly-spaced with about a 1000 points, but the knots depend on the smoothing---if i understand your question correctly).

BTW, if anyone else ever needs to compare results from matlab's `spaps` with smoothed splines using splrep or UnivariateSpline: Note there's a big difference in the **default** error calculations in matlab and scipy. 

If you want to match matlab's error calculation, you need to pass in weights ("w") that match the weights used for the trapezoidal rule. Also there's a subtle, but important difference between the error equations: in matlab "w" is outside the square of the differences; in scipy "w" is inside the square of the differences.

In short, to match matlab's error calculation, you need to pass "w" to splrep or UnivariateSpline, where w = np.sqrt(trapz_weights(x))

>>> def trapz_weights(x):
>>>     dx = np.diff(x)
>>>     w = np.empty(x.shape)
>>>     w[1:-1] = (dx[1:] + dx[:-1])/2.
>>>     w[0] = dx[0] / 2.
>>>     w[-1] = dx[-1] / 2.
>>>     return w

Unfortunately, the splines produced by matlab and scipy don't really match (not sure why---different smoothing algorithms?), but at least their errors are the same.

Cheers,
-T

> 
> See also http://projects.scipy.org/scipy/ticket/864
> "The documentation for class scipy.interpolate.UnivariateSpline? is
> misleading, and maybe completely wrong.
> UnivariateSpline? behaves in ways that are unpredictable ...
> (Fitpack is just a big dense package => big dense doc.)
> 
> cheers
>  -- denis
> 
> 
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user