python - 噪声数据中的梯度，python

Question

我有来自宇宙射线探测器的能谱。光谱遵循指数曲线，但其中会有宽(可能非常小)的肿 block 。显然，数据包含噪声元素。

我正在尝试平滑数据，然后绘制其梯度。到目前为止，我一直在使用 scipy sline 函数对其进行平滑处理，然后使用 np.gradient()。

从图中可以看出，梯度函数的方法是找出每个点之间的差异，并没有很清楚地表现出 block 状。

我基本上需要一个平滑的梯度图。任何帮助都会很棒!

我试过两种样条方法:

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

但是，它只是再次返回相同的图形!

Answer 1

在此发布了一个有趣的方法:Numerical Differentiation of Noisy Data .它应该可以很好地解决您的问题。更多细节在另一个accompanying paper中给出。 .作者还给出Matlab code that implements it ;替代品 implementation in Python也可用。

如果你想追求样条插值方法，我建议调整scipy.interpolate.UnivariateSpline() s.

另一种解决方案是通过卷积(例如使用高斯)来平滑您的函数。

我链接到的论文声称可以防止卷积方法产生的一些伪像(样条方法可能会遇到类似的困难)。