噪声数据中的渐变，python
我有一个宇宙射线探测器的能谱。频谱遵循指数曲线，但它将有广泛的(可能是非常轻微的)块。显然，这些数据包含了一些噪声元素。
我正在尝试平滑数据，然后绘制其梯度图。到目前为止，我一直使用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然而，它只是再次返回相同的图！
回答 3
Stack Overflow用户
发布于 2013-04-07 20:15:39
在上面发布了一个有趣的方法：Numerical Differentiation of Noisy Data。它会为你的问题提供一个很好的解决方案。更多细节在另一篇文章accompanying paper中给出。作者还给出了Matlab code that implements it；还有一个替代的implementation in Python。
如果你想用样条方法进行插值，我建议调整scipy.interpolate.UnivariateSpline()的平滑因子s。
另一种解决方案是通过高斯卷积(比如高斯卷积)平滑你的函数。
我链接到的论文声称要防止卷积方法出现的一些伪影(样条法可能会遇到类似的困难)。
Stack Overflow用户
发布于 2013-11-06 02:26:27
我不能保证这一点的数学正确性；看起来EOL引用的LANL的论文值得一看。无论如何，在使用splev时，我已经使用SciPy的splines的内置差异化得到了不错的结果。
%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()
Stack Overflow用户
发布于 2019-08-30 18:05:11
您也可以使用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()https://stackoverflow.com/questions/15862066
复制相似问题
腾讯云开发者
