10.4: Uso de R para limpiar datos

R tiene dos funciones útiles, filter () y fft (), que podemos usar para suavizar o filtrar el ruido y para eliminar las señales de fondo. Para explorar su uso, primero vamos a crear dos conjuntos de datos que podemos usar como ejemplos: una señal ruidosa y una señal pura superpuesta sobre un fondo exponencial. Para crear la señal ruidosa, primero creamos un vector de 256 valores que define el eje x; aunque no especificaremos una unidad aquí, estos podrían ser tiempos o frecuencias. A continuación usamos la función dnorm () de R para generar una señal gaussiana pura con una media de 125 y una desviación estándar de 10, y la función rnorm () de R para generar 256 puntos de ruido aleatorio con una media de cero y una desviación estándar de 10. Finalmente, agregamos la señal pura y el ruido para llegar a nuestra señal ruidosa y luego trazamos la señal ruidosa y superponemos la señal pura.

x = seq (1,256,1)
gaus_signal = 1250 * dnorm (x, media = 125, sd = 10)
ruido = rnorm (256, media = 0, sd = 10)
noisy_signal = gaus_signal +
gráfico de ruido (x = x, y = noisy_signal, type = “l”, lwd = 2, col = “azul”, xlab = “x”, ylab = “señal”)
líneas (x = x, y = gaus_signal, lwd = 2)

Figura\(\PageIndex{1}\): Ejemplo de una señal ruidosa con una relación señal/ruido de 5.1. La figura\(\PageIndex{3}\)\(\PageIndex{4}\), la figura y las figuras\(\PageIndex{8}\) muestran esta misma figura después de aplicar un filtro de promedio móvil de siete puntos, un filtro de suavizado Savitizky-Golay de siete puntos y un filtro de Fourier.

Para estimar la relación señal/ruido, utilizamos el máximo de la señal pura y la desviación estándar de la señal ruidosa determinada usando 100 puntos divididos uniformemente entre los dos extremos.

s_to_n = max (gaus_signal) /sd (noisy_signal [c (1:50, 201:250)])
s_to_n

[1] 5.14663

Para crear una señal superpuesta sobre un fondo exponencial, utilizamos la función exp () de R para generar 256 puntos para la señal del fondo, agregarlo a nuestra señal gaussiana pura y trazar el resultado.

exp_bkgd = 30*exp (-0.01 * x)
plot (x, exp_bkgd, type = “l”)
signal_bkgd = gaus_signal + exp_bkgd
plot (x = x, y = signal_bkgd, type = “l”, lwd = 2, col = “azul”, xlab = “x”, ylab = “señal”, ylim = c (0,60))
líneas (x = x, y = gaus_señal, lwd = 2, lty = 2)

Figura\(\PageIndex{2}\): Ejemplo de una señal gaussiana pura superpuesta sobre un fondo exponencial. La figura\(\PageIndex{5}\) muestra esta misma figura después de usar un filtro Savitzky-Golay derivado de siete puntos para eliminar el fondo.

Uso de R filter () Función para suavizar el ruido y eliminar las señales de fondo

La función filter () de R toma la forma general

filtro (x, filtro)

donde x es el objeto que se está filtrando y filter es un objeto que contiene los coeficientes del filtro. Para crear un filtro de promedio móvil de siete puntos, usamos la función rep () para crear un vector que tiene siete valores idénticos, cada uno igual a 1/7.

mov_avg_7 = rep (1/7, 7)

Aplicar este filtro a nuestra señal ruidosa devuelve el siguiente resultado

noisy_signal_movavg = filter (noisy_signal, mov_avg_7)
plot (x = x, y = noisy_signal_movavg, type = “l”, lwd = 2, col = “azul”, xlab = “x”, ylab = “señal”)
líneas (x = x, y = gaus_signal, lwd = 2)

con la relación señal/ruido mejorada para

s_to_n_movavg = max (gaus_signal) /sd (noisy_signal_movavg [c (1:50, 200:250)], na.rm = VERDADERO)
s_to_n_movavg

[1] 11.29943

Tenga en cuenta que debemos agregar na.rm = VERDADERO a la función sd () porque aplicar un filtro de promedio móvil de siete puntos reemplaza los tres primeros y los últimos tres puntos con valores de NA que debemos decirle a la función sd () que ignore.

Figura\(\PageIndex{3}\): El resultado de usar la función filter () de R para aplicar un filtro de promedio móvil de siete puntos a la señal ruidosa de la Figura\(\PageIndex{1}\).

Para crear un filtro de suavizado Savitzky-Golay de siete puntos, creamos un vector para almacenar los coeficientes, obteniendo los valores del papel original (Savitzky, A.; Golay, M. J. E. Anal Chem, 1964, 36, 1627-1639) y luego aplicarlo a nuestro ruidoso señal, obteniendo los resultados a continuación.

sg_smooth_7 = c (-2,3,6,7,6,5, -2) /21
noisy_signal_sg = filter (noisy_signal, sg_smooth_7)
plot (x = x, y = noisy_signal_sg, type = “l”, lwd = 2, col = “azul”, xlab = “x”, ylab = “señal”)
líneas (x = y, = gaus_signal, lwd = 2)
s_to_n_movavg = max ( gaus_signal) /sd (noisy_signal_sg [c (1:50, 200:250)], na.rm = VERDADERO)
s_to_n_movavg

[1] 7.177931

