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

Question

我需要（以数字方式）计算一个函数的一阶和二阶导数，我尝试使用 splrep 和 UnivariateSpline 来创建样条线以进行插值函数取导数。

然而，对于幅度为 10^-1 或更低的函数（快速）振荡，样条表示本身似乎存在固有问题。

作为示例，请考虑使用以下代码在区间 (0,6*pi) 上创建正弦函数的样条表示（因此该函数仅振荡 3 次）：

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 时 M 的结果A = 1.e-2

https://i.sstatic.net/HOmtf.png 幅度= 1

https://i.sstatic.net/W8r3l.png 振幅 = 1/100

显然由样条线创建的插值函数完全不正确！第二张图甚至没有振荡正确的频率。

有人对这个问题有任何见解吗？或者知道在 numpy/scipy 中创建样条线的另一种方法吗？

干杯， 罗里

Answer 1

<罢工>
我猜你的问题是由于别名造成的。

您的示例中的 x 是什么？

如果您插值的 x 值的间距小于原始点的间距，您本质上就会丢失频率信息。这完全独立于任何类型的插值。这是下采样所固有的。

别介意上面关于别名的事情。它不适用于这种情况（尽管我仍然不知道您的示例中的 x 是什么...

我刚刚意识到您正在评估您的观点at当您使用非零平滑因子 (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，我们仍然会看到平滑的效果...

Answer 2

问题在于为 s 参数选择合适的值。其值取决于数据的缩放比例。

仔细阅读文档，可以推断出应该围绕 s = len(y) * np.var(y) 来选择参数，即数据点数量 * 方差。以 s = 0.05 * len(y) * np.var(y) 为例，给出一个不依赖于数据缩放或数据点数量的平滑样条曲线。

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