Comment implémenter un filtre Butterworth passe-bande avec Scipy.signal.butter

Question

UPDATE :

J'ai trouvé une recette Scipy basée sur cette question ! Donc, pour toute personne intéressée, allez directement à : Contenu " Traitement du signal " Passe-bande de Butterworth

---

J'ai du mal à réaliser ce qui semblait initialement une tâche simple, à savoir la mise en œuvre d'un filtre passe-bande Butterworth pour un tableau numpy 1-D (séries temporelles).

Les paramètres que je dois inclure sont le sample_rate, les fréquences de coupure en HERTZ et éventuellement l'ordre (d'autres paramètres, comme l'atténuation, la fréquence naturelle, etc. sont plus obscurs pour moi, donc toute valeur "par défaut" ferait l'affaire).

Ce que j'ai maintenant est ceci, qui semble fonctionner comme un filtre passe-haut mais je ne suis pas sûr de le faire correctement :

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

La documentation et les exemples sont confus et obscurs, mais j'aimerais mettre en œuvre la forme présentée dans le commentaire marqué comme "pour le passe-bande". Les points d'interrogation dans les commentaires montrent les endroits où j'ai simplement copié-collé un exemple sans comprendre ce qui se passe.

Je ne suis ni ingénieur électricien ni scientifique, juste un concepteur d'équipement médical qui doit effectuer un filtrage passe-bande assez simple sur des signaux EMG.

Answer 1

Vous pourriez éviter l'utilisation de buttord, et à la place choisir simplement un ordre pour le filtre et voir s'il répond à votre critère de filtrage. Pour générer les coefficients de filtre pour un filtre passe-bande, donnez à butter() l'ordre du filtre, les fréquences de coupure Wn=[lowcut, highcut] la fréquence d'échantillonnage fs (exprimée dans les mêmes unités que les fréquences de coupure) et le type de bande btype="band" .

Voici un script qui définit quelques fonctions pratiques pour travailler avec un filtre passe-bande de Butterworth. Lorsqu'il est exécuté comme un script, il réalise deux tracés. L'un montre la réponse en fréquence à plusieurs ordres de filtre pour le même taux d'échantillonnage et les mêmes fréquences de coupure. L'autre tracé démontre l'effet du filtre (avec l'ordre=6) sur une série temporelle d'échantillons.

from scipy.signal import butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order=5):
    return butter(order, [lowcut, highcut], fs=fs, btype='band')

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, fs=fs, worN=2000)
        plt.plot(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.arange(0, nsamples) / fs
    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()

Voici les graphiques qui sont générés par ce script :

Answer 2

La méthode de conception des filtres dans réponse acceptée est correcte, mais elle a un défaut. Les filtres passe-bande SciPy conçus avec b, a sont instable et peut entraîner filtres erronés en ordres de filtre supérieurs .

Utilisez plutôt la sortie sos (sections de second ordre) de la conception du filtre.

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

Vous pouvez également tracer la réponse en fréquence en modifiant

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

à

sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)

Answer 3

Pour un filtre passe-bande, ws est un tuple contenant les fréquences de coin inférieure et supérieure. Celles-ci représentent la fréquence numérique où la réponse du filtre est inférieure de 3 dB à celle de la bande passante.

wp est un tuple contenant les fréquences numériques de la bande d'arrêt. Elles représentent l'endroit où commence l'atténuation maximale.

gpass est l'atténuation maximale dans la bande passante en dB tandis que gstop est l'atténuation dans les bandes d'arrêt.

Supposons, par exemple, que vous vouliez concevoir un filtre pour un taux d'échantillonnage de 8000 échantillons/seconde ayant des fréquences de coin de 300 et 3100 Hz. La fréquence de Nyquist est la fréquence d'échantillonnage divisée par deux, soit dans cet exemple, 4000 Hz. La fréquence numérique équivalente est de 1,0. Les deux fréquences de coin sont donc 300/4000 et 3100/4000.

Supposons maintenant que vous vouliez que les bandes d'arrêt soient réduites de 30 dB +/- 100 Hz par rapport aux fréquences de coin. Ainsi, vos bandes d'arrêt commenceraient à 200 et 3200 Hz, ce qui donnerait des fréquences numériques de 200/4000 et 3200/4000.

Pour créer votre filtre, vous devez appeler buttord comme suit

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 longueur du filtre résultant dépendra de la profondeur des bandes d'arrêt et de la pente de la courbe de réponse qui est déterminée par la différence entre la fréquence de coin et la fréquence de la bande d'arrêt.