使用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)

下面是M的结果,A = 1.e0,A = 1.e-2

http://i.imgur.com/uEIxq.png幅度= 1

http://i.imgur.com/zFfK0.png振幅= 1/100

很明显,样条函数创建的插值函数完全不正确!第二个图甚至没有振荡正确的频率.

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

干杯,罗里

Answer 1

我猜你的问题是由于别名造成的.

什么是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,我们仍然会看到平滑的效果......