在周期性边界条件的情况下,如何在 scipy.interpolate.splprep 中设置手动结?
How to set manual knots in scipy.interpolate.splprep in case of periodic boundary conditions?
或
为什么 splprep
不能使用它自己的结?
我想弄清楚如何在 scipy.interpolate.splprep
中设置结。
在非周期性情况下我成功了,即我可以重现 this SE example.
不过,在周期性边界条件 (PBC) 的情况下,我遇到了问题。这里 scipy.interpolate.splprep
甚至不能使用自己的结,如本例所示:
import numpy as np
from scipy.interpolate import splev, splrep
import scipy
print "my scipy version: ", scipy.__version__
srate = 192000.
freq = 1200.
timeList = [ i / srate for i in range( int( srate / freq + 1 ) ) ]
signal = np.fromiter( ( np.sin( 2 * np.pi * freq * t ) for t in timeList ), np.float )
spl = splrep( timeList, signal, per=1 )
knots = spl[0]
spl = splrep( timeList, signal, t=knots, per=1 )
给出:
my scipy version: 1.0.0
Traceback (most recent call last):
File "splrevtest.py", line 15, in <module>
spl = splrep( timeList, signal, t=knots, per=1 )
File "/somepath/python2.7/site-packages/scipy-1.0.0-py2.7-linux-/python2.7/site-packages/scipy-1.0.0-py2.7-linux-x86_64.egg/scipy/interpolate/fitpack.py", line 289, in splrep
res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
File "/somepath/python2.7/site-packages/scipy-1.0.0-py2.7-linux-x86_64.egg/scipy/interpolate/_fitpack_impl.py", line 514, in splrep
raise _iermess[ier][1](_iermess[ier][0])
ValueError: Error on input data
此外,如果您将第一个和最后一个结绘制成:
print timeList[-1]
print knots[:4]
print knots[-4:]
你得到
>> 0.000833333333333
>> [-1.56250000e-05 -1.04166667e-05 -5.20833333e-06 0.00000000e+00]
>> [0.00083333 0.00083854 0.00084375 0.00084896]
表示实际数据范围外有六个点,前三个,后三个。这可能没问题,因为我们有 PBC,结与数据范围内的结相同,但无论如何都很奇怪。 (顺便说一句,这个例子在设置 per=0
时甚至不起作用。)此外,如果我手动设置点,设置数据范围之外的点会失败。即使使用 per=1
,如果点正好在数据范围内,它也能工作。比如:
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import splev, splrep
import scipy
x = np.linspace( 0, 1, 120 )
timeList = np.linspace( 0, 1, 15 )
signal = np.fromiter( ( np.sin( 2 * np.pi * t +.3 ) for t in timeList ), np.float )
signal[ -1 ] -= .15
myKnots=np.linspace( .05, .95, 8 )
spl = splrep( timeList, signal, t=myKnots, per=1 )
fit = splev(x,spl)
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( timeList, signal, marker='o' )
ax.plot( x, fit , 'r' )
for i in myKnots:
ax.axvline( i )
plt.show()
提供:
我们还可以看到最后一点在 PBC 中被忽略了。
那么 scipy.interpolate.splprep
在这里实际做了什么,如果 per=1
为什么它不接受类似的手动设置结的结结构。这是一个错误吗?
我希望最后一个问题的答案是 'no',否则我在这里问这个问题就错了。
为了使 splrep
使用您自己的结,您应该更改:
spl = splrep(timeList, signal, t=knots, per=1)
使用内部结:
spl = splrep(timeList, signal, t=knots[4:-4], per=1)
这个界面不是很直观,但它出现在the documentation中(尽管在task
参数下而不是在t
下)(我的重点):
If task=-1 find the weighted least square spline for a given set of knots, t. These should be interior knots as knots on the ends will be added automatically.
那里也记录了额外的结。这与结数与系数数之间的一般联系是一致的:number_of_knots = number_of_coefficients + degree + 1
.
PBC 构造中的最后一个点被忽略(如图所示),因为周期定义要求周期中最后一个点的值等于第一个点。这是在 per
参数下记录的 "Values of y[m-1] and w[m-1] are not used."
(同样,我猜这不是一个非常直观的界面,但这可能是 Fortran 代码的实现方式)。
或
为什么 splprep
不能使用它自己的结?
我想弄清楚如何在 scipy.interpolate.splprep
中设置结。
在非周期性情况下我成功了,即我可以重现 this SE example.
不过,在周期性边界条件 (PBC) 的情况下,我遇到了问题。这里 scipy.interpolate.splprep
甚至不能使用自己的结,如本例所示:
import numpy as np
from scipy.interpolate import splev, splrep
import scipy
print "my scipy version: ", scipy.__version__
srate = 192000.
freq = 1200.
timeList = [ i / srate for i in range( int( srate / freq + 1 ) ) ]
signal = np.fromiter( ( np.sin( 2 * np.pi * freq * t ) for t in timeList ), np.float )
spl = splrep( timeList, signal, per=1 )
knots = spl[0]
spl = splrep( timeList, signal, t=knots, per=1 )
给出:
my scipy version: 1.0.0
Traceback (most recent call last):
File "splrevtest.py", line 15, in <module>
spl = splrep( timeList, signal, t=knots, per=1 )
File "/somepath/python2.7/site-packages/scipy-1.0.0-py2.7-linux-/python2.7/site-packages/scipy-1.0.0-py2.7-linux-x86_64.egg/scipy/interpolate/fitpack.py", line 289, in splrep
res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
File "/somepath/python2.7/site-packages/scipy-1.0.0-py2.7-linux-x86_64.egg/scipy/interpolate/_fitpack_impl.py", line 514, in splrep
raise _iermess[ier][1](_iermess[ier][0])
ValueError: Error on input data
此外,如果您将第一个和最后一个结绘制成:
print timeList[-1]
print knots[:4]
print knots[-4:]
你得到
>> 0.000833333333333
>> [-1.56250000e-05 -1.04166667e-05 -5.20833333e-06 0.00000000e+00]
>> [0.00083333 0.00083854 0.00084375 0.00084896]
表示实际数据范围外有六个点,前三个,后三个。这可能没问题,因为我们有 PBC,结与数据范围内的结相同,但无论如何都很奇怪。 (顺便说一句,这个例子在设置 per=0
时甚至不起作用。)此外,如果我手动设置点,设置数据范围之外的点会失败。即使使用 per=1
,如果点正好在数据范围内,它也能工作。比如:
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import splev, splrep
import scipy
x = np.linspace( 0, 1, 120 )
timeList = np.linspace( 0, 1, 15 )
signal = np.fromiter( ( np.sin( 2 * np.pi * t +.3 ) for t in timeList ), np.float )
signal[ -1 ] -= .15
myKnots=np.linspace( .05, .95, 8 )
spl = splrep( timeList, signal, t=myKnots, per=1 )
fit = splev(x,spl)
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( timeList, signal, marker='o' )
ax.plot( x, fit , 'r' )
for i in myKnots:
ax.axvline( i )
plt.show()
提供:
我们还可以看到最后一点在 PBC 中被忽略了。
那么 scipy.interpolate.splprep
在这里实际做了什么,如果 per=1
为什么它不接受类似的手动设置结的结结构。这是一个错误吗?
我希望最后一个问题的答案是 'no',否则我在这里问这个问题就错了。
为了使 splrep
使用您自己的结,您应该更改:
spl = splrep(timeList, signal, t=knots, per=1)
使用内部结:
spl = splrep(timeList, signal, t=knots[4:-4], per=1)
这个界面不是很直观,但它出现在the documentation中(尽管在task
参数下而不是在t
下)(我的重点):
If task=-1 find the weighted least square spline for a given set of knots, t. These should be interior knots as knots on the ends will be added automatically.
那里也记录了额外的结。这与结数与系数数之间的一般联系是一致的:number_of_knots = number_of_coefficients + degree + 1
.
PBC 构造中的最后一个点被忽略(如图所示),因为周期定义要求周期中最后一个点的值等于第一个点。这是在 per
参数下记录的 "Values of y[m-1] and w[m-1] are not used."
(同样,我猜这不是一个非常直观的界面,但这可能是 Fortran 代码的实现方式)。