Bandpass frecuencia de filtro de Butterworth en scipy

Estoy diseñando un filtro de paso de banda en Scipy siguiendo el libro de cocina . Sin embargo, si disminuyo demasiado las frecuencias de filtrado, termino con basura en los filtros de alto orden. ¿Qué estoy haciendo mal?

from scipy.signal import butter, lfilter def butter_bandpass(lowcut, highcut, fs, order=5): nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = butter(order, [low, high], btype='band') return b, a if __name__ == "__main__": import numpy as np import matplotlib.pyplot as plt from scipy.signal import freqz # Sample rate and desired cutoff frequencies (in Hz). fs = 25 # Plot the frequency response for a few different orders. plt.figure(1) plt.clf() for order in [1, 3, 5, 6, 9]: b, a = butter_bandpass(0.5, 4, fs, order=order) w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000)) plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order) plt.xlabel('Frequency (Hz)') plt.ylabel('Gain') plt.grid(True) plt.legend(loc='best') plt.figure(2) plt.clf() for order in [1, 3, 5, 6, 9]: b, a = butter_bandpass(0.05, 0.4, fs, order=order) w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000)) plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order) plt.xlabel('Frequency (Hz)') plt.ylabel('Gain') plt.grid(True) plt.legend(loc='best') plt.show() 

fs = 25, bajo = 0.5, alto = 4fs = 25, bajo = 0.05, alto = 0.4

Actualización: el problema se discutió y aparentemente se resolvió en Scipy 0.14. Sin embargo, la ttwig todavía se ve muy mal después de la actualización de Scipy. Que pasa

Después de Scipy 0.14, aún peor.

  1. No use b, a = butter para filtros de alto orden, ya sea en Matlab o SciPy u Octave. El formato de la función de transferencia tiene problemas de estabilidad numérica , porque algunos de los coeficientes son muy grandes mientras que otros son muy pequeños. Es por esto que cambiamos las funciones de diseño del filtro para usar el formato zpk internamente . Para ver los beneficios de esto, necesita usar z, p, k = butter(output='zpk') y luego trabajar con polos y ceros en lugar de numerador y denominador.
  2. No haga filtros digitales de alto orden en una sola etapa. Esta es una mala idea, sin importar en qué software o hardware los esté implementando. Normalmente, lo mejor es dividirlos en secciones de segundo orden . En Matlab, puedes usar zp2sos para generar estos automáticamente. Esto aún no se ha implementado en SciPy , aunque tengo un código GPL para hacerlo aquí: https://gist.github.com/endolith/4525003 (la licencia GPL es incompatible con scipy, por lo que debe volver a escribirse desde cero).

Aparentemente el problema es un error conocido:

Github

Este es un problema común en los filtros digitales. Los filtros de alto orden con frecuencias de corte muy por debajo de la frecuencia de nyquist tienden a tener coeficientes inestables debido a la precisión limitada de los números de punto flotante. La última vez que lo verifiqué (es cierto que hace un par de años), Matlab hizo un trabajo mucho mejor de conservar la precisión que el scipy, aunque todavía dará problemas con filtros suficientemente extremos.

Hay un par de opciones si no puedes usar matlab. La primera es dividir su filtro en secciones de segundo orden en cascada. Básicamente, calcula los polos y ceros deseados, los divide en pares conjugados complejos y calcula la función de transferencia para cada par.

La segunda opción es remuestrear a una frecuencia de muestreo más similar a la frecuencia del filtro. Por ejemplo, en su segundo ejemplo, su frecuencia de muestreo es de 25 y su frecuencia de corte más alta es .4. Puede usar un filtro anti-aliasing de paso bajo y luego diezmar por un factor de 10 a una frecuencia de muestreo de 2.5. Con la tasa de muestreo más baja, sus coeficientes de filtro de paso de banda serán menos sensibles a los errores de redondeo. Si hace esto, debe asegurarse de que el filtro de suavizado no tenga el mismo problema.

Lo que sucede es que los órdenes de los filtros de paso de banda (BP) creados en el script son, de hecho, el doble de los que se muestran en el gráfico . Recuerde que el orden del filtro es el orden del polinomio en el denominador de la función de transferencia. Los filtros de paso de banda Canonic son siempre de orden .

Esos números que se muestran son los órdenes del prototipo de paso bajo (LP) (normalmente normalizado a una frecuencia de corte de 1 rad / s) que se utiliza para aplicar la transformación de LP a BP que duplica el orden del filtro. Entonces, por ejemplo, si comenzamos con un LP de orden 1, terminamos con un paso de banda de segundo orden:

1 / (S + 1) => LP-2-BP transf. => ks / (s ^ 2 + a.s + b)

Donde k , ayb son constantes. El numerador de un filtro de paso de banda estándar es ks ^ (N / 2), por lo que el orden N del filtro debe ser par.

Este problema de orden para el paso de banda (que también sucede con las muescas o los filtros de rechazo de banda) no se menciona en la documentación de SciPy. De hecho, si imprime la longitud del denominador a (utilizando print(len(a)) ) antes de plt.show() , verá que tiene 19 coeficientes, lo que corresponde a un polinomio de 18º orden.