lowpass - signal firwin python

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

Question

ACTUALIZAR:

Para mi sorpresa, mientras buscaba el mismo tema casi dos años después, ¡encontré una Receta Scipy basada en esta pregunta! Entonces, para cualquier persona interesada, ve directamente a:

http://wiki.scipy.org/Cookbook/ButterworthBandpass

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

Los parámetros que tengo que incluir son sample_rate, cutoff frequency IN HERTZ y posiblemente orden (otros parámetros, como atenuación, frecuencia natural, etc. son más oscuros para mí, por lo que cualquier valor "predeterminado").

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 poco claros, pero me gustaría implementar el formulario presentado en el comentario marcado como "para paso de banda". Los signos de interrogación en los comentarios muestran dónde simplemente copié-pegué algún ejemplo sin entender 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 directo en señales EMG.

¡Gracias por cualquier ayuda!

Answer 1

Para un filtro pasa banda, ws es una tupla que contiene las frecuencias de esquina inferior y superior. 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 la banda de paro. 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 stop.

Digamos, por ejemplo, que quería diseñar un filtro para una frecuencia de muestreo de 8000 muestras / seg 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 quería que las bandas de parada estuvieran a 30 dB +/- 100 Hz de las frecuencias de esquina. Por lo tanto, sus bandas de stop comenzarían a 200 y 3200 Hz, lo que daría como resultado 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 detención y de la pendiente de la curva de respuesta, que está determinada por la diferencia entre la frecuencia de esquina y la frecuencia de banda de detención.

Answer 2

Puede omitir el uso de buttord, y en su lugar solo elegir un orden 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, proporcione a butter () el orden de filtro, las frecuencias de corte Wn=[low, high] (expresadas como la fracción de la frecuencia Nyquist, que es la mitad de la frecuencia de muestreo) y el tipo de banda btype="band" .

Aquí hay una secuencia de comandos que define un par de funciones de conveniencia para trabajar con un filtro pasa banda Butterworth. Cuando se ejecuta como un script, hace dos tramas. Uno muestra la respuesta de frecuencia en varias órdenes de filtro para la misma frecuencia de muestreo y frecuencia de corte. La otra gráfica muestra el efecto del filtro (con orden = 6) en una serie temporal 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()

Estas son las tramas que genera este script: