¿Cómo uso un filtro Savitzky Golay para encontrar máximos locales (entre muestras) en una señal 1D discretamente muestreada?

9

Tengo una señal sísmica y (i): ingrese la descripción de la imagen aquí

Aquí he encontrado un máximo: i = 152.54, y = 222.29 manualmente y lo tracé en rojo.

Quiero encontrar todos los máximos automáticamente.

Leí que el filtro Savitzky Golay (SGF) se puede utilizar para encontrar estimaciones uniformes de una señal y sus derivados, y que uno de los beneficios del SGF es que conserva los mínimos y máximos mucho mejor que otros filtros. Esto suena genial para mi uso.

Encontré un script de Matlab que genera coeficientes SGF. Y usé esto para encontrar que los coeficientes SGF de cuarto orden para la derivada. Codifiqué un pequeño script de Matlab que

  1. encuentra la derivada de la señal haciendo girar la señal con los coeficientes SGF de cuarto orden para la derivada
  2. encuentra un par de muestras (i, i + 1) donde la derivada cambia de signo
  3. encuentra el cruce cero de la derivada por interpolación lineal entre i e i + 1

Guión:

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

En mi script, tuve que probar si la derivada cambia de negativa a positiva para que la función dé el resultado deseado, sin embargo, esto me confunde. ¿No debería la derivada para un máximo pasar de positivo a negativo? ¿Hay alguna forma mejor de distinguir entre máximos y mínimos?

A continuación se muestra el resultado de usar esta función para encontrar los máximos en mi señal: ingrese la descripción de la imagen aquí

Los resultados se ven bien, pero noto que no se encuentran algunos máximos: i = 143.13, 190.88, 256.97.

¿Es esto porque están cerca de otros máximos?

¿Cómo puedo controlar los dos máximos más cercanos?

Gracias de antemano por cualquier respuesta!

smoothing interpolation peak-detection Andy
fuente
¿Puedes trazar la salida del filtro?
Jim Clay

Respuestas:

5

Si bien no estoy familiarizado con este tipo específico de filtro, según el gráfico que ha mostrado, supongo que los máximos que no se encuentran en su proceso solo están en contra de la resolución de tiempo inherente en el proceso. Cualquier tipo de "suavizado" implica que hay un poco de tiempo local de la señal de interés, de modo que si hay dos picos cercanos, es posible que se fusionen en uno. Es posible que un filtro de orden inferior presente menos comportamiento, probablemente a expensas de la cantidad de suavizado que obtenga.

Jason R
fuente