Cómo realizar un suavizado de datos con el filtro Savitzky-Golay y R

November 11, 2017

El filtro Savitzky-Golay se utiliza normalmente para el suavizado de datos, ajustando un polinomio con al menos k+1 puntos equidistantes para determinar el nuevo valor de cada punto.

En teledetección, cuando se trabaja con firmas espectrales se suele utilizar este filtro por conservar las características de la distribución inicial, los máximos y mínimos relativos y el ancho entre picos. El filtro simula una ventana que se desliza por los datos, suavizándolos a partir de los datos de la ventana.

Para poder emplearlo en R, necesitaremos el paquete signal, el cual contiene funciones de procesamiento compatibles con Matlab/Octvave, como filtros, funciones de filtrado, visualización de modelos de filtros, interpolación, etc. Las funciones que se van a emplear son:

  • sgolay, para calcular los coeficientes.
  • filter, para aplicar los coeficientes a los datos, con un ajuste de mínimos cuadrados (cuadrado de las diferencias entre los valores calculados y observados).

Cómo establecer los parámetros de filtrado

  • p es el orden del polinomio. Si es alto, gestiona mejor los rasgos espectrales estrechos, pero suaviza menos los anchos.
  • n es la longitud del filtro (siempre impar). Entre 1 y 2 la anchura del rasgo m?s estrecho en los datos.
  • m devuelve la derivada de los coeficientes de filtrado
  • ts es el factor de escala de tiempo

Establecer el directorio de trabajo

# Directorio de trabajo
dir <- setwd("/home/joan/TRATAMIENTO/")

Cargar los datos y crear un DF

data <- read.csv("FIRMAS_REF.csv", header = T, sep = ",", dec = ",")
# Nótese que se transforman los vectores a factores, de lo contrario el filtro no funcionará
df <- data.frame(data, stringsAsFactors = TRUE)

Visualizamos las gráficas para determinar la n

plot(df$WL, df$Firma_1, type = "l", main = "Gr?fica firma 1", xlab = "wavelength", ylab = "Reflectividad")

plot

Se escoge el ancho de entre la longitud de onda 737 - 748 como rasgo más estrecho, de manera que su doble nos da una n de 21.

#install.packages("signal")
library(signal)

# Calcular los coeficientes de filtro
p = 3
n = 21
m = 0
sg <- sgolay(p, n, m)
# Se aplica el filtro con los coeficientes a nuestros datos
library(signal)

df$Firma_1sg <- filter(sg, df$Firma_1)
# Visualizamos el resultado
library(ggplot2)
plotsg <- ggplot(df) +
  geom_point(aes(df$WL, df$Firma_1), size = 0.1) +
  geom_line(aes(df$WL, df$Firma_1sg), col = "red2", size = 1) +
  xlab("Wavelength") +
  ylab("Reflectancia") +
  theme_light()
print(plotsg)

golay

# Exportamos como DF
dfsg <- data.frame(df$WL, df$Firma_1sg)
write.csv(dfsg, file = "firmas_SG.csv")