我有宇宙射线探测器的能谱.光谱遵循指数曲线,但它会有宽的(也可能是非常轻微的)肿块.显然,数据包含噪声元素.
我正在尝试平滑数据,然后绘制其渐变.
到目前为止,我一直在使用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,factor): xnew=np.linspace(min(x),factor*len(x)) f=interpolate.UnivariateSpline(x,y) g=interpolate.interp1d(x,y) return g(xnew),xnew
编辑:尝试数值区分:
def smooth_data(y,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
但是,它只是再次返回相同的图形!
解决方法
有一个有趣的方法发表于此:
Numerical Differentiation of Noisy Data.它应该为您提供一个很好的解决方案来解决您的问题.更多细节见另一个,accompanying paper.作者还给出了
Matlab code that implements it;替代
implementation in Python也可用.
如果你想用样条方法进行插值,我建议调整scipy.interpolate.UnivariateSpline()的平滑因子s.