Jak użyć filtra Savitzky Golay, aby znaleźć lokalne maksima (pomiędzy próbkami) w dyskretnie próbkowanym sygnale 1D?

Mam sygnał sejsmiczny y (i):

Tutaj znalazłem jedno maksimum: i = 152,54, y = 222,29 ręcznie i narysowałem je na czerwono.

Chcę znaleźć wszystkie maksima automatycznie.

Czytałem, że filtr Savitzky Golay (SGF) może być wykorzystywany do wyszukiwania wygładzonych oszacowań zarówno sygnału, jak i jego pochodnych, oraz że jedną z korzyści SGF jest to, że zachowuje minima i maksima znacznie lepiej niż inne filtry. Brzmi świetnie na mój użytek.

Znalazłem skrypt Matlab, który generuje współczynniki SGF. I wykorzystał to do stwierdzenia, że ​​współczynniki SGF czwartego rzędu dla pochodnej. Kodowałem to mały skrypt Matlaba

1. znajduje pochodną sygnału przez zwoje sygnału ze współczynnikami SGF czwartego rzędu dla pochodnej
2. znajduje parę próbek (i, i + 1), w których znak pochodnej zmienia się
3. znajduje przecięcie zera pochodnej przez interpolację liniową między i i + 1

Scenariusz:

function [maxX,maxY] = findLocalMax(y)
% Kernel for 4th order Savitzky-Golay filter for finding derivative:
d4 = [0.0724 -0.1195 -0.1625 -0.1061 0 0.1061 0.1625 0.1195 -0.0724];

dy = conv(y,d4,'same'); % derivative

[m n] = size(dy);
maxX = [];
maxY = [];
for i = 1 : n - 1
  if dy(i) < 0 && dy(i+1) > 0 % max somewhere between i and i+1
    a = dy(i)/(dy(i) - dy(i+1)); % linear interpolation
    mx = i + a;
    maxX = [maxX mx];
    my = y(i)*(1-a) + y(i+1)*a; % linear interpolation
    maxY = [maxY my];
  end
end

W moim skrypcie musiałem sprawdzić, czy pochodna zmienia się z ujemnej na dodatnią, aby uzyskać funkcję zapewniającą pożądany wynik, jednak to mnie myli. Czy pochodna dla maksimum nie powinna zmieniać się z dodatniej na ujemną? Czy istnieje lepszy sposób na odróżnienie maksimów od minimów?

Poniżej znajduje się wynik użycia tej funkcji do znalezienia maksimów na moim sygnale:

Wyniki wyglądają dobrze, ale zauważam, że nie znaleziono niektórych maksimów: i = 143,13, 190,88, 256,97.

Czy to dlatego, że mają zbliżyć się do innych maksimów?

Jak mogę kontrolować najbliższe dwa maksima?

Z góry dziękuję za wszelkie odpowiedzi!

Chociaż nie jestem zaznajomiony z tym konkretnym typem filtra, w oparciu o pokazany przez ciebie wykres, zgaduję, że maksima, które nie są znalezione w twoim procesie, po prostu zbliżają się do rozdzielczości czasowej właściwej dla procesu. Każdy rodzaj „wygładzania” oznacza, że ​​występuje pewne lokalne rozmazanie sygnału zainteresowania, tak że jeśli w pobliżu znajdują się dwa szczyty, możliwe jest, że zostaną one połączone w jeden. Możliwe jest, że filtr niższego rzędu może wykazywać mniej tego zachowania, prawdopodobnie kosztem wygładzenia, które otrzymujesz.