r signal-processing frequency-analysis

¿Cómo ejecuto un filtro de paso alto o paso bajo en los puntos de datos en R?



signal-processing frequency-analysis (8)

Soy un principiante en R y he tratado de encontrar información sobre lo siguiente sin encontrar nada.

El gráfico verde en la imagen está compuesto por los gráficos rojo y amarillo. Pero digamos que solo tengo los puntos de datos de algo así como el gráfico verde. ¿Cómo extraigo las frecuencias bajas / altas (es decir, aproximadamente los gráficos rojos / amarillos) usando un filtro de paso bajo / paso alto ?

Actualización: el gráfico se generó con

number_of_cycles = 2 max_y = 40 x = 1:500 a = number_of_cycles * 2*pi/length(x) y = max_y * sin(x*a) noise1 = max_y * 1/10 * sin(x*a*10) plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5)) points(x, y + noise1, col="green", pch=20) points(x, noise1, col="yellow", pch=20)

Actualización 2: Utilizando el filtro Butterworth en el paquete de signal sugerido, recibo lo siguiente:

library(signal) bf <- butter(2, 1/50, type="low") b <- filter(bf, y+noise1) points(x, b, col="black", pch=20) bf <- butter(2, 1/25, type="high") b <- filter(bf, y+noise1) points(x, b, col="black", pch=20)

Los cálculos fueron un poco laboriosos, signal.pdf dio al lado de ningún indicio sobre qué valores W debería tener, pero la documentación original de octavas al menos mencionó los radians que me ayudaron. Los valores en mi gráfica original no fueron elegidos con ninguna frecuencia específica en mente, así que terminé con las siguientes frecuencias no tan simples: f_low = 1/500 * 2 = 1/250 , f_high = 1/500 * 2*10 = 1/25 y la frecuencia de muestreo f_s = 500/500 = 1 . Luego elegí un f_c en algún lugar entre las frecuencias baja y alta para los filtros de paso bajo / alto (1/100 y 1/50, respectivamente).

Me encontré con un problema similar recientemente y no encontré las respuestas aquí particularmente útiles. Aquí hay un enfoque alternativo.

Comencemos definiendo los datos de ejemplo de la pregunta:

number_of_cycles = 2 max_y = 40 x = 1:500 a = number_of_cycles * 2*pi/length(x) y = max_y * sin(x*a) noise1 = max_y * 1/10 * sin(x*a*10) y <- y + noise1 plot(x, y, type="l", ylim=range(-1.5*max_y,1.5*max_y,5), lwd = 5, col = "green")

Entonces, la línea verde es el conjunto de datos que queremos para el filtro de paso bajo y paso alto.

Nota al margen: La línea en este caso podría expresarse como una función mediante spline cúbico ( spline(x,y, n = length(x)) ), pero con datos del mundo real esto raramente sería el caso, así que supongamos que no es posible expresar el conjunto de datos como una función.

La forma más sencilla de suavizar los datos que he encontrado es usar loess o smooth.spline con span / spar apropiado. De acuerdo con los estadísticos, loess / smooth.spline probablemente no sea el enfoque correcto aquí , ya que en realidad no presenta un modelo definido de los datos en ese sentido. Una alternativa es usar modelos aditivos generalizados gam() función gam() del paquete mgcv ). Mi argumento para usar loess o suavizado spline aquí es que es más fácil y no hace la diferencia, ya que estamos interesados ​​en el patrón visible resultante. Los conjuntos de datos del mundo real son más complicados que en este ejemplo y encontrar una función definida para filtrar varios conjuntos de datos similares puede ser difícil. Si el ajuste visible es bueno, ¿por qué hacerlo más complicado con los valores de R2 y p? Para mí, la aplicación es visual para la cual loess / suavizado splines son métodos apropiados. Ambos métodos asumen relaciones polinomiales con la diferencia de que el loess es más flexible también utilizando polinomios de mayor grado, mientras que el spline cúbico es siempre cúbico (x ^ 2). Cuál usar depende de las tendencias en un conjunto de datos. Dicho esto, el siguiente paso es aplicar un filtro de paso bajo en el conjunto de datos utilizando loess() o smooth.spline() :

