[필터연재] 2차 디지털 Band Pass, Band Stop 필터

요즘 필터 연재를 진행중이네요. 이제 이 글이 발행되고 나면 남은 글은 두 개?^^ 오늘은 2차 디지털 Band Pass Filter와 Notch 필터로 알려진 Band Stop Filter인데요. 근데 연재를 처음 시작하기 전에 Python에서 Filter를 어떻게 다룰것인가라는 생각에 정리했던 글이 이미 있답니다. Band Pass Filter[바로가기], Band Stop Filter[바로가기]. 그 글에서 이미 다루었던 내용을 하나로 합치고, 연재의 전체 흐름에 맞게 살짝꿍 조절만하는 것이 이번 글이 되겠네요^^

Filter

[필터연재] 1차 디지털 저역/고역 통과필터

Filter

[필터연재] 2차 디지털 저역/고역 통과필터

Filter

[필터연재] 2차 디지털 Band Pass, Band Stop 필터

먼저 Band Pass Filter부터 이야기를 해야겠네요^^

먼저 Band의 Center값을 f_peak라고 보고...

각주파수 성분을 바꾼것을 w_0로 하죠~~

그리고... peak 주파수의 양 폭을 bandWidth로 설정하면 Q를 구할 수 있습니다.

H를 계산하고...

Q와 H로 H_0까지 계산해 두면~

위 식으로 2차 Band Pass Filter의 s-domain에서의 표현식이 됩니다.^^ Band Pass Filter의 Bode 선도는

이렇게 됩니다. Band의 양끝에서 -3dB 떨어지는 위로 볼록한 모양을 가지게 되죠.

import scipy.signal as sig
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

# Band Pass Filter in Continuous Time
f_peak = 1000
w0_peak = 2*np.pi*f_peak
bandWidth = 400
Q = f_peak/bandWidth
H = 1/w0_peak
H0 = H/Q

num = np.array([H0*w0_peak**2, 0])
den = np.array([1, w0_peak/Q, w0_peak**2])

s = sig.lti(num, den)
w, m, p = sig.bode(s, n=10000)


plt.figure(figsize=(12,5))
plt.semilogx(toHz(w), m, lw=2, label='2nd BPF')
plt.axvline(f_peak+bandWidth/2, color='k', lw=1)
plt.axvline(f_peak, color='k', lw=1)
plt.axvline(f_peak-bandWidth/2, color='k', lw=1)
plt.xlim(10**2, 10**4)
plt.ylim(-50, 10)
plt.ylabel('Amplitude [dB]')
plt.xticks([f_peak+bandWidth/2, f_peak, f_peak-bandWidth/2], ('','$f_{peak}$[Hz]',''), fontsize = 15)
plt.legend()
plt.grid()
plt.show()

Band Stop Filter(혹은 Notch 필터)는 H_0를

1로 하고...

위와 같이 표현가능합니다. 음... 위 식에서 좌변의 L_BPF가 잘못 되었네요... L_BSF인데..ㅠㅠ. 뭐 아무튼 위 식은 Band Stop Filter입니다. Bode 선도는 

위와 같이 생겼습니다.^^

import scipy.signal as sig
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

# Band Stop Filter in Continuous Time
f_peak = 1000
w0_peak = 2*np.pi*f_peak
bandWidth = 400
Q = f_peak/bandWidth
H0 = 1

num = np.array([H0,0,H0*w0_peak**2])
den = np.array([1, w0_peak/Q, w0_peak**2])

s = sig.lti(num, den)
w, m, p = sig.bode(s, n=10000)


plt.figure(figsize=(12,5))
plt.semilogx(toHz(w), m, lw=2, label='2nd BPF')
plt.axvline(f_peak+bandWidth/2, color='k', lw=1)
plt.axvline(f_peak, color='k', lw=1)
plt.axvline(f_peak-bandWidth/2, color='k', lw=1)
plt.xlim(2*10**2, 5000)
plt.ylim(-60, 10)
plt.ylabel('Amplitude [dB]')
plt.xticks([f_peak+bandWidth/2, f_peak, f_peak-bandWidth/2], ('','$f_{peak}$[Hz]',''), fontsize = 15)
plt.legend()
plt.grid()
plt.show()

여기서 bilinear transform을 이용하기 위해 s-domain의 s에 위 식을 대입하여 정리하면~~

Band Pass Filter의 z-domain에서의 표현식은 위와 같습니다.~

Band Stop Filter의 z-domain에서의 표현식은 위와 같구요~^^

위 그림의 Direct2Form에 적용하기 위한

2차 z-domain 방정식의 표현을 사용한다면... Band Pass Filter의 경우는

위와 같이 계수를 계산할 수 있고... Band Stop Filter의 경우는

의 공식으로 계산할 수 있습니다. 이를 계산하는 Python 코드는 Band Pass Filter의 경우는

def getBPFCoeffi(f_peak, bandWidth, Ts):
    from numpy import pi
    
    w0_peak = 2*np.pi*f_peak
    
    Q = f_peak/bandWidth
    H = 1/w0_peak
    H0 = H/Q
    
    b0_ = 2*H0*w0_peak**2/Ts
    b1_ = 0
    b2_ = -2*H0*w0_peak**2/Ts
    
    a0_ = 4/Ts**2+2*w0_peak/Q/Ts+w0_peak**2
    a1_ = -8/Ts**2+2*w0_peak**2
    a2_ = 4/Ts**2-2*w0_peak/Q/Ts+w0_peak**2
    
    return a1_/a0_, a2_/a0_, b0_/a0_, b1_/a0_, b2_/a0_

으로 Band Stop Filter의 경우는

def getBSFCoeffi(f_peak, bandWidth, Ts):
    from numpy import pi
    
    w0_peak = 2*np.pi*f_peak
    Q = f_peak/bandWidth
    H0 = 1
    
    b0_ = H0*4/Ts**2 + H0*w0_peak**2
    b1_ = -2*H0*4/Ts**2 + 2*w0_peak**2
    b2_ = H0*4/Ts**2 + H0*w0_peak**2
    
    a0_ = 4/Ts**2+2*w0_peak/Q/Ts+w0_peak**2
    a1_ = -8/Ts**2+2*w0_peak**2
    a2_ = 4/Ts**2-2*w0_peak/Q/Ts+w0_peak**2
    
    return a1_/a0_, a2_/a0_, b0_/a0_, b1_/a0_, b2_/a0_

으로 작성할 수 있습니다. 이제 이전 글에서도 사용한 시험신호...

에 적용해 볼까 합니다.

시험신호가 가지는 주파수 성분이었구요~~ 위 시험신호에 BPF를 적용하기 위해

a1, a2, b0, b1, b2 = getBPFCoeffi(200, 100, Ts)
dataBPF = direct2FormModel(inputSig, a1, a2, b0, b1, b2)
draw_FFT_Graph(dataBPF, Fs, title='dataBPF', xlim=(0, 600))

위 코드를 적용해서 얻은 결과가

입니다. 또 BSF에 적용하기 위해

a1, a2, b0, b1, b2 = getBSFCoeffi(200, 100, Ts)
dataBSF = direct2FormModel(inputSig, a1, a2, b0, b1, b2)
draw_FFT_Graph(dataBSF, Fs, title='dataBSF', xlim=(0, 600))

작성한 위 코드의 결과는

입니다. 두 결과 모두 Band Pass/Stop이 잘 적용되었네요^^