Algorytm znajdowania pików dla Python / SciPy

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:

2-wymiarowe piki:

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.

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:

Idea jest taka:

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

Test:

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

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.

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

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:

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).

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.

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.

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.