在周期性边界条件的情况下,如何在 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 代码的实现方式)。