Разработайте БИХ-фильтр верхних частот Баттерворта с использованием метода билинейного преобразования в Scipy - Python

Опубликовано: 6 Января, 2022

IIR расшифровывается как Infinite Impulse Response, это одна из ярких особенностей многих систем, инвариантных по линейному времени, которые отличаются наличием импульсной характеристики h (t) / h (n), которая не становится нулевой после некоторой точки, а вместо этого продолжается бесконечно. .

Что такое IIR Highpass Butterworth?

По сути, он ведет себя так же, как обычный цифровой высокочастотный фильтр Баттерворта с бесконечной импульсной характеристикой.

Технические характеристики следующие:

Частота полосы пропускания: 2-4 кГц
Частота стоп-полосы: 0-500 Гц
Пульсация полосы пропускания: 3 дБ
Затухание в полосе задерживания: 20 дБ
Частота дискретизации: 8 кГц
Мы построим график амплитуды, фазы, импульса, переходной характеристики фильтра.

Пошаговый подход:

Шаг 1: Импорт всех необходимых библиотек.

Python3

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

Шаг 2: Определение пользовательских функций mfreqz () и impz () . mfreqz - это функция для амплитуды и фазового графика, а impz - это функция для импульсной и ступенчатой характеристики.

Python3

 def mfreqz(b, a, Fs):
   
    # Compute frequency response of the filter
    # using signal.freqz function
    wz, hz = signal.freqz(b, a)
 
    # Calculate Magnitude from hz in dB
    Mag = 20 * np.log10( abs (hz))
 
    # Calculate phase angle in degree from hz
    Phase = np.unwrap(np.arctan2(np.imag(hz), np.real(hz))) * ( 180 / np.pi)
 
    # Calculate frequency in Hz from wz
    Freq = wz * Fs / ( 2 * np.pi) # START CODE HERE ### (≈ 1 line of code)
 
    # Plot filter magnitude and phase responses using subplot.
    fig = plt.figure(figsize = ( 10 , 6 ))
 
    # Plot Magnitude response
    sub1 = plt.subplot( 2 , 1 , 1 )
    sub1.plot(Freq, Mag, 'r' , linewidth = 2 )
    sub1.axis([ 1 , Fs / 2 , - 100 , 5 ])
    sub1.set_title( 'Magnitute Response' , fontsize = 20 )
    sub1.set_xlabel( 'Frequency [Hz]' , fontsize = 20 )
    sub1.set_ylabel( 'Magnitude [dB]' , fontsize = 20 )
    sub1.grid()
 
    # Plot phase angle
    sub2 = plt.subplot( 2 , 1 , 2 )
    sub2.plot(Freq, Phase, 'g' , linewidth = 2 )
    sub2.set_ylabel( 'Phase (degree)' , fontsize = 20 )
    sub2.set_xlabel(r 'Frequency (Hz)' , fontsize = 20 )
    sub2.set_title(r 'Phase response' , fontsize = 20 )
    sub2.grid()
 
    plt.subplots_adjust(hspace = 0.5 )
    fig.tight_layout()
    plt.show()
 
# Define impz(b,a) to calculate impulse response
 # and step response of a system input: b= an array
 # containing numerator coefficients,a= an array containing
 #denominator coefficients
 def impz(b, a):
     
    # Define the impulse sequence of length 60
    impulse = np.repeat( 0. , 60 )
    impulse[ 0 ] = 1.
    x = np.arange( 0 , 60 )
 
    # Compute the impulse response
    response = signal.lfilter(b, a, impulse)
 
    # Plot filter impulse and step response:
    fig = plt.figure(figsize = ( 10 , 6 ))
    plt.subplot( 211 )
    plt.stem(x, response, 'm' , use_line_collection = True )
    plt.ylabel( 'Amplitude' , fontsize = 15 )
    plt.xlabel(r 'n (samples)' , fontsize = 15 )
    plt.title(r 'Impulse response' , fontsize = 15 )
 
    plt.subplot( 212 )
    step = np.cumsum(response) # Compute step response of the system
    plt.stem(x, step, 'g' , use_line_collection = True )
    plt.ylabel( 'Amplitude' , fontsize = 15 )
    plt.xlabel(r 'n (samples)' , fontsize = 15 )
    plt.title(r 'Step response' , fontsize = 15 )
    plt.subplots_adjust(hspace = 0.5 )
 
    fig.tight_layout()
    plt.show()

Шаг 3: Определите переменные с заданными характеристиками фильтра.

Python3

 # Given specification
 Fs = 8000 # Sampling frequency in Hz
 fp = 2000 # Pass band frequency in Hz
 fs = 500 # Stop Band frequency in Hz
 Ap = 3 # Pass band ripple in dB
 As = 20 # Stop band attenuation in dB
 
# Compute Sampling parameter
 Td = 1 / Fs 

Шаг 4: Расчет частоты среза

Python3

 # Compute cut-off frequency in radian/sec
 wp = 2 * np.pi * fp # pass band frequency in radian/sec
 ws = 2 * np.pi * fs # stop band frequency in radian/sec 

Шаг 5: предварительная установка частоты среза

Python3

# Prewarp the analog frequency
Omega_p = (2/Td)*np.tan(wp*Td/2)  # Prewarped analog passband frequency
Omega_s = (2/Td)*np.tan(ws*Td/2)  # Prewarped analog stopband frequency

Шаг 6: вычисление фильтра Баттерворта

Python3

 # Compute Butterworth filter order and cutoff frequency
 N, wc = signal.buttord(Omega_p, Omega_s, Ap, As, analog = True )
 
# Print the values of order and cut-off frequency
 print ( 'Order of the filter=' , N)
 print ( 'Cut-off frequency=' , wc) 

Выход:

Шаг 7. Разработайте аналоговый фильтр Баттерворта, используя N и wc с помощью функции signal.butter ( ).

Python3

 # Design analog Butterworth filter using N and
 # wc by signal.butter function
 b, a = signal.butter(N, wc, 'high' , 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) 

Выход:

Шаг 8: Построение графика амплитуды и фазовой характеристики

Python3

 # Call mfreqz function to plot the
 # magnitude and phase response
 mfreqz(z, p, Fs)

Выход:

Шаг 9: Построение импульсной и ступенчатой характеристики

Python3

 # Call impz function to plot impulse and
 # step response of the filter
 impz(z, p)

Выход:

Ниже представлена реализация:

Python3

 # import required library
 import numpy as np
 import scipy.signal as signal
 import matplotlib.pyplot as plt
 
# User defined functions mfreqz for
 # Magnitude & Phase Response
 def mfreqz(b, a, Fs):
     
    # Compute frequency response of the filter
    # using signal.freqz function
    wz, hz = signal.freqz(b, a)
 
    # Calculate Magnitude from hz in dB
    Mag = 20 * np.log10( abs (hz))
 
    # Calculate phase angle in degree from hz
    Phase = np.unwrap(np.arctan2(np.imag(hz), np.real(hz))) * ( 180 / np.pi)
 
    # Calculate frequency in Hz from wz
    Freq = wz * Fs / ( 2 * np.pi) # START CODE HERE ### (≈ 1 line of code)
 
    # Plot filter magnitude and phase responses using subplot.
    fig = plt.figure(figsize = ( 10 , 6 ))
 
    # Plot Magnitude response
    sub1 = plt.subplot( 2 , 1 , 1 )
    sub1.plot(Freq, Mag, 'r' , linewidth = 2 )
    sub1.axis([ 1 , Fs / 2 , - 100 , 5 ])
    sub1.set_title( 'Magnitute Response' , fontsize = 20 )
    sub1.set_xlabel( 'Frequency [Hz]' , fontsize = 20 )
    sub1.set_ylabel( 'Magnitude [dB]' , fontsize = 20 )
    sub1.grid()
 
    # Plot phase angle
    sub2 = plt.subplot( 2 , 1 , 2 )
    sub2.plot(Freq, Phase, 'g' , linewidth = 2 )
    sub2.set_ylabel( 'Phase (degree)' , fontsize = 20 )
    sub2.set_xlabel(r 'Frequency (Hz)' , fontsize = 20 )
    sub2.set_title(r 'Phase response' , fontsize = 20 )
    sub2.grid()
 
    plt.subplots_adjust(hspace = 0.5 )
    fig.tight_layout()
    plt.show()
 
