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

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

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.

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:

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)

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.