Градиент в зашумленных данных, python

У меня есть энергетический спектр с детектора космических лучей. Спектр следует экспоненциальной кривой, но в нем будут большие (и, возможно, очень небольшие) глыбы. Данные, очевидно, содержат элемент шума.

Я пытаюсь сгладить данные, а затем построить их градиент. До сих пор я использовал функцию 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

Однако он просто снова возвращает тот же график!

Данные, сглаженные данные и градиенты


person Lucidnonsense    schedule 07.04.2013    source источник


Ответы (3)


По этому поводу опубликован интересный метод: Числовое дифференцирование зашумленных данных . Это должно дать вам хорошее решение вашей проблемы. Дополнительные сведения приведены в другом, сопроводительном документе. Автор также предоставляет код Matlab, реализующий это; также доступна альтернативная реализация на Python.

Если вы хотите использовать метод интерполяции со сплайнами, я бы посоветовал настроить коэффициент сглаживания s для scipy.interpolate.UnivariateSpline().

Другим решением было бы сгладить вашу функцию с помощью свертки (скажем, с помощью гауссова).

Документ, на который я ссылался, утверждает, что предотвращает некоторые артефакты, возникающие при использовании подхода свертки (подход сплайна может страдать от аналогичных трудностей).

person Eric O Lebigot    schedule 07.04.2013
comment
Я попробовал метод численного дифференцирования: см. новое вложение - person Lucidnonsense; 07.04.2013
comment
Строка part2 = simps(u) неверна: part2 вместо этого должен быть массивом, содержащим интеграл от u от 0 до каждой абсциссы. Затем вам следует попробовать экспоненциально изменяющийся коэффициент сглаживания, чтобы найти тот, который лучше всего соответствует вашим потребностям. Если вы действительно вычисляете производную, принимая во внимание размер шага по x, то я ожидаю, что хороший коэффициент сглаживания будет около 1e6 для производной, которая находится в основном между -1 и +1 - хотя я могу ошибаться, вы можете захотеть в любом случае попробовать это значение, на всякий случай, если мой предварительный расчет верен. - person Eric O Lebigot; 08.04.2013
comment
PS: также нет необходимости предварительно сглаживать данные с помощью метода численного дифференцирования зашумленных данных. Вроде интересный метод. Я с нетерпением жду ваших результатов! Удачи… - person Eric O Lebigot; 08.04.2013
comment
PPS: Хорошим исходным предположением для решения u может быть производная (зашумленная) (а не сама функция). Фактически процедура, описанная в статье, находит производную без сглаживания функции. - person Eric O Lebigot; 09.04.2013
comment
Я не уверен, что понимаю вашу версию части 2. Я бы реализовал это как цикл for, вычисляющий каждый интеграл. Но тогда как этот новый массив будет работать с остальной частью процедуры подгонки? Кроме того, я сгладил данные, чтобы дать начальное предположение, будет ли это в корне неверным? Я не понимаю, как производная с шумом может дать хорошее предположение - person Lucidnonsense; 09.04.2013
comment
Вы действительно можете использовать циклы для интегралов part2 (или использовать правило трапеций или прямоугольников и вычислять быстрее с cumsum()). Согласно уравнению (2) в документе, на который я ссылался, вы просто используете его, как и раньше, интегрируя abs(part2-f) (разницу двух массивов одинакового размера). Что касается начального предположения, алгоритм вычисляет производную u ваших данных, а не сглаженную версию ваших данных. - person Eric O Lebigot; 10.04.2013
comment
PS: я добавил ссылку на более подробную версию алгоритма и на код Matlab, который его реализует. Научное сообщество Python, скорее всего, выиграет от реализации Python, доступной через индекс пакетов Python! - person Eric O Lebigot; 10.04.2013
comment
Спасибо за вашу помощь, но я сейчас очень занят и вряд ли буду реализовывать алгоритм. Однако, когда мой текущий проект будет завершен, я попытаюсь закончить его и опубликовать результаты. Спасибо еще раз! - person Lucidnonsense; 13.04.2013
comment
Может ли кто-нибудь посоветовать мне, как справиться с аналогичной проблемой, но с неравномерно распределенными данными? У меня есть наборы измерений для X и Y. Подойдет ли описанный в статье алгоритм? Я пытаюсь пройти через реализацию Python, но пока не нашел способа ввести в игру измерения X. - person Spu; 20.10.2015
comment
@Spu, это работает только для равномерно распределенных данных. Функция, содержащаяся в реализации Python, опубликованной здесь @EOL, принимает в качестве аргумента интервал сетки dx, который является скаляром - person David; 08.08.2017
comment
@Spu Обратите внимание, что сам алгоритм может учитывать неравномерно расположенные значения X. Ограничена только реализация. Обобщить его не составит труда. - person Eric O Lebigot; 09.08.2017

Я не буду ручаться за математическую достоверность этого; похоже, что документ из LANL, на который ссылается EOL, заслуживает внимания. Во всяком случае, я получил приличные результаты, используя встроенную дифференциацию сплайнов SciPy при использовании 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()

вывод matplotlib

person billyjmc    schedule 05.11.2013
comment
Можно ли использовать этот метод для неравномерно распределенных данных? Могу ли я добавить свои наборы измерений как для X, так и для Y? - person Spu; 20.10.2015
comment
В документации к используемой здесь функции (scipy.interpolate.splrep()) не упоминается никаких ограничений для неравномерно распределенных данных. Помимо простого просмотра документации, вы также можете попробовать сами, изменив значение x в коде. В более общем плане приветствуется, что вы прикладываете видимые усилия для ответа на свой вопрос в Stack Overflow, чтобы сэкономить время другим (и повысить вероятность того, что они потратят время на ответ на ваш вопрос). - person Eric O Lebigot; 20.10.2015
comment
@Спу, да! Я использовал splrep всего два дня назад, чтобы выполнить интерполяцию кубических b-сплайнов выборочных данных, полученных с неравномерными интервалами, чтобы я мог выполнить БПФ. - person billyjmc; 21.10.2015

Вы также можете использовать 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()
person Roald    schedule 30.08.2019