Filtro Butterworth de rechazo de banda digital en Python

En este artículo, vamos a discutir cómo diseñar un filtro Butterworth de rechazo de banda digital usando Python. El filtro Butterworth es un tipo de filtro de procesamiento de señal diseñado para tener una respuesta de frecuencia lo más plana posible en la banda de paso. Tomemos las siguientes especificaciones para diseñar el filtro y observar la respuesta de magnitud, fase e impulso del filtro digital Butterworth.

¿Qué es el filtro de rechazo de banda digital?

Un filtro de paso de banda es un filtro que pasa frecuencias dentro de un rango y rechaza frecuencias fuera de ese rango.

En qué se diferencia de Highpass, Lowpass y Bandpass:

La principal diferencia se puede detectar observando la respuesta de magnitud del filtro de paso de banda. En el filtro Band-Reject, todas las señales entre el rango de frecuencia especificado son rechazadas por el filtro. 

Las especificaciones son las siguientes:  

• La frecuencia de muestreo de 12 kHz.
• Las frecuencias de borde de la banda de paso son 2100 Hz y 4500 Hz.
• Las frecuencias de borde de la banda de parada son 2700 Hz y 3900 Hz.
• Ondulación de banda de paso de 0,6 dB.
• Atenuación mínima de la banda de parada de 45 dB.

Graficaremos la respuesta de magnitud, fase e impulso del filtro.

Enfoque paso a paso :

Antes de comenzar, primero, crearemos una función definida por el usuario para convertir las frecuencias de borde, la estamos definiendo como el método convert().

# explicit function to convert
# edge frequencies
  
  
def convertX(f_sample, f):
    w = []
  
    for i in range(len(f)):
        b = 2*((f[i]/2) / (f_sample/2))
        w.append(b)
  
    omega_mine = []
  
    for i in range(len(w)):
        c = (2/Td)*np.tan(w[i]/2)
        omega_mine.append(c)
  
    return omega_mine

Paso 1: Importación de todas las bibliotecas necesarias.

# import required modules
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math

Paso 2: Definir variables con las especificaciones dadas del filtro.

# Specifications of Filter
  
# sampling frequency
f_sample = 12000
  
# pass band frequency
f_pass = [2100, 4500]
  
# stop band frequency
f_stop = [2700, 3900]
  
# pass band ripple
fs = 0.5
  
# Sampling Time
Td = 1
  
# pass band ripple
g_pass = 0.6
  
# stop band attenuation
g_stop = 45

Paso 3: construir el filtro usando la función signal.buttord().

# Conversion to prewrapped analog
# frequency
omega_p = convertX(f_sample, f_pass)
omega_s = convertX(f_sample, f_stop)
  
# Design of Filter using signal.buttord
# function
N, Wn = signal.buttord(omega_p, omega_s,
                       g_pass, g_stop,
                       analog=True)
  
  
# Printing the values of order & cut-off frequency
# N is the order
print("Order of the Filter=", N)
  
# Wn is the cut-off freq of the filter
print("Cut-off frequency= {:} rad/s ".format(Wn))
  
  
# Conversion in Z-domain
  
# b is the numerator of the filter & a is
# the denominator
b, a = signal.butter(N, Wn, 'bandpass', True)
z, p = signal.bilinear(b, a, fs)
  
# w is the freq in z-domain & h is the
# magnitude in z-domain
w, h = signal.freqz(z, p, 512)

Producción:

Paso 4: Trazado de la Respuesta de Magnitud.

# Magnitude Response
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xscale('log')
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green')
plt.show()

Producción:

Paso 5: Representación gráfica de la respuesta al impulso.

# Impulse Response
imp = signal.unit_impulse(40)
c, d = signal.butter(N, 0.5)
response = signal.lfilter(c, d, imp)
  
plt.stem(np.arange(0, 40), imp, markerfmt='D', use_line_collection=True)
plt.stem(np.arange(0, 40), response, use_line_collection=True)
plt.margins(0, 0.1)
  
plt.xlabel('Time [samples]')
plt.ylabel('Amplitude')
plt.grid(True)
plt.show()

Producción:

Paso 6: Trazado de la respuesta de fase.

# Frequency Response
fig, ax1 = plt.subplots()
ax1.set_title('Digital filter frequency response')
ax1.set_ylabel('Angle(radians)', color='g')
ax1.set_xlabel('Frequency [Hz]')
  
angles = np.unwrap(np.angle(h))
  
ax1.plot(w/2*np.pi, angles, 'g')
ax1.grid()
ax1.axis('tight')
plt.show()

Producción:

A continuación se muestra el programa completo basado en el enfoque anterior:

# import required modules
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math
  
# explicit function to convert
# edge frequencies
def convertX(f_sample, f):
    w = []
  
    for i in range(len(f)):
        b = 2*((f[i]/2)/(f_sample/2))
        w.append(b)
  
    omega_mine = []
  
    for i in range(len(w)):
        c = (2/Td)*np.tan(w[i]/2)
        omega_mine.append(c)
  
    return omega_mine
  
# Specifications of Filter
# sampling frequency
f_sample = 12000
  
# pass band frequency
f_pass = [2100, 4500]
  
# stop band frequency
f_stop = [2700, 3900]
  
# pass band ripple
fs = 0.5
  
# Sampling Time
Td = 1
  
# pass band ripple
g_pass = 0.6
  
# stop band attenuation
g_stop = 45
  
# Conversion to prewrapped analog
# frequency
omega_p = convertX(f_sample, f_pass)
omega_s = convertX(f_sample, f_stop)
  
# Design of Filter using signal.buttord
# function
N, Wn = signal.buttord(omega_p, omega_s, g_pass,
                       g_stop, analog=True)
  
# Printing the values of order & cut-off frequency
# N is the order
print("Order of the Filter=", N)
  
# Wn is the cut-off freq of the filter
print("Cut-off frequency= {:} rad/s ".format(Wn))
  
  
# Conversion in Z-domain
  
# b is the numerator of the filter & a is
# the denominator
b, a = signal.butter(N, Wn, 'bandpass', True)
z, p = signal.bilinear(b, a, fs)
  
# w is the freq in z-domain & h is the
# magnitude in z-domain
w, h = signal.freqz(z, p, 512)
  
# Magnitude Response
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xscale('log')
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green')
plt.show()
  
# Impulse Response
imp = signal.unit_impulse(40)
c, d = signal.butter(N, 0.5)
response = signal.lfilter(c, d, imp)
  
plt.stem(np.arange(0, 40), imp,
         markerfmt='D',
         use_line_collection=True)
  
plt.stem(np.arange(0, 40), response,
         use_line_collection=True)
  
plt.margins(0, 0.1)
  
plt.xlabel('Time [samples]')
plt.ylabel('Amplitude')
plt.grid(True)
plt.show()
  
# Frequency Response
fig, ax1 = plt.subplots()
ax1.set_title('Digital filter frequency response')
ax1.set_ylabel('Angle(radians)', color='g')
ax1.set_xlabel('Frequency [Hz]')
  
angles = np.unwrap(np.angle(h))
  
ax1.plot(w/2*np.pi, angles, 'g')
ax1.grid()
ax1.axis('tight')
plt.show()

Producción: