Diseñe el filtro IIR Lowpass Butterworth usando el método de transformación bilineal en Scipy-Python

IIR significa Infinite Impulse Response. Es una de las características sorprendentes de muchos sistemas invariantes en el tiempo lineal que se distinguen por tener una respuesta de impulso h(t)/h(n) que no se vuelve cero después de algún punto, sino que continúa infinitamente . .

¿Qué es IIR Lowpass Butterworth?

Básicamente se comporta como un filtro Butterworth de paso bajo digital ordinario con una respuesta de impulso infinita. 

Las especificaciones son las siguientes:  

• Tasa de muestreo de 8 kHz
• Orden de Filtro 2
• Frecuencia de corte 3400Hz

Trazaremos la respuesta de magnitud y fase del filtro.

Enfoque paso a paso:

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

# import required library
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

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

# Given specification
N = 2  # Order of the filter
Fs = 8000  # Sampling frequency in Hz
fc = 3400  # Cut-off frequency in Hz
 
# Compute Design Sampling parameter
Td = 1/Fs

Paso 3: Cálculo de la frecuencia de corte

# Compute cut-off frequency in radian/sec
wd = 2*np.pi*fc
print(wd)  # Cut-off frequency in radian/sec

Producción:

Paso 4: Preenvolver la frecuencia analógica

# Prewarp the analog frequency
 
wc = (2/Td)*np.tan(wd*Td/2)
print('Order of the filter=', N)  # Order
 
# Prewarped analog cut-off frequency
print('Cut-off frequency (in rad/s)=', wc)

Producción:

Paso 5: diseñar el filtro usando la función signal.butter() y luego realizar una transformación bilineal usando la función signal.bilinear()

# Design analog Butterworth filter using signal.butter function
 
b, a = signal.butter(N, wc, 'low', analog='True')
# Perform bilinear Transformation
 
z, p = signal.bilinear(b, a, fs=Fs)
 
# Print numerator and denomerator coefficients of the filter
print('Numerator Coefficients:', z)
print('Denominator Coefficients:', p)

Producción:

Paso 6: Calcular la respuesta de frecuencia del filtro usando la función signal.freqz() y trazar la respuesta de magnitud y fase

# Compute frequency response of the filter using signal.freqz function
wz, hz = signal.freqz(z, p, 512)
 
# Plot filter magnitude and phase responses using subplot.
# Convert digital frequency wz into analog frequency in Hz
fig = plt.figure(figsize=(12, 10))
 
# Calculate Magnitude from hz in dB
Mag = 20*np.log10(abs(hz))
 
# Calculate frequency in Hz from wz
Freq = wz*Fs/(2*np.pi)
 
# Plot Magnitude response
sub1 = plt.subplot(2, 1, 1)
sub1.plot(Freq, Mag, 'r', linewidth=2)
sub1.axis([1, Fs/2, -60, 5])
sub1.set_title('Magnitude Response', fontsize=15)
sub1.set_xlabel('Frequency [Hz]', fontsize=15)
sub1.set_ylabel('Magnitude [dB]', fontsize=15)
sub1.grid()
 
 
# Plot phase angle
sub2 = plt.subplot(2, 1, 2)
 
# Calculate phase angle in degree from hz
Phase = np.unwrap(np.angle(hz))*180/np.pi
sub2.plot(Freq, Phase, 'g', linewidth=2)
sub2.set_ylabel('Phase (degree)', fontsize=15)
sub2.set_xlabel(r'Frequency (Hz)', fontsize=15)
sub2.set_title(r'Phase response', fontsize=15)
sub2.grid()
 
plt.subplots_adjust(hspace=0.5)
fig.tight_layout()
plt.show()

Producción:

A continuación se muestra la implementación:

# import required library
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
 
# Given specification
N = 2  # Order of the filter
Fs = 8000  # Sampling frequency in Hz
fc = 3400  # Cut-off frequency in Hz
 
# Compute Design Sampling parameter
Td = 1/Fs
 
# Compute cut-off frequency in radian/sec
wd = 2*np.pi*fc
print(wd)  # Cut-off frequency in radian/sec
 
# Prewarp the analog frequency
wc = (2/Td)*np.tan(wd*Td/2)
print('Order of the filter=', N)  # Order
 
# Prewarped analog cut-off frequency
print('Cut-off frequency (in rad/s)=', wc)
 
# Design analog Butterworth filter using signal.butter function
b, a = signal.butter(N, wc, 'low', analog='True')
 
# Perform bilinear Transformation
z, p = signal.bilinear(b, a, fs=Fs)
 
# Print numerator and denomerator coefficients of the filter
print('Numerator Coefficients:', z)
print('Denominator Coefficients:', p)
 
# Compute frequency response of the filter using signal.freqz function
wz, hz = signal.freqz(z, p, 512)
 
# Plot filter magnitude and phase responses using subplot.
#Convert digital frequency wz into analog frequency in Hz
fig = plt.figure(figsize=(10, 8))
 
# Calculate Magnitude from hz in dB
Mag = 20*np.log10(abs(hz))
 
# Calculate frequency in Hz from wz
Freq = wz*Fs/(2*np.pi)
 
# Plot Magnitude response
sub1 = plt.subplot(2, 1, 1)
sub1.plot(Freq, Mag, 'r', linewidth=2)
sub1.axis([1, Fs/2, -60, 5])
sub1.set_title('Magnitude Response', fontsize=15)
sub1.set_xlabel('Frequency [Hz]', fontsize=15)
sub1.set_ylabel('Magnitude [dB]', fontsize=15)
sub1.grid()
 
# Plot phase angle
sub2 = plt.subplot(2, 1, 2)
 
# Calculate phase angle in degree from hz
Phase = np.unwrap(np.angle(hz))*180/np.pi
sub2.plot(Freq, Phase, 'g', linewidth=2)
sub2.set_ylabel('Phase (degree)', fontsize=15)
sub2.set_xlabel(r'Frequency (Hz)', fontsize=15)
sub2.set_title(r'Phase response', fontsize=15)
sub2.grid()
 
plt.subplots_adjust(hspace=0.5)
fig.tight_layout()
plt.show()

Producción: