Algorytm znajdowania pików dla Python / SciPy


136

Mogę coś napisać samodzielnie, znajdując przecięcia zerowe pierwszej pochodnej lub czegoś w tym rodzaju, ale wydaje się, że jest to dość powszechna funkcja, którą można włączyć do standardowych bibliotek. Czy ktoś o jednym wie?

Moja konkretna aplikacja to tablica 2D, ale zwykle byłaby używana do znajdowania pików w FFT itp.

W szczególności w tego rodzaju problemach występuje wiele silnych szczytów, a następnie wiele mniejszych „szczytów”, które są spowodowane przez szum, który należy zignorować. To są tylko przykłady; nie moje rzeczywiste dane:

1-wymiarowe piki:

Wyjście FFT ze szczytami

2-wymiarowe piki:

Wyjście transformacji radonu z zakreślonym pikiem

Algorytm znajdowania pików znalazłby lokalizację tych pików (nie tylko ich wartości), a idealnie by znalazłby prawdziwy pik między próbkami, a nie tylko indeks z wartością maksymalną, prawdopodobnie przy użyciu interpolacji kwadratowej lub czegoś podobnego.

Zwykle zależy Ci tylko na kilku mocnych pikach, więc zostaną wybrane albo dlatego, że przekraczają pewien próg, albo dlatego, że są pierwszymi n pikami uporządkowanej listy, uszeregowanymi według amplitudy.

Jak powiedziałem, sam wiem, jak napisać coś takiego. Pytam tylko, czy istnieje wcześniej istniejąca funkcja lub pakiet, o którym wiadomo, że działa dobrze.

Aktualizacja:

I przetłumaczony skrypt MATLAB i działa przyzwoicie dla przypadku 1-D, ale mogłoby być lepiej.

Zaktualizowana aktualizacja:

sixtenbe stworzył lepszą wersję dla obudowy 1-D.


@endolith Czy masz oryginalny plik MATLAB, który przetłumaczyłeś w tym celu na Pythona? Dzięki!
Spacey



1
@endolith Wiem, że to pytanie jest dość stare, ale jest całkiem przydatne;) Dziś rano spędziłem kilka godzin find_peaks, więc dodałem tę odpowiedź, która może być przydatna w przyszłości. (Jestem pewien, że już to znalazłeś od 2009 roku, ale to dla innych ludzi + dla siebie, kiedy zadam sobie to pytanie ponownie za kilka lat!)
Basj

Odpowiedzi:


74

Funkcja scipy.signal.find_peaks, jak sama nazwa wskazuje, jest do tego przydatna. Ale ważne jest, aby dobrze zrozumieć jego parametry width, threshold, distance a przede wszystkimprominence , aby uzyskać dobry ekstrakcji szczytową.

Zgodnie z moimi testami i dokumentacją, koncepcja wyeksponowania jest „użyteczną koncepcją” do utrzymania dobrych szczytów i odrzucenia zaszumionych szczytów.

Co to jest widoczność (topograficzna) ? Jest to „minimalna wysokość niezbędna do zejścia ze szczytu na wyższy teren” , jak widać tutaj:

wprowadź opis obrazu tutaj

Idea jest taka:

Im wyższa widoczność, tym ważniejszy jest szczyt.

Test:

wprowadź opis obrazu tutaj

Celowo użyłem (hałaśliwej) sinusoidy o zmiennej częstotliwości, ponieważ wykazuje wiele trudności. Widzimy, że widthparametr nie jest tutaj zbyt przydatny, ponieważ jeśli ustawisz minimum widthza wysokie, nie będzie w stanie śledzić bardzo bliskich pików w części o wysokiej częstotliwości. Jeśli ustawisz widthzbyt nisko, będziesz miał wiele niechcianych szczytów w lewej części sygnału. Ten sam problem z distance. thresholdporównuje się tylko z bezpośrednimi sąsiadami, co nie jest tutaj przydatne. prominencejest tym, który daje najlepsze rozwiązanie. Pamiętaj, że możesz łączyć wiele z tych parametrów!

Kod:

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

x = np.sin(2*np.pi*(2**np.linspace(2,10,1000))*np.arange(1000)/48000) + np.random.normal(0, 1, 1000) * 0.15
peaks, _ = find_peaks(x, distance=20)
peaks2, _ = find_peaks(x, prominence=1)      # BEST!
peaks3, _ = find_peaks(x, width=20)
peaks4, _ = find_peaks(x, threshold=0.4)     # Required vertical distance to its direct neighbouring samples, pretty useless
plt.subplot(2, 2, 1)
plt.plot(peaks, x[peaks], "xr"); plt.plot(x); plt.legend(['distance'])
plt.subplot(2, 2, 2)
plt.plot(peaks2, x[peaks2], "ob"); plt.plot(x); plt.legend(['prominence'])
plt.subplot(2, 2, 3)
plt.plot(peaks3, x[peaks3], "vg"); plt.plot(x); plt.legend(['width'])
plt.subplot(2, 2, 4)
plt.plot(peaks4, x[peaks4], "xk"); plt.plot(x); plt.legend(['threshold'])
plt.show()

To jest to, czego szukam. Ale czy znasz jakąś implementację, która jest ważna w macierzy 2D?
Jason

43

Patrzę na podobny problem i odkryłem, że niektóre z najlepszych odniesień pochodzą z chemii (z pików znalezionych w danych masowych). Aby uzyskać dokładny przegląd algorytmów wyszukiwania szczytów, przeczytaj to . To jedna z najlepszych i najbardziej przejrzystych recenzji technik znajdowania szczytów, z jakimi się spotkałem. (Wavelety są najlepsze do znajdowania pików tego rodzaju w zaszumionych danych).

Wygląda na to, że twoje szczyty są wyraźnie określone i nie są ukryte w szumie. W takim przypadku polecam użycie gładkich pochodnych savtizky-golay do znalezienia szczytów (jeśli po prostu rozróżnisz powyższe dane, będziesz miał bałagan fałszywych alarmów). Jest to bardzo skuteczna technika i dość łatwa do wdrożenia (potrzebujesz klasy macierzy z podstawowymi operacjami). Jeśli po prostu znajdziesz przejście przez zero pierwszej pochodnej SG, myślę, że będziesz szczęśliwy.


2
Szukałem rozwiązania ogólnego przeznaczenia, a nie takiego, które działa tylko na tych konkretnych obrazach. Dostosowałem skrypt MATLAB do Pythona i działa przyzwoicie.
endolit

1
Tak jest. Matlab jest dobrym źródłem algorytmów. Jakiej techniki używa skrypt? (BTW, SG to technika bardzo ogólnego przeznaczenia).
Paul

2
Podłączyłem to powyżej. Po prostu wyszukuje lokalne maksima, które są większe niż pewien próg powyżej ich sąsiadów. Z pewnością są lepsze metody.
endolit

1
@Paul Dodałem tę stronę do zakładek. IYO i podsumowując, jaka konkretna technika, według ciebie, najlepiej sprawdziła się w tej branży zbierania szczytów?
Spacey

dlaczego zera pochodnej są lepsze niż po prostu testowanie, jeśli środek z trzech punktów jest większy lub mniejszy od pozostałych dwóch. już zastosowałem sg transfor, wydaje się, że to dodatkowy koszt.
kirill_igum

20

W Scipy jest funkcja o nazwie, scipy.signal.find_peaks_cwtktóra wydaje się być odpowiednia dla Twoich potrzeb, jednak nie mam z nią doświadczenia, więc nie mogę jej polecić.

http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks_cwt.html


12
Tak, to nie istniało, kiedy o to zapytałem, i nadal nie jestem pewien, jak go użyć
endolith

1
Dodałeś to jakiś czas temu, ale zadziałało świetnie. Korzystanie z niego jest proste. Po prostu podaj tablicę i inną tablicę (np.arange (1,10)), która zawiera listę wszystkich szerokości pików, jakie chcesz; niezłą korzyścią jest filtrowanie pod kątem wąskich lub szerokich szczytów, jeśli trzeba. Dzięki jeszcze raz!
Miles

15

Dla tych, którzy nie są pewni, które algorytmy znajdowania szczytów mają być używane w Pythonie, tutaj szybki przegląd alternatyw: https://github.com/MonsieurV/py-findpeaks

Chcąc być odpowiednikiem findpeaksfunkcji MatLab , odkryłem, że funkcja Discover_peaks Marcosa Duarte to dobry chwyt.

Całkiem łatwy w użyciu:

import numpy as np
from vector import vector, plot_peaks
from libs import detect_peaks
print('Detect peaks with minimum height and distance filters.')
indexes = detect_peaks.detect_peaks(vector, mph=7, mpd=2)
print('Peaks are: %s' % (indexes))

Co da ci:

Wykryj wyniki


1
Ponieważ ten post został napisany, find_peaksfunkcja została dodana do scipy.
onewhaleid

6

Wykrywanie pików w widmie w wiarygodny sposób było przedmiotem wielu badań, na przykład wszystkie prace nad modelowaniem sinusoidalnym dla sygnałów muzycznych / audio w latach 80-tych. Szukaj w literaturze tematu „Modelowanie sinusoidalne”.

Jeśli twoje sygnały są tak czyste, jak na przykładzie, proste „daj mi coś o amplitudzie większej niż N sąsiadów” powinno działać dość dobrze. Jeśli masz zakłócone sygnały, prostym, ale skutecznym sposobem jest spojrzenie na swoje piki w czasie, aby je śledzić: następnie wykrywasz linie widmowe zamiast szczytów widmowych. IOW, obliczasz FFT na przesuwanym oknie twojego sygnału, aby otrzymać zestaw widma w czasie (zwany również spektrogramem). Następnie przyjrzyj się ewolucji piku widmowego w czasie (tj. W kolejnych oknach).


Spójrz na szczyty w czasie? Wykrywać linie widmowe? Nie jestem pewien, co to oznacza. Czy zadziała na fale kwadratowe?
endolith

Och, mówisz o używaniu STFT zamiast FFT. To pytanie nie dotyczy konkretnie FFT; to tylko przykład. Chodzi o znalezienie pików w dowolnej ogólnej tablicy 1D lub 2D.
endolith

4

Nie sądzę, że to, czego szukasz, jest dostarczane przez SciPy. W takiej sytuacji sam napisałbym kod.

Interpolacja splajnu i wygładzanie z scipy.interpolate są całkiem przyjemne i mogą być całkiem pomocne w dopasowywaniu szczytów, a następnie znajdowaniu lokalizacji ich maksimum.


16
Przepraszam, ale myślę, że to powinien być komentarz, a nie odpowiedź. Sugeruje po prostu napisanie go samodzielnie, z niejasną sugestią dotyczącą funkcji, które mogą być przydatne (nawiasem mówiąc, te w odpowiedzi Pawła są znacznie bardziej odpowiednie).
Ami Tavory

1

Istnieją standardowe funkcje statystyczne i metody wyszukiwania wartości odstających względem danych, co prawdopodobnie jest tym, czego potrzebujesz w pierwszym przypadku. Używanie pochodnych rozwiązałoby twoją sekundę. Nie jestem jednak pewien metody, która rozwiązuje zarówno funkcje ciągłe, jak i próbkowane dane.


0

Po pierwsze, definicja „wartości szczytowej” jest niejasna, jeśli nie ma dalszych specyfikacji. Na przykład, w następującej serii, czy nazwałbyś 5-4-5 jeden szczyt czy dwa?

1-2-1-2-1-1-5-4-5-1-1-5-1

W takim przypadku będziesz potrzebować co najmniej dwóch progów: 1) tylko wysoki próg, powyżej którego może rejestrować wartość ekstremalną jako wartość szczytową; i 2) niski próg, tak aby skrajne wartości oddzielone małymi wartościami poniżej niego stały się dwoma szczytami.

Wykrywanie pików jest dobrze zbadanym tematem w literaturze dotyczącej teorii wartości ekstremalnych, znanym również jako „deklasteryzacja wartości ekstremalnych”. Jego typowe zastosowania obejmują identyfikację zdarzeń niebezpiecznych w oparciu o ciągłe odczyty zmiennych środowiskowych, np. Analizę prędkości wiatru w celu wykrycia zdarzeń burzowych.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.