嘈雜的數據中的漸變，python

Question

我有一個宇宙射線探測器的能譜。光譜遵循指數曲線，但它會有廣泛的（也許非常輕微）腫塊。顯然，這些數據包含了一些噪音。

我想平滑數據，然後繪製其漸變。 到目前爲止，我一直在使用scipy sline函數來平滑它，然後使用np.gradient（）。

正如您從圖片中看到的，梯度函數的方法是找出每個點之間的差異，並且它不會非常清楚地顯示塊。

我基本上需要一個平滑的梯度圖。任何幫助將是驚人的！

我已經試過2樣條方法：

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。

另一種解決方案是通過卷積（比如用高斯表示）來平滑你的函數。

我鏈接到的論文聲稱可以防止某些與卷積方法（樣條方法可能會遇到類似困難）相關的工件。

Answer 2

我不會擔保此數學有效性;它看起來像是LANL的論文，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()