Como encontro picos em um conjunto de dados?

Se eu tiver um conjunto de dados que produza um gráfico como o seguinte, como determinaria algoritmicamente os valores x dos picos mostrados (neste caso, três deles):

Uma abordagem geral é suavizar os dados e encontrar picos, comparando um filtro máximo local ao suave . Em R:

argmax <- function(x, y, w=1, ...) {
  require(zoo)
  n <- length(y)
  y.smooth <- loess(y ~ x, ...)$fitted
  y.max <- rollapply(zoo(y.smooth), 2*w+1, max, align="center")
  delta <- y.max - y.smooth[-c(1:w, n+1-1:w)]
  i.max <- which(delta <= 0) + w
  list(x=x[i.max], i=i.max, y.hat=y.smooth)
}

Seu valor de retorno inclui os argumentos do máximo local ( x) - que responde à pergunta - e os índices nas matrizes x e y onde esses máximos locais ocorrem ( i).

Existem dois parâmetros a serem ajustados às circunstâncias: w é a meia largura da janela usada para calcular o máximo local. (Seu valor deve ser substancialmente menor que a metade do comprimento da matriz de dados.) Valores pequenos captarão pequenos saliências locais, enquanto valores maiores passarão sobre eles. Outro - não explícito neste código - é o spanargumento do loessmais suave. (Geralmente, está entre zero e um; reflete a largura da janela como uma proporção do intervalo de valores x.) Valores maiores suavizam os dados de maneira mais agressiva, fazendo com que os solavancos locais desapareçam completamente.

Para ver esse ajuste em vigor, vamos criar uma pequena função de teste para plotar os resultados:

test <- function(w, span) {
  peaks <- argmax(x, y, w=w, span=span)

  plot(x, y, cex=0.75, col="Gray", main=paste("w = ", w, ", span = ", span, sep=""))
  lines(x, peaks$y.hat,  lwd=2) #$
  y.min <- min(y)
  sapply(peaks$i, function(i) lines(c(x[i],x[i]), c(y.min, peaks$y.hat[i]),
         col="Red", lty=2))
  points(x[peaks$i], peaks$y.hat[peaks$i], col="Red", pch=19, cex=1.25)
}

Aqui estão algumas experiências aplicadas a alguns dados sintéticos e levemente barulhentos.

x <- 1:1000 / 100 - 5
y <- exp(abs(x)/20) * sin(2 * x + (x/5)^2) + cos(10*x) / 5 + rnorm(length(x), sd=0.05)
par(mfrow=c(3,1))
test(2, 0.05)
test(30, 0.05)
test(2, 0.2)

Uma janela larga (plotagem do meio) ou suave mais agressivo (plotagem inferior) elimina os máximos locais detectados na plotagem superior. A melhor combinação aqui é provavelmente uma janela larga e apenas suavização suave, porque a suavização agressiva parece alterar esses picos (consulte os pontos médio e direito no gráfico inferior e compare sua localização com os picos aparentes dos dados brutos). Neste exemplo, w=50e span=0.05faz um ótimo trabalho (não mostrado).

Observe que os máximos locais nos pontos de extremidade não são detectados. Estes podem ser inspecionados separadamente. (Para suportar isso, argmaxretorna os valores y suavizados.)

Essa abordagem tem várias vantagens sobre a modelagem mais formal para trabalhos de uso geral:

• Não adota nenhum modelo pré-concebido dos dados.

• Pode ser adaptado às características dos dados.

• Pode ser adaptado para detectar os tipos de picos nos quais está interessado.

Como mencionei no comentário, se a série temporal parece ser periódica, um modelo de regressão harmônica fornece uma maneira de suavizar a função e identificar o pico aplicando o primeiro e o segundo testes derivados. Huber apontou um teste não paramétrico que apresenta vantagens quando há vários picos e a função não é necessariamente periódica. Mas não há almoço grátis. Embora haja vantagens no método que ele menciona, pode haver desvantagens se um modelo paramétrico for apropriado. Esse é sempre o outro lado do uso de técnicas não paramétricas. Embora evite suposições paramétricas, a abordagem paramétrica é melhor quando as suposições paramétricas são apropriadas. Seu procedimento também não tira o máximo proveito da estrutura de séries temporais nos dados.

Penso que, embora seja apropriado apontar vantagens de um procedimento sugerido, também é importante apontar as possíveis desvantagens. Tanto minha abordagem quanto a de Huber encontram os picos de maneira eficiente. No entanto, acho que o procedimento dele exige um pouco mais de trabalho quando um máximo local é menor que o pico mais alto determinado anteriormente.

Uma abordagem clássica de detecção de pico no processamento de sinal é a seguinte:

1. Filtre o sinal para uma faixa razoável razoável, dependendo da taxa de amostragem e das propriedades do sinal, por exemplo, para ECG, um filtro passa-banda IIR a 0,5-20Hz, um filtro de fase zero garantirá que nenhuma mudança de fase (e atraso de tempo associado) seja introduzida
2. Uma transformação de hilbert ou uma abordagem wavelet pode então ser usada para enfatizar os picos
3. Um limite estático ou dinâmico pode ser aplicado, onde todas as amostras acima do limite são consideradas picos. No caso de um limite dinâmico, geralmente é definido como um desvio padrão do limite N acima ou abaixo de uma estimativa da média móvel da média.

Outra abordagem que funciona é comparar um sinal filtrado de alta precisão com um filtro altamente suavizado (filtro passa-baixo ou mediano) e aplicar a etapa 3.

Espero que isto ajude.