使用scipy.interpolate的样条表示：低振幅快速振荡函数的插值效果差

Question

我需要（通过数值方法）计算一个函数的第一和第二导数。为此，我尝试使用 splrep 和 UnivariateSpline 来创建样条曲线，以便对这个函数进行插值，从而得到导数。

但是，似乎在样条曲线的表示上存在一个固有的问题，特别是对于那些幅度在10^-1或更低，并且（快速）振荡的函数。

举个例子，下面的代码是用来在区间（0,6*pi）上创建正弦函数的样条表示（所以这个函数只振荡三次）：

import scipy
from scipy import interpolate
import numpy
from numpy import linspace
import math
from math import sin

k = linspace(0, 6.*pi, num=10000) #interval (0,6*pi) in 10'000 steps
y=[]
A = 1.e0 # Amplitude of sine function

for i in range(len(k)):

  y.append(A*sin(k[i]))

tck =interpolate.UnivariateSpline(x, y, w=None, bbox=[None, None], k=5, s=2)
M=tck(k)

下面是当A = 1.e0和A = 1.e-2时M的结果

https://i.stack.imgur.com/HOmtf.png 振幅 = 1

https://i.stack.imgur.com/W8r3l.png 振幅 = 1/100

很明显，通过样条曲线创建的插值函数完全不正确！第二张图甚至没有振荡出正确的频率。

有没有人对这个问题有什么见解？或者知道在numpy/scipy中还有其他创建样条曲线的方法吗？

谢谢，
Rory

Answer 1

问题在于如何选择合适的 s 参数值。这个值取决于数据的缩放情况。

仔细阅读文档后，可以推测出这个参数应该选择在 s = len(y) * np.var(y) 附近，也就是数据点的数量乘以方差。举个例子，选择 s = 0.05 * len(y) * np.var(y) 会得到一个平滑样条曲线，这个曲线不受数据缩放或数据点数量的影响。

编辑：当然，s 的合理值还要考虑数据中的噪声水平。文档似乎建议选择 s 在范围 (m - sqrt(2*m)) * std**2 <= s <= (m + sqrt(2*m)) * std**2 之间，其中 std 是与您想要平滑的“噪声”相关的标准差。

Answer 2

我猜你的问题可能是因为别名效应。

在你的例子中，x是什么？

如果你在插值时使用的x值之间的间距比你原始点的间距要大，你就会失去一些频率信息。这和插值的方式无关，而是因为你在降低采样率。

算了，别管上面关于别名效应的事。这在这个情况下不适用（不过我还是不明白你例子中的x是什么……）

我刚意识到，当你使用非零的平滑因子（s）时，你是在原始输入点上评估你的数据。

根据定义，平滑处理不会完全符合数据。试试把s=0放进去。

举个简单的例子：

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

x = np.linspace(0, 6.*np.pi, num=100) #interval (0,6*pi) in 10'000 steps
A = 1.e-4 # Amplitude of sine function
y = A*np.sin(x)

fig, axes = plt.subplots(nrows=2)
for ax, s, title in zip(axes, [2, 0], ['With', 'Without']):
    yinterp = interpolate.UnivariateSpline(x, y, s=s)(x)
    ax.plot(x, yinterp, label='Interpolated')
    ax.plot(x, y, 'bo',label='Original')
    ax.legend()
    ax.set_title(title + ' Smoothing')

plt.show()

你之所以只在低幅度下清楚地看到平滑效果，是因为平滑因子的定义方式。想了解更多细节，可以查看scipy.interpolate.UnivariateSpline的文档。

即使幅度更高，如果使用平滑处理，插值后的数据也不会和原始数据完全一致。

例如，如果我们在上面的代码示例中把幅度（A）改成1.0，我们仍然会看到平滑的效果……