# Define impz(b,a) to calculate impulse
 # response and step response of a system
 # input: b= an array containing numerator
 # coefficients,a= an array containing
 #denominator coefficients
 def impz(b, a):
     
    # Define the impulse sequence of length 60
    impulse = np.repeat( 0. , 60 )
    impulse[ 0 ] = 1.
    x = np.arange( 0 , 60 )
 
    # Compute the impulse response
    response = signal.lfilter(b, a, impulse)
 
    # Plot filter impulse and step response:
    fig = plt.figure(figsize = ( 10 , 6 ))
    plt.subplot( 211 )
    plt.stem(x, response, 'm' , use_line_collection = True )
    plt.ylabel( 'Amplitude' , fontsize = 15 )
    plt.xlabel(r 'n (samples)' , fontsize = 15 )
    plt.title(r 'Impulse response' , fontsize = 15 )
 
    plt.subplot( 212 )
    step = np.cumsum(response) # Compute step response of the system
    plt.stem(x, step, 'g' , use_line_collection = True )
    plt.ylabel( 'Amplitude' , fontsize = 15 )
    plt.xlabel(r 'n (samples)' , fontsize = 15 )
    plt.title(r 'Step response' , fontsize = 15 )
    plt.subplots_adjust(hspace = 0.5 )
 
    fig.tight_layout()
    plt.show()
 
 
# Given specification
 Fs = 8000 # Sampling frequency in Hz
 fp = 2000 # Pass band frequency in Hz
 fs = 500 # Stop Band frequency in Hz
 Ap = 3 # Pass band ripple in dB
 As = 20 # Stop band attenuation in dB
 
# Compute Sampling parameter
 Td = 1 / Fs
 
# Compute cut-off frequency in radian/sec
 wp = 2 * np.pi * fp # pass band frequency in radian/sec
 ws = 2 * np.pi * fs # stop band frequency in radian/sec
 
# Prewarp the analog frequency
 Omega_p = ( 2 / Td) * np.tan(wp * Td / 2 ) # Prewarped analog passband frequency
 Omega_s = ( 2 / Td) * np.tan(ws * Td / 2 ) # Prewarped analog stopband frequency
 
# Compute Butterworth filter order and cutoff frequency
 N, wc = signal.buttord(Omega_p, Omega_s, Ap, As, analog = True )
 
# Print the values of order and cut-off frequency
 print ( 'Order of the filter=' , N)
 print ( 'Cut-off frequency=' , wc)
 
# Design analog Butterworth filter using N and
 # wc by signal.butter function
 b, a = signal.butter(N, wc, 'high' , 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)
 
# Call mfreqz function to plot the magnitude
 # and phase response
 mfreqz(z, p, Fs)
 
# Call impz function to plot impulse and step
 # response of the filter
 impz(z, p)

Выход:

Внимание компьютерщик! Укрепите свои основы с помощью базового курса программирования Python и изучите основы.

Для начала подготовьтесь к собеседованию. Расширьте свои концепции структур данных с помощью курса Python DS. А чтобы начать свое путешествие по машинному обучению, присоединяйтесь к курсу Машинное обучение - базовый уровень.

Python Визуализация данных

Разработайте БИХ-фильтр верхних частот Баттерворта с использованием метода билинейного преобразования в Scipy - Python

Что такое IIR Highpass Butterworth?

Python3

Python3

Python3

Python3

Python3

Python3

Python3

Python3

Python3

Python3

РЕКОМЕНДУЕМЫЕ СТАТЬИ