Cómo implementar el filtro Butterworth de paso de banda con Scipy.signal.butter

83

ACTUALIZAR:

¡Encontré una Receta Scipy basada en esta pregunta! Entonces, para cualquier persona interesada, vaya directamente a: Contenido »Procesamiento de señales» Butterworth Bandpass


Estoy teniendo dificultades para lograr lo que inicialmente parecía una tarea simple de implementar un filtro de paso de banda de Butterworth para una matriz numpy 1-D (serie temporal).

Los parámetros que tengo que incluir son el sample_rate, las frecuencias de corte EN HERTZ y posiblemente el orden (otros parámetros, como atenuación, frecuencia natural, etc. son más oscuros para mí, por lo que cualquier valor "predeterminado" sería suficiente).

Lo que tengo ahora es esto, que parece funcionar como un filtro de paso alto, pero no estoy seguro de si lo estoy haciendo bien:

def butter_highpass(interval, sampling_rate, cutoff, order=5):
    nyq = sampling_rate * 0.5

    stopfreq = float(cutoff)
    cornerfreq = 0.4 * stopfreq  # (?)

    ws = cornerfreq/nyq
    wp = stopfreq/nyq

    # for bandpass:
    # wp = [0.2, 0.5], ws = [0.1, 0.6]

    N, wn = scipy.signal.buttord(wp, ws, 3, 16)   # (?)

    # for hardcoded order:
    # N = order

    b, a = scipy.signal.butter(N, wn, btype='high')   # should 'high' be here for bandpass?
    sf = scipy.signal.lfilter(b, a, interval)
    return sf

ingrese la descripción de la imagen aquí

Los documentos y ejemplos son confusos y oscuros, pero me gustaría implementar el formulario presentado en el elogio marcado como "para paso de banda". Los signos de interrogación en los comentarios muestran dónde simplemente copié y pegué un ejemplo sin comprender lo que está sucediendo.

No soy ingeniero eléctrico ni científico, solo un diseñador de equipos médicos que necesita realizar un filtrado de paso de banda bastante sencillo en las señales EMG.

heltonbiker
fuente
Probé algo en dsp.stackexchange, pero se enfocan demasiado (más de lo que puedo manejar) en cuestiones conceptuales de ingeniería y no tanto en el uso de funciones scipy.
heltonbiker

Respuestas:

117

Puede omitir el uso de buttord y, en su lugar, simplemente elegir un pedido para el filtro y ver si cumple con su criterio de filtrado. Para generar los coeficientes de filtro para un filtro de paso de banda, asigne a butter () el orden del filtro, las frecuencias de corte Wn=[low, high](expresadas como la fracción de la frecuencia de Nyquist, que es la mitad de la frecuencia de muestreo) y el tipo de banda btype="band".

Aquí hay un script que define un par de funciones de conveniencia para trabajar con un filtro de paso de banda de Butterworth. Cuando se ejecuta como un script, crea dos gráficos. Uno muestra la respuesta de frecuencia en varios órdenes de filtro para la misma frecuencia de muestreo y frecuencias de corte. El otro gráfico demuestra el efecto del filtro (con orden = 6) en una serie de tiempo de muestra.

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


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


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 = 5000.0
    lowcut = 500.0
    highcut = 1250.0

    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [3, 6, 9]:
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        w, h = freqz(b, a, worN=2000)
        plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)

    plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    # Filter a noisy signal.
    T = 0.05
    nsamples = T * fs
    t = np.linspace(0, T, nsamples, endpoint=False)
    a = 0.02
    f0 = 600.0
    x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    x += a * np.cos(2 * np.pi * f0 * t + .11)
    x += 0.03 * np.cos(2 * np.pi * 2000 * t)
    plt.figure(2)
    plt.clf()
    plt.plot(t, x, label='Noisy signal')

    y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
    plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
    plt.xlabel('time (seconds)')
    plt.hlines([-a, a], 0, T, linestyles='--')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc='upper left')

    plt.show()

Estos son los gráficos que genera este script:

Respuesta de frecuencia para varias órdenes de filtro

ingrese la descripción de la imagen aquí

Warren Weckesser
fuente
1
¿Sabes por qué la salida filtrada siempre comienza en el valor cero? ¿Es posible emparejarlo con el valor de entrada real x[0]? Probé cosas similares con el filtro de paso bajo Cheby1 y tuve el mismo problema.
LWZ
2
@LWZ: usa la función scipy.signal.lfilter_ziy el ziargumento para lfilter. Para obtener más información, consulte la cadena de documentación lfilter_zi. TL; DR? Simplemente cambie y = lfilter(b, a, data)a zi = lfilter_zi(b, a); y, zo = lfilter(b, a, data, zi=zi*data[0]). (Pero esto podría no marcar la diferencia con un filtro de paso de banda o de paso alto.)
Warren Weckesser
1
Noté que hay un cambio de fase de 180 grados en la salida de scipy.signal.lfiter()wrt a la señal original y la signal.filtfilt()salida, ¿por qué? ¿Debería usarlo filtfilt()en su lugar si el tiempo es importante para mí?
Jason
1
Ese es el retardo de fase del filtro a esa frecuencia. El retardo de fase de una sinusoide a través de un filtro Butterworth depende de forma no lineal de la frecuencia. Para retardo de fase cero, sí, se puede utilizar filtfilt(). Mi respuesta aquí incluye un ejemplo de uso filtfilt()para evitar un retraso inducido por el filtro.
Warren Weckesser
1
Hola Jason, recomiendo hacer preguntas sobre la teoría del procesamiento de señales en dsp.stackexchange.com . Si tiene una pregunta sobre algún código que escribió que no funciona como se esperaba, puede comenzar una nueva pregunta aquí en stackoverflow.
Warren Weckesser
37

