¿Cómo encuentro picos en un conjunto de datos?

Si tengo un conjunto de datos que produce un gráfico como el siguiente, ¿cómo determinaría algorítmicamente los valores de x de los picos mostrados (en este caso, tres de ellos):

Un enfoque general es suavizar los datos y luego encontrar picos comparando un filtro máximo local con el suavizado . En 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)
}

Su valor de retorno incluye los argumentos de los máximos locales ( x), que responden a la pregunta, y los índices en las matrices x e y donde ocurren esos máximos locales ( i).

Hay dos parámetros para ajustarse a las circunstancias: w es el ancho medio de la ventana utilizada para calcular el máximo local. (Su valor debe ser sustancialmente menor que la mitad de la longitud de la matriz de datos). Los valores pequeños recogerán pequeñas protuberancias locales, mientras que los valores más grandes pasarán por encima de ellos. Otro, no explícito en este código, es el spanargumento del loesssuavizador. (Por lo general, está entre cero y uno; refleja el ancho de una ventana como una proporción del rango de valores de x). Los valores más grandes suavizarán los datos de manera más agresiva, haciendo desaparecer por completo los baches locales.

Para ver este ajuste en efecto, creemos una pequeña función de prueba para trazar los 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)
}

Aquí hay algunos experimentos aplicados a algunos datos sintéticos, un poco ruidosos.

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)

Ya sea una ventana ancha (gráfico central) o una suavidad más agresiva (gráfico inferior) eliminan los máximos locales detectados en el gráfico superior. La mejor combinación aquí es probablemente una ventana amplia y solo un suavizado suave, porque el suavizado agresivo parece cambiar estos picos (vea los puntos medio y derecho en el gráfico inferior y compare sus ubicaciones con los picos aparentes de los datos sin procesar). En este ejemplo, w=50y span=0.05hace un gran trabajo (no se muestra).

Observe que los máximos locales en los puntos finales no se detectan. Estos pueden ser inspeccionados por separado. (Para respaldar esto, argmaxdevuelve los valores de y suavizados).

Este enfoque tiene varias ventajas sobre el modelado más formal para el trabajo de propósito general:

• No adopta ningún modelo preconcebido de los datos.

• Se puede adaptar a las características de los datos.

• Se puede adaptar para detectar los tipos de picos que le interesan.

Como mencioné en el comentario, si la serie temporal parece ajustarse periódicamente, un modelo de regresión armónica proporciona una forma de suavizar la función e identificar el pico aplicando las pruebas derivadas primera y segunda. Huber ha señalado una prueba no paramétrica que tiene ventajas cuando hay múltiples picos y la función no es necesariamente periódica. Pero no hay almuerzo gratis. Si bien existen las ventajas de su método que menciona, puede haber desventajas si un modelo paramétrico es apropiado. Esa es siempre la otra cara del uso de técnicas no paramétricas. Aunque evita suposiciones paramétricas, el enfoque paramétrico es mejor cuando las suposiciones paramétricas son apropiadas. Su procedimiento tampoco aprovecha al máximo la estructura de series temporales en los datos.

Creo que si bien es apropiado señalar las ventajas de un procedimiento sugerido, también es importante señalar las posibles desventajas. Tanto mi enfoque como el de Huber encuentran los picos de manera eficiente. Sin embargo, creo que su procedimiento requiere un poco más de trabajo cuando un máximo local es más bajo que el pico más alto previamente determinado.

Un enfoque clásico de detección de picos en el procesamiento de señales es el siguiente:

1. Filtre la señal a un rango razonable razonable, dependiendo de la frecuencia de muestreo y las propiedades de la señal, por ejemplo, para ECG, un filtro de paso de banda IIR @ 0.5-20Hz, un filtro de fase cero asegurará que no se introduzca un cambio de fase (y el retraso de tiempo asociado)
2. Una transformación hilbert o un enfoque wavelet se puede utilizar para enfatizar los picos
3. Entonces se puede aplicar un umbral estático o dinámico, donde todas las muestras por encima del umbral se consideran picos. En el caso de un umbral dinámico, generalmente se define como un umbral N desviaciones estándar por encima o por debajo de una estimación promedio móvil de la media.

Otro enfoque que funciona es comparar una señal filtrada de paso alto agudo con una señal altamente suavizada (paso bajo o filtrado medio) y aplicar el paso 3.

Espero que esto ayude.