我有一个宇宙射线探测器的能谱。光谱遵循指数曲线,但会有很宽(可能很小)的肿块。显然,这些数据包含了一个噪声元素。
我试着把数据平滑,然后画出它的梯度。 到目前为止,我一直在使用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
但是,它只是返回相同的图形!
我不保证这篇文章的数学有效性;看起来EOL引用的LANL的论文值得一看。不管怎样,在使用
splev
时,我使用SciPy样条的内置微分得到了不错的结果。在这个:Numerical Differentiation of Noisy Data上发布了一个有趣的方法。它应该能给你一个解决问题的好办法。更多细节在另一个accompanying paper中给出。作者还提供了Matlab code that implements it;另一个implementation in Python也可用。
如果您想使用样条函数进行插值,我建议您调整
scipy.interpolate.UnivariateSpline()
的平滑因子s
。另一个解决方案是通过卷积来平滑函数(比如用高斯函数)。
我在论文中联系了一些声明,以防止出现卷积方法(样条方法可能会遇到类似的困难)的一些伪影。
您也可以使用^{} 。
结果
示例
相关问题 更多 >
编程相关推荐