Como uso um filtro Savitzky Golay para encontrar máximos locais (entre amostras) em um sinal 1D com amostragem discreta?

9

Eu tenho um sinal sísmico y (i): insira a descrição da imagem aqui

Aqui eu encontrei um máximo: i = 152,54, y = 222,29 manualmente e o plotei em vermelho.

Quero encontrar todos os máximos automaticamente.

Eu li que o Savitzky Golay Filter (SGF) pode ser usado para encontrar estimativas suavizadas tanto de um sinal quanto de seus derivados, e que um dos benefícios do SGF é que ele preserva mínimos e máximos muito melhores do que outros filtros. Isso parece ótimo para o meu uso.

Encontrei um script Matlab que gera coeficientes SGF. E usei isso para descobrir que os coeficientes SGF de 4ª ordem para o derivado. Eu codifiquei um pequeno script Matlab que

  1. encontra a derivada do sinal convolvendo o sinal com os coeficientes SGF de 4ª ordem para a derivada
  2. encontra par de amostras (i, i + 1) onde a derivada muda de sinal
  3. encontra cruzamento zero da derivada por interpolação linear entre i e i + 1

Roteiro:

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

No meu script, eu tive que testar se a derivada muda de negativa para positiva para obter a função para fornecer o resultado desejado, porém isso me confunde. A derivada para um máximo não deve passar de positiva para negativa? Existe alguma maneira melhor de distinguir entre máximos e mínimos?

Abaixo está o resultado do uso desta função para encontrar os máximos no meu sinal: insira a descrição da imagem aqui

Os resultados parecem bons, mas noto que alguns máximos não foram encontrados: i = 143,13, 190,88, 256,97.

Isso é porque eles estão perto de outros máximos?

Como posso controlar os dois máximos mais próximos?

Agradecemos antecipadamente por qualquer resposta!

Andy
fonte
Você pode plotar a saída do filtro?
Jim Clay

Respostas:

5

Embora eu não esteja familiarizado com esse tipo específico de filtro, com base no gráfico que você mostrou, eu acho que os máximos que não são encontrados pelo seu processo estão apenas contra a resolução de tempo inerente ao processo. Qualquer tipo de "suavização" implica que haja algum borrão no tempo local do sinal de interesse, de modo que, se houver dois picos próximos, é possível que eles se unam em um. É possível que um filtro de ordem inferior possa apresentar menos desse comportamento, provavelmente à custa da quantidade de suavização que você recebe.

Jason R
fonte