[SciPy-User] Splines in scipy.signal vs scipy.interpolation
Tony S Yu
tonyyu@MIT....
Wed Jan 20 16:56:33 CST 2010
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). I believe the main difference between the two spline representations (assuming cubic splines with no smoothing) is their boundary conditions (right?). Is there any way to condition the inputs such that I can feed in a spline from `cspline1d` and get a sensible result from `splev`? (Below is an example of what I mean by "conditioning the inputs").
Alternatively, is there another way to get functionality similar to matlab's `spaps` function?
Thanks,
-Tony
#---- Failed attempt to get cspline1d/splev roundtrip
import numpy as np
from scipy import signal, interpolate
x = np.linspace(1, 10, 20)
y = np.cos(x)
tck_interp = interpolate.splrep(x, y)
c_signal = signal.cspline1d(y, 0) # set lambda to zero to eliminate smoothing
# knots and coefficients from splrep have more values at boundaries
t_match = np.hstack(([x[0]]*4, x[2:-2], [x[-1]]*4))
c_match = np.hstack((c_signal, [0]*4))
tck_signal = [t_match, c_match, 3]
y_signal = signal.cspline1d_eval(c_signal, x, dx=x[1]-x[0], x0=x[0])
y_signal_interp = interpolate.splev(x, tck_signal, der=0)
y_interp = interpolate.splev(x, tck_interp, der=0)
print 'spline orders match? ', np.allclose(tck_signal[2], tck_interp[2]) #True
print 'knots match? ', np.allclose(tck_signal[0], tck_interp[0]) #True
print 'spline coefficients match? ', np.allclose(tck_signal[1], tck_interp[1]) #False
print 'y (signal roundtrip) matches? ', np.allclose(y, y_signal) #True
print 'y (interp roundtrip) matches? ', np.allclose(y, y_interp) #True
print 'y (signal in/interp out) matches? ', np.allclose(y, y_signal_interp) #False
More information about the SciPy-User
mailing list