Как найти оптимизированные или правильные пики

#python #scipy #signal-processing

Вопрос:

У меня есть следующий график

введите описание изображения здесь

Я использую python scipy.signal.find_peaks, чтобы найти пики. Но я не уверен, как мне это сделать. Я сделал следующее :

 per = np.percentile(x,[70])
peaks_control = findPeaks(x, per[0])
 

где x — сигнал

массив([1.07541259 е 09, 1.13851049 е 09, 1.19241492 е 09, 1.23706527 е 09, 1.27240093 е 09, 1.29836131 е 09, 1.31217483 е 09, 1.32037296 е 09, 1.31908858 е 09, 1.30896503 е 09, 1.29216550 е 09, 1.26958042 е 09, 1.24561632 е 09, 1.21202121 е 09, 1.16869371 е 09, 1.11054499 е 09, 1.04006154 е 09, 9.65663403 е 08, 8.87706760 е 08, 8.09340093 е 08, 7.37568765 е 08, 6.79736364 е 08, 6.38576457 е 08, 6.06062937 е 08, 5.80650350 е 08, 5.55089744 е 08, 5.36334499 е 08, 5.20236597 е 08, 5.06529837 е 08, 4.91825175 е 08, 4.77937063 е 08, 4.65475058 е 08, 4.56520513 е 08, 4.48393240 е 08, 4.41944988 е 08, 4.34822844 е 08, 4.33688578e 08, 4.33451049 е 08, 4.36256177 е 08, 4.33553613 е 08, 4.29191142 е 08, 4.28492541 е 08, 4.24465967 е 08, 4.20074825 е 08, 4.19935897 е 08, 4.16652681 е 08, 4.12419580 е 08, 4.11747552 е 08, 4.08801166 е 08, 4.02351981 е 08, 3.99620513 е 08, 3.98716550 е 08, 3.46023077 е 08, 3.53969464 е 08, 4.17131235 е 08, 5.19363869 е 08, 6.50956410 е 08, 8.01530303 е 08, 9.50162937 е 08, 1.08249790 е 09, 1.18242378 е 09, 1.22732168 е 09, 1.20123077 е 09, 1.21067599 е 09, 1.21556410 е 09, 1.21272261 е 09, 1.20310023 е 09, 1.18692774 е 09, 1.16694033 е 09, 1.14330117 е 09, 1.11635338 е 09, 1.07947529 е 09, 1.03222145e 09, 9.73427972 е 08, 9.08558974 е 08, 8.39966200 е 08, 7.70457343 е 08, 7.04976224 е 08, 6.49436131 е 08, 6.02085548 е 08, 5.68915385 е 08, 5.41638928 е 08, 5.18758741 е 08, 5.01973660 е 08, 4.88766667 е 08, 4.77643823 е 08, 4.65681818 е 08, 4.56193240 е 08, 4.46851515 е 08, 4.36135198 е 08, 4.32282984 е 08, 4.27913520 е 08, 4.23408625 е 08, 4.24119580 е 08, 4.22399068 е 08, 4.22415385 е 08, 4.20193939 е 08, 4.17638462 е 08, 4.14822378 е 08, 4.10636364 е 08, 4.08388345 е 08, 4.04844522 е 08, 4.00571562 е 08, 4.00841026 е 08, 4.00764802 е 08, 4.00432867 е 08, 4.00336364 е 08, 4.00724709 е 08, 4.03048019e 08, 3.57437995 е 08, 3.62371096 е 08, 4.16658741 е 08, 5.10148019 е 08, 6.31750117 е 08, 7.65175991 е 08, 8.96832168 е 08, 1.01666597 е 09, 1.10373263 е 09, 1.14380816 е 09, 1.11629790 е 09, 1.12228904 е 09, 1.12378788 е 09, 1.11974825 е 09, 1.10812774 е 09, 1.09125035 е 09, 1.07033566 е 09, 1.04667389 е 09, 1.02016830 е 09, 9.86036830 е 08, 9.42176457 е 08, 8.88900233 е 08, 8.27962005 е 08, 7.64362238 е 08, 7.00755245 е 08, 6.42390909 е 08, 5.92395338 е 08, 5.52426107 е 08, 5.26319114 е 08, 5.03317249 е 08, 4.85524942 е 08, 4.70421911 е 08, 4.59389510 е 08, 4.51644988 е 08, 4.46288578 е 08, 4.41076923e 08, 4.37533566 е 08, 4.31993007 е 08, 4.28625641 е 08, 4.25406294 е 08, 4.21161538 е 08, 4.19049650 е 08, 4.16719347 е 08, 4.13124242 е 08, 4.08404429 е 08, 4.06154545 е 08, 4.03386014 е 08, 4.00980420 е 08, 3.99442657 е 08, 3.97792075 е 08, 3.95606527 е 08, 3.97922378 е 08, 3.98345221 е 08, 3.96253613 е 08, 3.95703030 е 08, 3.96108392 е 08, 3.67136830 е 08, 3.58382051 е 08, 3.95844289 е 08, 4.70853846 е 08, 5.76629837 е 08, 6.97682284 е 08, 8.21169930 е 08, 9.32588112 е 08, 1.01885804 е 09, 1.06315152 е 09, 1.05128159 е 09, 1.03944545 е 09, 1.03769580 е 09, 1.03132145 е 09, 1.02008601 е 09, 1.00327389e 09, 9.85387646 e 08, 9.66403030 e 08, 9.44620746 e 08, 9.18596737 e 08, 8.82269697 e 08, 8.37750816 e 08, 7.84877156 e 08, 7.27590443 e 08, 6.70183217 e 08, 6.14567832 e 08, 5.67404895 e 08, 5.30862471 e 08, 5.03108625 e 08, 4.84348718 e 08, 4.68116550 e 08, 4.55809907 e 08, 4.46616783 e 08, 4.39725175 e 08, 4.34323077 e 08])