Figura\(\PageIndex{4}\): El resultado de usar la función filter () de R para aplicar un filtro de smooting Savitizky-Golay de siete puntos a la señal ruidosa de la Figura\(\PageIndex{1}\).

Para eliminar un fondo de una señal, utilizamos el mismo enfoque, sustituyendo un filtro Savitxky-Golay de primera derivada (o de orden superior).

sg_fd_7 = c (22, -67, -58, 0, 58, 67, -22) /252
signal_bkgd_sg = filter (signal_bkgd, sg_fd_7)
plot (x = x, y = signal_bkgd_sg, type = “l”, lwd = 2, col = “azul”, xlab = “x”, ylab = “señal”)

Figura\(\PageIndex{5}\): El resultado de usar la función filter () de R para aplicar un filtro Savitzky-Golay de siete puntos de primera derivada a la señal ruidosa de la Figura\(\PageIndex{1}\).

Uso de la función fft () de R para el filtrado de Fourier

Para completar una transformada de Fourier en R utilizamos la función fft (), que toma la forma fft (z, inverse = FALSE) donde z es el objeto que contiene los valores a los que deseamos aplicar la transformada de Fourier y donde establecer inverse = TRUE permite una inversa Transformada de Fourier. Antes de aplicar el filtrado de Fourier a nuestra señal ruidosa, primero apliquemos la función fft () a un vector que contenga los enteros del 1 al 8. Primero creamos un vector para mantener nuestros valores y aplicar la función fft () al vector, obteniendo los siguientes resultados

vector_prueba = seq (1, 8, 1)
test_vector_ft = fft (test_vector)
test_vector_ft

[1] 36+0.000000i -4+9.656854i -4+4.000000i -4+1.656854i -4+0.000000i -4-1.656854i

[7] -4-4.000000i -4-9.656854i

Cada uno de los ocho resultados es un número complejo con un componente real y uno imaginario. Tenga en cuenta que el componente real del primer valor es 36, que es la suma de los elementos en nuestro vector de prueba. Obsérvese, también, la simetría en los valores restantes donde los valores segundo y octavo, los valores tercero y séptimo, y los valores cuarto y sexto son idénticos a excepción de un cambio de signo para el componente imaginario.

Tomando la transformada inversa de Fourier devuelve los ocho valores originales (tenga en cuenta que los términos imaginarios son ahora cero), pero cada uno es ocho veces mayor en valor que en nuestro vector original.

test_vector_ifft = fft (test_vector_ft, inverso = VERDADERO)
test_vector_ifft

[1] 8+0i 16-0i 24+0i 32+0i 40+0i 48+0i 56-0i 64+0i

Para compensar esto, dividimos por la longitud de nuestro vector

test_vector_ifft = fft (test_vector_ft, inverso = VERDADERO) /longitud (test_vector)
test_vector_ifft

[1] 1+0i 2-0i 3+0i 4+0i 5+0i 6+0i 7-0i 8+0i

que devuelve nuestro vector original.

Con este fondo en su lugar, usemos R para completar un filtrado de Fourier de nuestra señal ruidosa. Primero, completamos la transformada de Fourier de la señal ruidosa y examinamos los valores para el componente real, usando la función Re () de R para extraerlos. Debido a la simetría señalada anteriormente, solo necesitamos mirar la primera mitad de los componentes reales (x = 1 a x = 128).

noisy_signal_ft = fft (noisy_signal)
gráfico (x = x [1:128], y = Re (noisy_signal_ft) [1:128], type = “l”, col = “azul”, xlab = “”, ylab = “intensidad”, lwd = 2)

Figura\(\PageIndex{6}\): Gráfica que muestra los primeros 128 componentes reales de los datos de la Figura\(\PageIndex{1}\) después de completar una transformada de Fourier.

A continuación, buscamos dónde la magnitud de la señal ha decaído a lo que parece ser ruido aleatorio y establecemos estos valores a cero. En este ejemplo, conservamos los primeros 24 puntos (y los últimos 24 puntos; recordemos la simetría señalada anteriormente) y establecemos tanto el componente real como el imaginario en 0 + 0i.

noisy_signal_ft [25:232] = 0 + 0i
parcela (x = x, y = Re (noisy_signal_ft), tipo = “l”, col = “azul”, xlab = “”, ylab = “intensidad”, lwd = 2)

Figura\(\PageIndex{7}\): Gráfica que muestra los componentes reales de los datos de la Figura\(\PageIndex{1}\) después de completar una transformada de Fourier y poner a cero todos los valores distintos de los primeros 24 y los últimos 24 puntos. Aquí estamos asumiendo que estos puntos están dominados por la señal original, mientras que los puntos restantes son en su mayoría del ruido.

Finalmente, tomamos la transformada inversa de Fourier y mostramos la señal filtrada resultante y reportamos la relación señal/ruido.

noisy_signal_ifft = fft (noisy_signal_ft, inverso = VERDADERO) /longitud (noisy_signal_ft)
gráfico (x = x, y = Re (noisy_signal_ifft), tipo = “l”, col = “azul”, xlab = “”, ylab = “intensidad”, ylim = c (-20,60), lwd = 3)
líneas (x = x, y = gaus_signal, lwd =2, col = “negro”)

s_to_n = 50/sd (Re (noisy_signal_ifft) [c (1:50, 200:250)], na.rm = VERDADERO)
s_to_n

[1] 9.695329

Figura\(\PageIndex{8}\): El resultado final de Fourier filtrando los datos de la Figura\(\PageIndex{1}\).