El método de diseño de filtro en la respuesta aceptada es correcto, pero tiene un defecto. Los filtros de paso de banda SciPy diseñados con b, a son inestables y pueden resultar en filtros erróneos en órdenes de filtro superiores .

En su lugar, utilice la salida sos (secciones de segundo orden) del diseño del filtro.

from scipy.signal import butter, sosfilt, sosfreqz

def butter_bandpass(lowcut, highcut, fs, order=5):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        sos = butter(order, [low, high], analog=False, btype='band', output='sos')
        return sos

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        y = sosfilt(sos, data)
        return y

Además, puede trazar la respuesta de frecuencia cambiando

b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)

a

sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)
user13107
fuente
+1 porque ahora es la mejor manera de hacerlo en muchos casos. Al igual que en los comentarios sobre la respuesta aceptada, también es posible eliminar el retraso de fase mediante el filtrado hacia adelante y hacia atrás. Simplemente reemplácelo sosfiltcon sosfiltfilt.
Mike
@Mike y user13107 ¿El mismo error afecta también a los filtros Butterworth de paso alto y paso bajo? ¿Y la solución es la misma?
dewarrn1
3
@ dewarrn1 En realidad, no es correcto llamarlo "error"; el algoritmo está implementado correctamente, pero es inherentemente inestable, por lo que es una mala elección de algoritmo. Pero sí, afecta a cualquier filtro de orden superior, no solo a los filtros de paso alto o bajo y no solo a los de Butterworth, sino también a otros como Chebyshev, etc. De todos modos, en general, es mejor elegir siempre la sossalida, porque eso siempre evitará la inestabilidad. Y a menos que necesite procesamiento en tiempo real, siempre debe usar sosfiltfilt.
Mike
¡Lamento no haber notado esta respuesta hace mucho tiempo! @ user13107, sí, la representación de la función de transferencia (o 'ba') de un filtro lineal tiene serios problemas numéricos cuando el orden del filtro es grande. Incluso los filtros de orden relativamente bajo pueden tener problemas cuando el ancho de banda deseado es pequeño en comparación con la frecuencia de muestreo. Mi respuesta original se escribió antes de que se agregara la representación SOS a SciPy, y antes de fsque se agregara el argumento a muchas de las funciones en scipy.signal. La respuesta se debe actualizar desde hace mucho tiempo.
Warren Weckesser
alguna ayuda para esto? stackoverflow.com/q/60866193/5025009
seralouk
4

Para un filtro de paso de banda, ws es una tupla que contiene las frecuencias de esquina superior e inferior. Estos representan la frecuencia digital donde la respuesta del filtro es 3 dB menor que la banda de paso.

wp es una tupla que contiene las frecuencias digitales de banda de parada. Representan la ubicación donde comienza la atenuación máxima.

gpass es la atenuación máxima en la banda de paso en dB mientras que gstop es la atenuación en las bandas de parada.

Digamos, por ejemplo, que desea diseñar un filtro para una frecuencia de muestreo de 8000 muestras / segundo con frecuencias de esquina de 300 y 3100 Hz. La frecuencia de Nyquist es la frecuencia de muestreo dividida por dos, o en este ejemplo, 4000 Hz. La frecuencia digital equivalente es 1.0. Las dos frecuencias de esquina son entonces 300/4000 y 3100/4000.

Ahora digamos que desea que las bandas de parada estén por debajo de 30 dB +/- 100 Hz desde las frecuencias de esquina. Por lo tanto, sus bandas de parada comenzarían en 200 y 3200 Hz, lo que resultaría en las frecuencias digitales de 200/4000 y 3200/4000.

Para crear su filtro, llamaría a buttord como

fs = 8000.0
fso2 = fs/2
N,wn = scipy.signal.buttord(ws=[300/fso2,3100/fso2], wp=[200/fs02,3200/fs02],
   gpass=0.0, gstop=30.0)

La longitud del filtro resultante dependerá de la profundidad de las bandas de parada y de la inclinación de la curva de respuesta que está determinada por la diferencia entre la frecuencia de esquina y la frecuencia de la banda de parada.

sizzzzlerz
fuente
Traté de implementarlo, pero todavía falta algo. Una cosa es que gpass=0.0genera una división por error cero, así que lo cambié a 0.1 y el error se detuvo. Además de eso, los documentos butterdicen: Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample).Tengo dudas de si su respuesta hizo los cálculos correctamente, así que todavía estoy trabajando en eso y pronto daré algunos comentarios.
heltonbiker
(también, aunque my wsy wptienen dos elementos cada uno, el filtro solo realiza paso alto o bajo (a través del btypeargumento), pero no paso de banda)
heltonbiker
1
De acuerdo con la documentación en docs.scipy.org/doc/scipy/reference/generated/… , buttord diseña filtros de paso de banda, alto y bajo. En cuanto a gpass, supongo que buttord no permite una atenuación de 0 dB en la banda de paso. Ajústelo a un valor distinto de cero entonces.
sizzzzlerz