Filtro de paso de banda digital Butterworth en Python

En este artículo, vamos a discutir cómo diseñar un filtro Butterworth de paso bajo 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 paso 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 y Lowpass:

La principal diferencia se puede detectar observando la respuesta de magnitud del filtro de paso de banda. La banda de paso del filtro es de un rango específico, lo que significa que el filtro de paso de banda puede pasar la única señal dentro de este rango. El filtro rechaza cualquier señal que no se encuentre dentro del rango especificado.

Las especificaciones son las siguientes:  

• Tasa de muestreo de 40 kHz
• Las frecuencias de borde de banda de paso son 1400 Hz y 2100 Hz
• Las frecuencias de borde de la banda de parada son 1050 Hz y 2450 Hz
• Ondulación de banda de paso de 0,4 dB
• Atenuación mínima de la banda de parada de 50 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 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

Ahora a continuación están los pasos:

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 = 7000
  
# pass band frequency
f_pass = [1400, 2100]
  
# stop band frequency
f_stop = [1050, 2450]
  
# pass band ripple
fs = 0.5
  
# Sampling Time
Td = 1
  
# pass band ripple
g_pass = 0.4
  
# stop band attenuation
g_stop = 50

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)

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()

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()

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()

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

# User-defined function to convert the 
# values of 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
  
# Importing Libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math
   
# Specifications of Filter 
    
# sampling frequency 
f_sample =7000
    
# pass band frequency 
f_pass =[1400,2100] 
    
# stop band frequency 
f_stop =[1050,2450]  
    
# pass band ripple 
fs = 0.5
    
# Sampling Time 
Td = 1  
    
 # pass band ripple 
g_pass = 0.4 
    
# stop band attenuation 
g_stop = 50  
  
# 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()