Вершины, которые я получаю, примыкают друг к другу, так как я вижу, что на втором, третьем и четвертом участках пиков есть небольшие выпуклости. введите описание изображения здесь

Как мне это рассчитать и игнорировать такие соседние? Чтобы рассчитать ширину, проминекне и т. Д., Мне Нужно рассчитать пики. Если бы я уже знал это, я мог бы установить какой-то порог.

Комментарии:

1. Если вы используете find_peaks scipy (см. Документацию , которую вы можете установить distance . Если вы не знаете расстояния, вы можете выполнить периодограмму и установить период 1-й гармоники в качестве расстояния

2. Я новичок в этом деле. не могли бы вы привести мне, пожалуйста, какой-нибудь пример

Ответ №1:

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

Итак, давайте создадим некоторые зашумленные данные

 import numpy as np
from scipy.signal import find_peaks, periodogram
import matplotlib.pyplot as plt

size = 100
a = np.linspace(1, .5, size)
x = np.linspace(0, 50, size)
np.random.seed(0)
y = a * np.sin(x)   np.random.normal(0, .1, size)   5
 

введите описание изображения здесь

теперь давайте попробуем найти вершины с find_peaks помощью scipy.signal

 peaks = find_peaks(y)[0]
plt.plot(x, y)
plt.plot(x[peaks], y[peaks], marker='o', ls='none')
plt.show()
 

введите описание изображения здесь

как вы можете видеть, есть некоторые «неправильные» пики. Нам нужно задать distance аргумент в find_peaks (см. Документацию).

Давайте предположим, что мы не знаем расстояния между вершинами. В этом случае мы видим, что данные являются периодическими. Таким образом, мы можем найти период с помощью периодограммы и использовать период в качестве расстояния в find_peaks

 _f, _p = periodogram(y, nfft=2**6)
# calculate the sample rate of x
sample_rate = 1 / np.median(np.diff(x))
periods = 1 / _f[1:] / sample_rate
density = _p[1:] / _p[1:].max()
max_density_idx = density.argmax()

period = periods[max_density_idx]

plt.semilogx(periods, density)
plt.scatter(period, density[max_density_idx], color='r')
plt.title(f"period {period:.2f}")
plt.show()
 

введите описание изображения здесь

Теперь мы можем использовать точку в качестве distance аргумента в find_peaks

 peaks = find_peaks(y, distance=period)[0]
plt.plot(x, y)
plt.plot(x[peaks], y[peaks], marker='o', ls='none')
plt.show()
 

введите описание изображения здесь


Обновить

В вашем конкретном случае все немного по-другому.

Определите сигнал (я буду называть переменные X и Y )

 Y = np.array([1.07541259e 09, 1.13851049e 09, 1.19241492e 09, 1.23706527e 09, 1.27240093e 09, 1.29836131e 09, 1.31217483e 09, 1.32037296e 09, 1.31908858e 09, 1.30896503e 09, 1.29216550e 09, 1.26958042e 09, 1.24561632e 09, 1.21202121e 09, 1.16869371e 09, 1.11054499e 09, 1.04006154e 09, 9.65663403e 08, 8.87706760e 08, 8.09340093e 08, 7.37568765e 08, 6.79736364e 08, 6.38576457e 08, 6.06062937e 08, 5.80650350e 08, 5.55089744e 08, 5.36334499e 08, 5.20236597e 08, 5.06529837e 08, 4.91825175e 08, 4.77937063e 08, 4.65475058e 08, 4.56520513e 08, 4.48393240e 08, 4.41944988e 08, 4.34822844e 08, 4.33688578e 08, 4.33451049e 08, 4.36256177e 08, 4.33553613e 08, 4.29191142e 08, 4.28492541e 08, 4.24465967e 08, 4.20074825e 08, 4.19935897e 08, 4.16652681e 08, 4.12419580e 08, 4.11747552e 08, 4.08801166e 08, 4.02351981e 08, 3.99620513e 08, 3.98716550e 08, 3.46023077e 08, 3.53969464e 08, 4.17131235e 08, 5.19363869e 08, 6.50956410e 08, 8.01530303e 08, 9.50162937e 08, 1.08249790e 09, 1.18242378e 09, 1.22732168e 09, 1.20123077e 09, 1.21067599e 09, 1.21556410e 09, 1.21272261e 09, 1.20310023e 09, 1.18692774e 09, 1.16694033e 09, 1.14330117e 09, 1.11635338e 09, 1.07947529e 09, 1.03222145e 09, 9.73427972e 08, 9.08558974e 08, 8.39966200e 08, 7.70457343e 08, 7.04976224e 08, 6.49436131e 08, 6.02085548e 08, 5.68915385e 08, 5.41638928e 08, 5.18758741e 08, 5.01973660e 08, 4.88766667e 08, 4.77643823e 08, 4.65681818e 08, 4.56193240e 08, 4.46851515e 08, 4.36135198e 08, 4.32282984e 08, 4.27913520e 08, 4.23408625e 08, 4.24119580e 08, 4.22399068e 08, 4.22415385e 08, 4.20193939e 08, 4.17638462e 08, 4.14822378e 08, 4.10636364e 08, 4.08388345e 08, 4.04844522e 08, 4.00571562e 08, 4.00841026e 08, 4.00764802e 08, 4.00432867e 08, 4.00336364e 08, 4.00724709e 08, 4.03048019e 08, 3.57437995e 08, 3.62371096e 08, 4.16658741e 08, 5.10148019e 08, 6.31750117e 08, 7.65175991e 08, 8.96832168e 08, 1.01666597e 09, 1.10373263e 09, 1.14380816e 09, 1.11629790e 09, 1.12228904e 09, 1.12378788e 09, 1.11974825e 09, 1.10812774e 09, 1.09125035e 09, 1.07033566e 09, 1.04667389e 09, 1.02016830e 09, 9.86036830e 08, 9.42176457e 08, 8.88900233e 08, 8.27962005e 08, 7.64362238e 08, 7.00755245e 08, 6.42390909e 08, 5.92395338e 08, 5.52426107e 08, 5.26319114e 08, 5.03317249e 08, 4.85524942e 08, 4.70421911e 08, 4.59389510e 08, 4.51644988e 08, 4.46288578e 08, 4.41076923e 08, 4.37533566e 08, 4.31993007e 08, 4.28625641e 08, 4.25406294e 08, 4.21161538e 08, 4.19049650e 08, 4.16719347e 08, 4.13124242e 08, 4.08404429e 08, 4.06154545e 08, 4.03386014e 08, 4.00980420e 08, 3.99442657e 08, 3.97792075e 08, 3.95606527e 08, 3.97922378e 08, 3.98345221e 08, 3.96253613e 08, 3.95703030e 08, 3.96108392e 08, 3.67136830e 08, 3.58382051e 08, 3.95844289e 08, 4.70853846e 08, 5.76629837e 08, 6.97682284e 08, 8.21169930e 08, 9.32588112e 08, 1.01885804e 09, 1.06315152e 09, 1.05128159e 09, 1.03944545e 09, 1.03769580e 09, 1.03132145e 09, 1.02008601e 09, 1.00327389e 09, 9.85387646e 08, 9.66403030e 08, 9.44620746e 08, 9.18596737e 08, 8.82269697e 08, 8.37750816e 08, 7.84877156e 08, 7.27590443e 08, 6.70183217e 08, 6.14567832e 08, 5.67404895e 08, 5.30862471e 08, 5.03108625e 08, 4.84348718e 08, 4.68116550e 08, 4.55809907e 08, 4.46616783e 08, 4.39725175e 08, 4.34323077e 08])

X = np.arange(Y.size)
 

Поскольку Y.size это 200, а на вашем графике 200 секунд, я предполагаю, что частота дискретизации составляет 1 сек.

Если мы будем искать пики с расстоянием по умолчанию, мы найдем много нежелательных пиков

 peaks = find_peaks(Y)[0]
plt.plot(X, Y)
plt.plot(X[peaks], Y[peaks], marker='o', ls='none')
plt.show()
 

введите описание изображения здесь

Давайте сделаем периодограмму

 _f, _p = periodogram(Y, nfft=2**12)
# the sample rate of your signal
sample_rate = 1 
periods = 1 / _f[1:] / sample_rate
density = _p[1:] / _p[1:].max()
max_density_idx = density.argmax()

period = periods[max_density_idx]

p_peaks_idx = find_peaks(density)[0]

plt.semilogx(periods, density)
plt.scatter(period, density[max_density_idx], color='r')
period_peaks = []
for p_peak in p_peaks_idx:
    if density[p_peak] < .1:
        continue
    period_peaks.append(periods[p_peak])
    plt.scatter(periods[p_peak], density[p_peak])
    plt.text(periods[p_peak], density[p_peak], f"{periods[p_peak]:.1f}  ", ha='right', va='center')
plt.title('periodogram')
plt.show()
 

введите описание изображения здесь

Мы обнаружили два основных периода

 period_peaks
[56.888888888888886, 28.444444444444443]
 

Если мы используем период с более высокой плотностью (56,9, фундаментальная или 1-я гармоника), мы пропустим пик

 peaks = find_peaks(Y, distance=period_peaks[0])[0]
plt.plot(X, Y)
plt.plot(X[peaks], Y[peaks], marker='o', ls='none')
plt.show()
 

введите описание изображения здесь

Это может быть потому, что

  • у вас слишком мало наблюдений
  • периодичность не является постоянной

Если мы эмпирически вычтем величину (скажем, 10) из периода, мы найдем все пики

 peaks = find_peaks(Y, distance=period_peaks[0] - 10)[0]
plt.plot(X, Y)
plt.plot(X[peaks], Y[peaks], marker='o', ls='none')
plt.show()
 

введите описание изображения здесь

Итак, у нас есть пики на

 X[peaks]
array([  7,  61, 118, 174])
 

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

 np.diff(X[peaks])
array([54, 57, 56])