Solución:
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] (expresado 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
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 una falla. 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 más altas.
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)
para
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 las esquinas 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 las esquinas. Por lo tanto, sus bandas de parada comenzarían en 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 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.