parámetros para filtro de abeto de paso bajo usando scipy

Question

Estoy tratando de escribir un filtro de paso bajo simple usando scipy, pero necesito ayuda para definir los parámetros.

Tengo 3.5 millones de registros en los datos de series de tiempo que necesitan ser filtrados, y los datos son muestreados a 1000 hz.

Estoy usando signal.firwin y signal.lfilter desde la biblioteca scipy.

Los parámetros que elijo en el siguiente código no filtran mis datos. En cambio, el siguiente código simplemente produce algo que se ve gráficamente como los mismos datos exactos excepto por una distorsión de fase temporal que desplaza el gráfico hacia la derecha en un poco menos de 1000 puntos de datos (1 segundo).

En otro programa de software, ejecutando un filtro de baja abeto pase a través de comandos de interfaz gráfica de usuario produce una salida que tiene medios similares para cada segmento de 10 segundos (10.000 puntos de datos), pero que tiene drásticamente menores desviaciones estándar de modo que esencialmente pierde la ruido en este archivo de datos en particular y reemplácelo con algo que retiene el valor medio mientras que muestra tendencias a más largo plazo que no están contaminadas por ruido de frecuencia más alta. El cuadro de diálogo de parámetros del otro software contiene una casilla de verificación que le permite seleccionar el número de coeficientes para que "se optimice en función del tamaño de muestra y la frecuencia de muestreo". (Los míos son de 3,5 millones de muestras recogidas a 1000 Hz, pero me gustaría una función que utiliza estas entradas como variables.)

* Puede alguien mostrará cómo ajustar el código de abajo de modo que elimina todas las frecuencias por encima de 0,05 hz? * Me gustaría ver ondas suaves en el gráfico en lugar de solo la distorsión temporal del mismo gráfico idéntico que obtengo ahora del código siguiente.

class FilterTheZ0(): 
    def __init__(self,ZSmoothedPylab): 
     #------------------------------------------------------ 
     # Set the order and cutoff of the filter 
     #------------------------------------------------------ 
     self.n = 1000 
     self.ZSmoothedPylab=ZSmoothedPylab 
     self.l = len(ZSmoothedPylab) 
     self.x = arange(0,self.l) 
     self.cutoffFreq = 0.05 

     #------------------------------------------------------ 
     # Run the filter 
     #------------------------------------------------------ 
     self.RunLowPassFIR_Filter(self.ZSmoothedPylab, self.n, self.l 
             , self.x, self.cutoffFreq) 

    def RunLowPassFIR_Filter(self,data, order, l, x, cutoffFreq): 
     #------------------------------------------------------ 
     # Set a to be the denominator coefficient vector 
     #------------------------------------------------------ 
     a = 1 
     #---------------------------------------------------- 
     # Create the low pass FIR filter 
     #---------------------------------------------------- 
     b = signal.firwin(self.n, cutoff = self.cutoffFreq, window = "hamming") 

     #--------------------------------------------------- 
     # Run the same data set through each of the various 
     # filters that were created above. 
     #--------------------------------------------------- 
     response = signal.lfilter(b,a,data) 
     responsePylab=p.array(response) 

     #-------------------------------------------------- 
     # Plot the input and the various outputs that are 
     # produced by running each of the various filters 
     # on the same inputs. 
     #-------------------------------------------------- 

     plot(x[10000:20000],data[10000:20000]) 
     plot(x[10000:20000],responsePylab[10000:20000]) 
     show() 

     return 

Answer 1

Las unidades de frecuencia de corte son probablemente [0,1) donde 1,0 es equivalente a FS (la frecuencia de muestreo). Entonces, si realmente quiere decir 0.05Hz y FS = 1000Hz, querría pasar cutoffFreq/1000. Es posible que necesite un filtro más largo para obtener un límite tan bajo.

(por cierto está de paso algunos argumentos, pero a continuación, utilizando los atributos del objeto en su lugar, pero no ve que la introducción de cualquier error obvias aún ...)

Answer 2

de corte se normaliza a la frecuencia de Nyquist, que es la mitad la tasa de muestreo. Entonces, con FS = 1000 y FC = 0.05, quiere un corte = 0.05/500 = 1e-4.

from scipy import signal 

FS = 1000.0           # sampling rate 
FC = 0.05/(0.5*FS)         # cutoff frequency at 0.05 Hz 
N = 1001            # number of filter taps 
a = 1            # filter denominator 
b = signal.firwin(N, cutoff=FC, window='hamming') # filter numerator 

M = FS*60           # number of samples (60 seconds) 
n = arange(M)          # time index 
x1 = cos(2*pi*n*0.025/FS)       # signal at 0.025 Hz 
x = x1 + 2*rand(M)         # signal + noise 
y = signal.lfilter(b, a, x)       # filtered output 

plot(n/FS, x); plot(n/FS, y, 'r')     # output in red 
grid() 

La salida del filtro se retrasa medio un segundo (el filtro está centrado en el grifo 500). Tenga en cuenta que el desplazamiento de CC añadido por el ruido se conserva mediante el filtro de paso bajo. Además, 0.025 Hz está dentro del rango de pase, por lo que el balanceo de salida de pico a pico es aproximadamente 2.