Gradient trong dữ liệu nhiễu, python

Question

Tôi có phổ năng lượng từ thiết bị dò tia vũ trụ. Quang phổ theo một đường cong theo cấp số nhân nhưng nó sẽ có các khối u rộng (và có thể rất nhỏ) trong đó. Các dữ liệu, rõ ràng, có chứa một yếu tố của tiếng ồn.

Tôi đang cố gắng làm mịn dữ liệu và sau đó vẽ đường viền của nó. Cho đến nay tôi đã sử dụng hàm sline scipy để làm mịn nó và sau đó là np.gradient().

Như bạn có thể thấy từ hình, phương pháp của hàm gradient là tìm sự khác biệt giữa mỗi điểm và không hiển thị các khối rất rõ ràng.

Tôi về cơ bản cần có biểu đồ độ dốc mượt mà. Bất kỳ trợ giúp sẽ là tuyệt vời!

Tôi đã thử phương pháp 2 spline:

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 

chỉnh sửa: Cố gắng số khác biệt:

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 

Tuy nhiên, nó chỉ trả về cùng một đồ thị một lần nữa!

Answer 1

Có một phương pháp thú vị công bố về vấn đề này: Numerical Differentiation of Noisy Data. Nó sẽ cung cấp cho bạn một giải pháp tốt cho vấn đề của bạn. Chi tiết khác được cung cấp trong một số khác, accompanying paper. Tác giả cũng cung cấp cho Matlab code that implements it; một thay thế implementation in Python cũng có sẵn.

Nếu bạn muốn theo dõi nội suy với phương pháp spline, tôi khuyên bạn nên điều chỉnh hệ số làm mịn s của scipy.interpolate.UnivariateSpline().

Một giải pháp khác là làm mịn chức năng của bạn thông qua convolution (nói với Gaussian).

Bài báo tôi liên kết với các khiếu nại để ngăn chặn một số hiện vật đi kèm với phương pháp tiếp cận chập chững (cách tiếp cận spline có thể gặp khó khăn tương tự).

Answer 2

Tôi sẽ không xác minh tính hợp lệ của thuật toán này; có vẻ như bài báo từ LANL mà EOL trích dẫn sẽ đáng xem xét. Dù sao, tôi đã nhận được kết quả tốt bằng cách sử dụng tính năng phân biệt được tích hợp sẵn của SpiPy khi sử dụng splev.

%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()