嘈杂数据中的渐变，Python

Question

嘈杂数据中的渐变，Python

pythongradientsmoothing

17

我有一个宇宙射线探测器的能谱。这个谱遵循指数曲线，但其中会有宽广（可能非常轻微）的波峰。显然，数据包含一定程度的噪音。

我试图平滑数据并绘制其梯度。目前我一直在使用scipy sline函数来平滑数据，然后使用np.gradient()函数来计算梯度。

正如您从图片中所看到的那样，梯度函数的方法是找出每个点之间的差异，并且它没有清晰地显示出波峰。

基本上，我需要一个平滑的梯度图。任何帮助都将非常惊人！

我已经尝试了两种样条方法：

def smooth_data(y,x,factor):
    print "smoothing data by interpolation..."
    xnew=np.linspace(min(x),max(x),factor*len(x))
    smoothy=spline(x,y,xnew)
    return smoothy,xnew

def smooth2_data(y,x,factor):
    xnew=np.linspace(min(x),max(x),factor*len(x))
    f=interpolate.UnivariateSpline(x,y)
    g=interpolate.interp1d(x,y)
    return g(xnew),xnew

编辑：尝试了数值微分：

def smooth_data(y,x,factor):
    print "smoothing data by interpolation..."
    xnew=np.linspace(min(x),max(x),factor*len(x))
    smoothy=spline(x,y,xnew)
    return smoothy,xnew

def minim(u,f,k):
    """"functional to be minimised to find optimum u. f is original, u is approx"""
    integral1=abs(np.gradient(u))
    part1=simps(integral1)
    part2=simps(u)
    integral2=abs(part2-f)**2.
    part3=simps(integral2)
    F=k*part1+part3
    return F


def fit(data_x,data_y,denoising,smooth_fac):
    smy,xnew=smooth_data(data_y,data_x,smooth_fac)
    y0,xnnew=smooth_data(smy,xnew,1./smooth_fac)
    y0=list(y0)
    data_y=list(data_y)
    data_fit=fmin(minim, y0, args=(data_y,denoising), maxiter=1000, maxfun=1000)
    return data_fit

然而，它只是再次返回相同的图形！

数据、平滑数据和梯度

- Lucidnonsense

3个回答

10

我不能保证这个方法在数学上的有效性；看起来EOL引用的洛斯阿拉莫斯国家实验室的论文可能值得研究。无论如何，当使用splev时，使用SciPy的样条内置微分功能，我已经得到了不错的结果。

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.interpolate import splrep, splev

x = np.arange(0,2,0.008)
data = np.polynomial.polynomial.polyval(x,[0,2,1,-2,-3,2.6,-0.4])
noise = np.random.normal(0,0.1,250)
noisy_data = data + noise

f = splrep(x,noisy_data,k=5,s=3)
#plt.plot(x, data, label="raw data")
#plt.plot(x, noise, label="noise")
plt.plot(x, noisy_data, label="noisy data")
plt.plot(x, splev(x,f), label="fitted")
plt.plot(x, splev(x,f,der=1)/10, label="1st derivative")
#plt.plot(x, splev(x,f,der=2)/100, label="2nd derivative")
plt.hlines(0,0,2)
plt.legend(loc=0)
plt.show()

matplotlib output

- billyjmc

这种方法适用于非均匀分布的数据吗？我可以添加我的X和Y测量集吗？ - Spu

这里使用的函数（scipy.interpolate.splrep()）的文档没有提到非均匀分布数据的任何限制。除了查看文档外，您还可以通过更改代码中的 x 值来尝试。更一般地说，在 Stack Overflow 上，我们欣赏您做出一些可见的努力来回答自己的问题，以节省他人的时间（并使他们更有可能花时间回答您的问题）。 - Eric O. Lebigot

1

@Spu，是的！我就在两天前使用了splrep来对非均匀间隔采集的样本数据进行三次B样条插值，以便我可以执行FFT。 - billyjmc

5

你也可以使用scipy.signal.savgol_filter。

结果

示例

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

# generate data
x = np.array(range(100))/10
y = np.sin(x) + np.array([random()*0.25 for _ in x])
dydx = scipy.signal.savgol_filter(y, window_length=11, polyorder=2, deriv=1)

# Plot result
plt.plot(x, y, label='Original signal')
plt.plot(x, dydx*10, label='1st Derivative')
plt.plot(x, np.cos(x), label='Expected 1st Derivative')
plt.legend()
plt.show()

- Roald

更新的文档链接：https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html - zabop

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Eric O Lebigot · Accepted Answer

这里有一种有趣的方法，对于嘈杂数据的数值微分。它将为您的问题提供良好的解决方案。更多细节在另一篇相关论文中给出。作者还提供了Matlab代码来实现它; 还有一个Python的可选实现。

如果您想采用样条插值的方法，我建议调整scipy.interpolate.UnivariateSpline()的平滑因子s 。

另一种解决方法是通过卷积（例如高斯卷积）平滑您的函数。

我链接的论文声称可以避免一些出现在卷积方法中的伪影效应（样条方法可能也会遇到类似的困难）。