lowpass.spline <- smooth.spline(x,y, spar = 0.6) ## Control spar for amount of smoothing lowpass.loess <- loess(y ~ x, data = data.frame(x = x, y = y), span = 0.3) ## control span to define the amount of smoothing lines(predict(lowpass.spline, x), col = "red", lwd = 2) lines(predict(lowpass.loess, x), col = "blue", lwd = 2)

La línea roja es el filtro spline suavizado y el filtro azul loess. Como ve, los resultados difieren ligeramente. Supongo que uno de los argumentos para usar GAM sería encontrar la mejor opción, si las tendencias fueran así de claras y constantes entre los conjuntos de datos, pero para esta aplicación estos dos ajustes son lo suficientemente buenos para mí.

Después de encontrar un filtro de paso bajo apropiado, el filtro de paso alto es tan simple como restar los valores de filtro de paso bajo de y :

highpass <- y - predict(lowpass.loess, x) lines(x, highpass, lwd = 2)

Esta respuesta llega tarde, pero espero que ayude a alguien más a lidiar con un problema similar.

No estoy seguro de si algún filtro es la mejor manera para ti. Un instrumento más útil para ese objetivo es la transformación rápida de Fourier.

Por solicitud de OP:

El cran.r-project.org/web/packages/signal/signal.pdf contiene todo tipo de filtros para el procesamiento de señales. La mayor parte es comparable / compatible con las funciones de procesamiento de señal en Matlab / Octave.

También tuve problemas para descubrir cómo el parámetro W en la función de mantequilla se asigna al corte del filtro, en parte porque la documentación para filter y filtfilt es incorrecta a partir de la publicación (sugiere que W = .1 resultaría en un 10 Filtro Hz lp cuando se combina con filtfilt cuando la frecuencia de muestreo de la señal es Fs = 100, pero en realidad, es solo un filtro de 5 Hz lp: el corte de media amplitud es de 5 Hz cuando se usa filtfilt, pero el corte de media potencia es 5 Hz cuando solo aplica el filtro una vez, usando la función de filtro). Estoy publicando un código de demostración que escribí a continuación que me ayudó a confirmar cómo funciona todo esto, y que podría utilizar para verificar que un filtro esté haciendo lo que desea.

#Example usage of butter, filter, and filtfilt functions #adapted from https://rdrr.io/cran/signal/man/filtfilt.html library(signal) Fs <- 100; #sampling rate bf <- butter(3, 0.1); #when apply twice with filtfilt, #results in a 0 phase shift #5 Hz half-amplitude cut-off LP filter # #W * (Fs/2) == half-amplitude cut-off when combined with filtfilt # #when apply only one time, using the filter function (non-zero phase shift), #W * (Fs/2) == half-power cut-off t <- seq(0, .99, len = 100) # 1 second sample #generate a 5 Hz sine wave x <- sin(2*pi*t*5) #filter it with filtfilt y <- filtfilt(bf, x) #filter it with filter z <- filter(bf, x) #plot original and filtered signals plot(t, x, type=''l'') lines(t, y, col="red") lines(t,z,col="blue") #estimate signal attenuation (proportional reduction in signal amplitude) 1 - mean(abs(range(y[t > .2 & t < .8]))) #~50% attenuation at 5 Hz using filtfilt 1 - mean(abs(range(z[t > .2 & t < .8]))) #~30% attenuation at 5 Hz using filter #demonstration that half-amplitude cut-off is 6 Hz when apply filter only once x6hz <- sin(2*pi*t*6) z6hz <- filter(bf, x6hz) 1 - mean(abs(range(z6hz[t > .2 & t < .8]))) #~50% attenuation at 6 Hz using filter #plot the filter attenuation profile (for when apply one time, as with "filter" function): hf <- freqz(bf, Fs = Fs); plot(c(0, 20, 20, 0, 0), c(0, 0, 1, 1, 0), type = "l", xlab = "Frequency (Hz)", ylab = "Attenuation (abs)") lines(hf$f[hf$f<=20], abs(hf$h)[hf$f<=20]) plot(c(0, 20, 20, 0, 0), c(0, 0, -50, -50, 0), type = "l", xlab = "Frequency (Hz)", ylab = "Attenuation (dB)") lines(hf$f[hf$f<=20], 20*log10(abs(hf$h))[hf$f<=20]) hf$f[which(abs(hf$h) - .5 < .001)[1]] #half-amplitude cutoff, around 6 Hz hf$f[which(20*log10(abs(hf$h))+6 < .2)[1]] #half-amplitude cutoff, around 6 Hz hf$f[which(20*log10(abs(hf$h))+3 < .2)[1]] #half-power cutoff, around 5 Hz

Un método es usar la fast fourier transform implementada en R como fft . Aquí hay un ejemplo de un filtro de paso alto. De las tramas anteriores, la idea implementada en este ejemplo es obtener la serie en amarillo a partir de la serie en verde (sus datos reales).

# I''ve changed the data a bit so it''s easier to see in the plots par(mfrow = c(1, 1)) number_of_cycles = 2 max_y = 40 N <- 256 x = 0:(N-1) a = number_of_cycles * 2 * pi/length(x) y = max_y * sin(x*a) noise1 = max_y * 1/10 * sin(x*a*10) plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5)) points(x, y + noise1, col="green", pch=20) points(x, noise1, col="yellow", pch=20) ### Apply the fft to the noisy data y_noise = y + noise1 fft.y_noise = fft(y_noise) # Plot the series and spectrum par(mfrow = c(1, 2)) plot(x, y_noise, type=''l'', main=''original serie'', col=''green4'') plot(Mod(fft.y_noise), type=''l'', main=''Raw serie - fft spectrum'')

### The following code removes the first spike in the spectrum ### This would be the high pass filter inx_filter = 15 FDfilter = rep(1, N) FDfilter[1:inx_filter] = 0 FDfilter[(N-inx_filter):N] = 0 fft.y_noise_filtered = FDfilter * fft.y_noise

par(mfrow = c(2, 1)) plot(x, noise1, type=''l'', main=''original noise'') plot(x, y=Re( fft( fft.y_noise_filtered, inverse=TRUE) / N ) , type=''l'', main = ''filtered noise'')

Utilice la función filtfilt en lugar del filtro (señal del paquete) para deshacerse del cambio de señal.

library(signal) bf <- butter(2, 1/50, type="low") b1 <- filtfilt(bf, y+noise1) points(x, b1, col="red", pch=20)

Vea este enlace donde hay un código R para filtrar (señales médicas). Es de Matt Shotwell y el sitio está lleno de información R / stats interesante con una inclinación médica:

biostattmat.com

El paquete fftfilt contiene muchos algoritmos de filtrado que también deberían ser de ayuda.

hay un paquete en CRAN llamado FastICA , este calcula la aproximación de las señales de fuente independientes; sin embargo, para calcular ambas señales necesita una matriz de al menos 2xn observaciones mixtas (para este ejemplo), este algoritmo no puede determinar las dos señales independientes con solo 1xn vector. Vea el ejemplo a continuación. Espero que esto le pueda ayudar.

number_of_cycles = 2 max_y = 40 x = 1:500 a = number_of_cycles * 2*pi/length(x) y = max_y * sin(x*a) noise1 = max_y * 1/10 * sin(x*a*10) plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5)) points(x, y + noise1, col="green", pch=20) points(x, noise1, col="yellow", pch=20) ###################################################### library(fastICA) S <- cbind(y,noise1)#Assuming that "y" source1 and "noise1" is source2 A <- matrix(c(0.291, 0.6557, -0.5439, 0.5572), 2, 2) #This is a mixing matrix X <- S %*% A a <- fastICA(X, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1, method = "R", row.norm = FALSE, maxit = 200, tol = 0.0001, verbose = TRUE) par(mfcol = c(2, 3)) plot(S[,1 ], type = "l", main = "Original Signals", xlab = "", ylab = "") plot(S[,2 ], type = "l", xlab = "", ylab = "") plot(X[,1 ], type = "l", main = "Mixed Signals", xlab = "", ylab = "") plot(X[,2 ], type = "l", xlab = "", ylab = "") plot(a$S[,1 ], type = "l", main = "ICA source estimates", xlab = "", ylab = "") plot(a$S[, 2], type = "l", xlab = "", ylab = "")