Интеллектуальный метод обнаружения пиков
Я хотел бы обнаружить пики из этих данных с помощью python:
data = [1.0, 0.35671858559485703, 0.44709399319470694, 0.29438948200831194, 0.5163825635166547, 0.3036363865322419, 0.34031782308777747, 0.2869558046065574, 0.28190537831716, 0.2807516154537239, 0.34320479518313507, 0.21117275536958913, 0.30304626765388043, 0.4972542099530442, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.18200891715227194, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.28830608331168983, 0.057156776746163526, 0.043418555819326035, 0.022527521866967784, 0.035414574439784685, 0.062273775107322626, 0.04569227783752021, 0.04978915781132807, 0.0599089458581528, 0.05692515997545401, 0.05884619933405206, 0.0809943356922021, 0.07466587894671428, 0.08548458657792352, 0.049216679971411645, 0.04742180324984401, 0.05822208549398862, 0.03465282733964001, 0.014005094192867372, 0.052004161876744344, 0.061297263734617496, 0.01867087951563289, 0.01390993522118277, 0.021515814095838564, 0.025260618727204275, 0.0157022555745128, 0.041999490119172936, 0.0441231248537558, 0.03079711140612242, 0.04177946154195037, 0.047476050325192885, 0.05087930020034335, 0.03889899267688956, 0.02114033158686702, 0.026726959895528927, 0.04623461918879543, 0.05426474524591766, 0.04421866212189775, 0.041911901968304605, 0.019982199103543322, 0.026520396430805435, 0.03952286472888431, 0.03842652984978244, 0.02779682035551695, 0.02043518392128019, 0.07706934170969436]
Вы можете построить это:
import matplotlib.pyplot as plt
plt.plot(data)
Пики, которые я хотел бы автоматически обнаружить, я обвел красным.
ХАРАКТЕРИСТИКА ПИКОВ:
Меня интересует поиск пиков, после которых для некоторых точек данных (например, 3-4) сигнал становится относительно гладким. Под плавностью я подразумеваю, что изменения амплитуд сравнимы между точками данных после пика. Я предполагаю, что в более математических терминах это означает: пики, после которых для некоторых точек данных, если бы вы подогнали линейную линию, наклон был бы близок к 0.
Что я пробовал до сих пор:
Я думал, что разница между элементами (добавив 0, чтобы иметь одинаковую длину) намного лучше покажет пики:
diff_list = []
# Append 0 to have the same length as data
data_d = np.append(data,0)
for i in range(len(data)):
diff = data_d[i]-data_d[i+1]
# If difference is samller than 0, I set it to 0 -> Just interested in "falling" peaks
if diff < 0:
diff = 0
diff_list= np.append(diff_list,diff)
Когда я рисую diff_list, это выглядит уже намного лучше:
Однако простой алгоритм обнаружения пика порогового значения не работает, потому что шум в первой секции имеет ту же амплитуду, что и пик позже.
Итак, мне нужен алгоритм, который надежно найдет пики, или метод резкого уменьшения шума без сильного демпфирования пиков и, что наиболее важно, без их смещения. У кого-нибудь есть идея?
РЕДАКТИРОВАТЬ 1:
Я наткнулся на этот блог и попробовал метод это:
peaks_d = detect_peaks(diff_list, mph=None, mpd=4, threshold=0.1, edge='falling', kpsh=False, valley=False, show=False, ax=None)
plt.plot(diff_list)
plt.plot(peaks_d[:-1], diff_list[peaks_d[:-1]], "x")
plt.show()
... но я получил:
... так что на самом деле я считаю, что мне нужна дополнительная предварительная обработка.
Обновлено еще раз:
Итак, я попытался вычислить градиент:
plt.plot(np.gradient(data))
Однако градиент внутри шума сравним с одним из пиков:
Что можно использовать:
-> Шум: существует множество точек с одинаковой амплитудой, расположенных рядом друг с другом. Возможно, можно было бы обнаружить эти области и отфильтровать их (т.е. установить для них значение 0)
РЕДАКТИРОВАТЬ 3:
Я пытался следовать Этот метод:
# Data
y = diff_list.tolist()
# Settings: lag = 30, threshold = 5, influence = 0
lag = 10
threshold = 0.1
influence = 1
# Run algo with settings from above
result = thresholding_algo(y, lag=lag, threshold=threshold, influence=influence)
# Plot result
plt.plot(result["signals"])
Однако я получаю:
РЕДАКТИРОВАТЬ 4:
На основе комментария @Jussi Nurminen:
compute the absolute value of the derivative, average it for some samples after the peak and see if the resulting value is "small enough". Of course you have to detect all candidate peaks first. For that, you could use scipy.signal.argrelextrema which detects all local maxima.
import scipy.signal as sg
max_places = (np.array(sg.argrelmax(diff_list))[0]).tolist()
plt.plot(diff_list)
plt.plot(max_places, diff_list[max_places], "x")
plt.show()
peaks = []
for check in max_places:
if check+5 < len(diff_list):
gr = abs(np.average(np.gradient(diff_list[check+1: check+5])))
if gr < 0.01:
peaks.append(check)
plt.plot(diff_list)
plt.plot(peaks[:-1], diff_list[peaks[:-1]], "x")
plt.show()
РЕДАКТИРОВАТЬ 5:
Вот аналогичные данные для проверки любого алгоритма:
data2 = [1.0, 0.4996410902399043, 0.3845950995707942, 0.38333441505960125, 0.3746384799687852, 0.28956967636700215, 0.31468441185494306, 0.5109048238958792, 0.5041481423190644, 0.41629226772762024, 0.5817609846838199, 0.3072152962171569, 0.5870564826981163, 0.4233247394608264, 0.5943712016644392, 0.4946091070102793, 0.36316740988182716, 0.4387555870158762, 0.45290920032442744, 0.48445358617984213, 0.8303387875295111, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.29678306715530073, 0.10146278147135124, 0.10120143287506084, 0.10330143251114839, 0.0802259786323741, 0.06858944745608002, 0.04600545347437729, 0.014440053029463367, 0.019023393725625705, 0.045201054387436344, 0.058496635702267374, 0.05656947149500993, 0.0463696266116956, 0.04903205756575247, 0.02781307505224703, 0.044280150764466876, 0.03746976646628557, 0.021526918040025544, 0.0038244080425488013, 0.008617907527160991, 0.0112760689575489, 0.009157686770957874, 0.013043259260489413, 0.01621417695776057, 0.016502269315028423, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3210019708643843, 0.11441868790191953, 0.12862935834434436, 0.08790971283197381, 0.09127615787146504, 0.06360039847679771, 0.032247149009635476, 0.07225952295002563, 0.095632185243862, 0.09171396569135751, 0.07935726217072689, 0.08690487354356599, 0.08787369092132288, 0.04980466729311508, 0.05675819557118429, 0.06826614158574265, 0.08491084598657253, 0.07037944101030547, 0.06549710463329293, 0.06429902857281444, 0.07282805735716101, 0.0667027178198566, 0.05590329380937183, 0.05189048980041104, 0.04609913889901785, 0.01884014489167378, 0.02782496113905073, 0.03343588833365329, 0.028423168106849694, 0.028895130687196867, 0.03146961123393891, 0.02287127937400026, 0.012173655214339595, 0.013332601407407033, 0.014040309216796854, 0.003450677642354792, 0.010854992025496528, 0.011804042414950701, 0.008100266690771957, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.18547803170164875, 0.008457776819382444, 0.006607607749756658, 0.008566964920042127, 0.024793283595437438, 0.04334031667011553, 0.012330921737457376, 0.00994343436054472, 0.008003962298473758, 0.0025523166577987263, 0.0009309499302016907, 0.0027602202618852126, 0.0034442123857338675, 0.0006448449815386562, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Используя ответ @jojo и выбрав соответствующие параметры (dy_lim = 0.1 и di_lim = 10, результат близок, но были добавлены некоторые точки, которые не должны быть пиками.
РЕДАКТИРОВАТЬ 5:
Тем не менее, другой случай.
data = [1.0, 0.0, -0.0, 0.014084507042253521, 0.0, -0.0, 0.028169014084507043, 0.0, -0.0, 0.014084507042253521, 0.0, 0.0, 0.39436619718309857, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.7887323943661971, 0.11267605633802817, 0.2535211267605634, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.4084507042253521, -0.0, 0.04225352112676056, 0.014084507042253521, 0.014084507042253521, 0.0, 0.28169014084507044, 0.04225352112676056, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.5633802816901409, -0.0, -0.0, -0.0, -0.0, 0.0, 0.08450704225352113, -0.0, -0.0, -0.0, -0.0, 0.0, 0.30985915492957744, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.9295774647887324]
Здесь почти все пики детектируются правильно, кроме одного.
Ответы 2
Возможно, вы захотите попробовать scipy.signal.find_peaks, который позволяет указывать различные критерии (заметность, ширина, высота и т. д.). Однако сначала вы должны четко определить, каковы ваши критерии «пика». Недостаточно сказать, что вам нужны некоторые пики, но не некоторые другие пики - между ними должна быть какая-то разница, которую алгоритм может обнаружить.
Привет ! Хороший ответ, пожалуйста, всегда думайте о предоставлении примера примера и работе над решением с данными ОП.
Ну почему data[13] == 0.497 пик, а data[4] == 0.516 нет? Пока OP не сможет ответить на этот фундаментальный вопрос таким образом, чтобы его можно было сформулировать алгоритмически, примеры кода бесполезны (IMO).
Спасибо за ваш ответ/комментарий. Чтобы уточнить: мне интересно найти пики, после которых для некоторых точек данных (т.е. 3-4) сигнал относительно гладкий. Под плавностью я подразумеваю, что изменения амплитуд сравнимы между точками данных.
@henry Хорошо, но вы отметили data[0] как пик, и после этого сигнал не выглядит очень гладким.
Хорошая точка зрения. Я забыл упомянуть, что первая точка данных по умолчанию считается пиком.... но это просто сбивает с толку, поэтому я быстро отредактирую изображение.
@henry сигнал после последнего пика data[122] тоже не выглядит супергладким. Но если вы хотите использовать «гладкость после пика» в качестве критерия, вы можете просто вычислить абсолютное значение производной, усреднить его для некоторых выборок после пика и посмотреть, является ли полученное значение «достаточно маленьким». Конечно, вы должны сначала обнаружить все пики-кандидаты. Для этого вы можете использовать scipy.signal.argrelextrema, который обнаруживает все локальные максимумы.
@JussiNurminen Хорошая идея. Попробую это.
Это прагматичное решение, поскольку, как я это вижу (пожалуйста, поправьте меня, если я ошибаюсь), вы хотите найти каждый пик после/до «гладкого» или нулевого периода.
Вы можете сделать это, просто проверив такие периоды и сообщив об их начале и конце.
Вот очень простая реализация, позволяющая указать, что квалифицируется как smooth период (здесь я использую изменение менее 0,001 в качестве условия):
dy_lim = 0.001
targets = []
in_lock = False
i_l, d_l = 0, data[0]
for i, d in enumerate(data[1:]):
if abs(d_l - d) > dy_lim:
if in_lock:
targets.append(i_l)
targets.append(i + 1)
in_lock = False
i_l, d_l = i, d
else:
in_lock = True
А затем рисуем targets:
plt.plot(range(len(data)), data)
plt.scatter(targets, [data[t] for t in targets], c='red')
plt.show()
Ничего особенно сложного, но он находит указанный вами пик.
Увеличение значения dy_lim позволит вам найти больше пиков. Также вы можете указать минимальную длину гладкого периода, вот как это может выглядеть (опять же, просто грубая реализация):
dy_lim = 0.001
di_lim = 50
targets = []
in_lock = False
i_l, d_l = 0, data[0]
for i, d in enumerate(data[1:]):
if abs(d_l - d) > dy_lim:
if in_lock:
in_lock = False
if i - i_l > di_lim:
targets.append(i_l)
targets.append(i + 1)
i_l, d_l = i, d
else:
in_lock = True
При этом вы не получите первое очко, так как разница между первым и вторым больше, чем di_lim=50.
Обновление для 2-го набора данных:
Это становится немного сложнее, так как теперь после пика происходит постепенное снижение, что приводит к медленному агрегированию различий, достаточному для достижения dy_lim, что приводит к ложному сообщению алгоритма о новой цели. Итак, вам нужно проверить, действительно ли эта цель является пиком, и сообщить только тогда.
Вот грубая реализация того, как этого добиться:
dy_lim = 0.1
di_lim = 5
targets = []
in_lock = False
i_l, d_l = 0, data[0]
for i, d in enumerate(data[1:]):
if abs(d_l - d) > dy_lim:
if in_lock:
in_lock = False
if i - i_l > di_lim:
# here we check whether the start of the period was a peak
if abs(d_l - data[i_l]) > dy_lim:
# assure minimal distance if previous target exists
if targets:
if i_l - targets[-1] > di_lim:
targets.append(i_l)
else:
targets.append(i_l)
# and here whether the end is a peak
if abs(d - data[i]) > dy_lim:
targets.append(i + 1)
i_l, d_l = i, d
else:
in_lock = True
В итоге вы получите следующее:
Общее примечание: Здесь мы используем восходящий подход: у вас есть конкретная функция, которую вы хотите обнаружить, поэтому вы пишете для этого специальный алгоритм.
Это может быть очень эффективно для простых задач, однако уже в этом простом примере мы понимаем, что если есть новые функции, с которыми алгоритм должен справляться, нам нужно его адаптировать. Если текущая сложность — это все, что есть, то все в порядке. Но если данные представляют собой другие шаблоны, то вы снова окажетесь в ситуации, когда вам нужно добавить дополнительные условия, и алгоритм будет становиться все более и более сложным, так как ему нужно справляться с дополнительной сложностью. Если вы окажетесь в такой ситуации, вы можете подумать о переключении передач и применить более искренний подход. В этом случае есть много вариантов, один из способов — работать с разницей исходных данных с отфильтрованной версией Савицкий-Голай, но это лишь одно из многих предложений, которые можно здесь сделать.
Спасибо за Ваш ответ. Это хорошая идея, но я на самом деле искал что-то более надежное. Это работает для этих данных, но у меня также есть другие данные, где есть, например, 4-й пик после немного более грубого растяжения данных (т.е. ненулевого)
вы можете справиться с шероховатостью, указав dy_lim, я добавлю пример.
Хорошо, отлично! Может быть, вы можете попробовать свой алгоритм с похожим, но другим набором данных, который я только что добавил. (см. обновленный вопрос)
Можешь попробовать еще раз сейчас? Теперь все должно быть в порядке. Спасибо большое.
Пожалуйста, взгляните на мой обновленный ответ. Может быть, если добавить минимальную высоту пика, то таких проблем (как видно из моего вопроса) не возникнет.
вы можете просто удалить оператор targets.append(i_l). Вы пропустите первую точку, но если присмотреться, первая точка также отличается от других, так как раньше у нее не было гладкой области.
Хм, у вас есть идеи, как это можно исправить?
Отличный ответ. +1 У вас похожая проблема, и я пытался понять ваш код, но немного запутался с переменной in_lock. Не могли бы вы добавить еще несколько комментариев, чтобы объяснить ваш код? Извините за беспокойство, если у вас нет времени, ничего страшного.
Большое спасибо за ваше исправление!!
@james Вот краткое объяснение идеи, лежащей в основе on_lock: поскольку алгоритм работает с набором данных, и мы хотим иметь возможность определить, находится ли он в настоящее время в гладкой области или нет, мы устанавливаем флаг in_lock, который вы сохраняете. на True до тех пор, пока разница между текущей точкой данных и точкой данных после последнего «большого изменения» меньше некоторого порогового значения, dy_lim. Если вы столкнулись с еще одним «большим изменением» и у вас установлен флаг True, это означает, что вы были в «гладком» периоде, который вы сейчас покидаете, поэтому вы хотите сообщить об этом. надеюсь, что помог, дайте мне знать, если нет!
@jojo Привет! Извините, что беспокою. Я столкнулся с набором данных, который, как мне кажется, будет легко обрабатывать вашим кодом. К сожалению, он обнаруживает небольшой пик там, где его быть не должно. Вы видите какой-нибудь способ, как это исправить?
@henry хм, все будет зависеть от того, как именно вы бы описали фактические пики и, как следствие, как бы вы описали пик, который не должен быть найден. Существует ли минимальный размер пика? Вы всегда можете настроить минимальный размер пика, настроив dy_lim. Просто взглянув на ваш сюжет, может быть, dy_lim=0.15 подойдет?
Спасибо за ваш быстрый ответ. Я думаю, что в этом случае можно было бы ввести минимальное расстояние между пиками. Например, я совершенно уверен, что в моем наборе данных пики разнесены более чем на 5 (значение x).
@ Генри, действительно, это должно сработать. Я добавлю небольшую модификацию в последнюю версию, чтобы установить ограничение на значение x (здесь это индекс, поэтому я называю его di_lim). Обратите внимание, что использование dy_lim=0.15, кажется, работает!
Большое спасибо ! С нетерпением жду. :) Да, dy_lim=0.15 работает, но у меня есть другие датасеты, где пики выражены гораздо меньше по высоте. Большое спасибо за Вашу помощь.
@henry хорошо, добавлено. теперь вы можете поэкспериментировать с di_lim, чтобы установить условие для x_values, и с dy_lim для условий по размеру пика. :)
Очень хорошо !! Большое спасибо !
Есть хороший фильтр, который использует Z-показатели для обнаружения пиков или, возможно, выбросов. Он сигнализирует, когда следующая точка превышает статистический порог, и хорошо объясняется многочисленными программными реализациями: stackoverflow.com/a/22640362/7621907