Comment mettre en œuvre band-pass Butterworth filter avec Scipy.signal.beurre
mise à jour:
j'ai trouvé une recette Scipy basée sur cette question! Donc, pour toute personne intéressée, allez directement à:
http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
j'ai du mal à réaliser ce qui semblait initialement une simple tâche de mise en œuvre D'un filtre Butterworth band-pass pour 1-d numpy array (time-series).
le les paramètres que je dois inclure sont le taux d'échantillonnage, les fréquences de coupure en HERTZ et éventuellement l'ordre (autres paramètres, comme l'atténuation, la fréquence naturelle, etc. sont plus obscurs pour moi, donc n'importe quelle valeur "par défaut" ferait).
ce que j'ai maintenant est ceci, qui semble fonctionner comme un filtre passe-haut mais je ne suis pas sûr si je le fais 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
les documents et les exemples portent à confusion et obscur, mais j'aimerais mettre en œuvre le formulaire présenté dans la mention élogieuse marquée comme "pour bandpass". Les points d'interrogation dans les commentaires montrent où je viens de copier-collé un exemple sans comprendre ce qui se passe.
Je ne suis pas un ingénieur électricien ou un scientifique, juste un concepteur d'équipement médical qui doit effectuer un filtrage bandpass assez simple sur les signaux EMG.
Merci pour votre aide!
3 réponses
vous pouvez sauter l'utilisation de buttord, et à la place, il suffit de choisir une commande pour le filtre et voir si elle répond à votre critère de filtrage. Pour générer les coefficients de filtre pour un filtre bandpass, donner à butter() l'ordre du filtre, les fréquences de coupure Wn=[low, high] (exprimées en fraction de la fréquence Nyquist, qui est la moitié de la fréquence d'échantillonnage) et le type de bande btype="band" .
voici un script qui définit quelques fonctions pratiques pour travailler avec un filtre Butterworth bandpass. Lorsque vous exécutez un script, il fait deux parcelles. L'une 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 graphique montre l'effet du filtre (ordre=6) sur une série temporelle d'échantillon.
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()
Voici les tracés qui sont générés par ce script:
la méthode de conception du filtre dans réponse acceptée est correcte, mais elle présente un défaut. Les filtres scipy bandpass conçus avec b, a sont unstable et peuvent avoir pour résultat filtres erronés à ordres de filtres supérieurs .
utilisez plutôt la sortie sos (sections de deuxième 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
Aussi, vous pouvez tracer réponse de fréquence en changeant
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)
pour un filtre bandpass, ws est un tuple contenant les fréquences des coins inférieur et supérieur. Ceux-ci représentent la fréquence numérique où la réponse du filtre est inférieure de 3 dB à la bande passante.
wp est un tuple contenant les fréquences numériques de bande d'arrêt. Ils représentent l'endroit où l'atténuation maximale commence.
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.
dire, par exemple, vous vouliez concevoir un filtre pour une fréquence d'échantillonnage de 8000 échantillons/sec ayant des fréquences d'angle de 300 et 3100 Hz. La fréquence de Nyquist est la fréquence d'échantillonnage divisée par deux, ou dans cet exemple, 4000 Hz. La fréquence numérique équivalente est de 1,0. Les deux fréquences d'angle sont alors 300/4000 et 3100/4000.
permet maintenant de dire que vous vouliez que les stopbands soient en bas de 30 dB +/- 100 Hz des fréquences d'angle. Ainsi, vos stopbands commenceraient à 200 et 3200 Hz résultant en des fréquences numériques de 200/4000 et 3200/4000.
pour créer votre filtre, vous appelleriez buttord comme
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 de butée 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 bande de butée.