Python으로 구현해보는 Band Pass Filter

MATLAB에서 보통 많이 하는 필터 설계나 확인을 Python으로도 할 수 있다는 걸 살짝 보여줄려고 시작한 글이 이제 세 번째네요. 처음 1차 저역통과필터[바로가기]였구요. 그 다음 1차 고역통과필터[바로가기]였습니다. 이제 이번에는 Band Pass Filter를 이야기할려고 합니다.

연속시간 영역에서의 Band Pass Filter

연속시간 영역의 s-domain에서 표현된 Band Pass Filter는

와 같이 2차로 나타납니다. 그 중에서 분자(num)에 s 일차항이 있으면 band pass입니다. 여러가지 접근법이 있지만, 위 수식처럼 표현하는 방법이 그 중 하나입니다.

Band pass filter의 경우 cut-off frequency라고 하지 않고, 통과시킬 band의 폭의 가운데를 peak frequency라고 합니다. 그 peak frequency가 정해지만, omega_0가 정해집니다.^^

그리고, 대역폭 band width는 peak frequency를 중심으로 통과시키는 주파수의 폭을 의미하며 위 수식으로 정해지고,

Q를 위와 같이 계산할 수 있습니다.

그리고, H와 H-0를 계산하면 최초 수식을 완성할 수 있습니다.^^

만약 peak frequency를 500Hz로 두고, 대역폭을 200Hz로 두면 위와 같을 결과를 얻을 수 있습니다.

그러면 위 수식이 만들어지고,

보드 선도를 그려보면 위와 같이 나타납니다. 500Hz가 peak이고, 그 양옆으로 100Hz씩 빼고 더한 곳에서 -3dB 떨어지게 됩니다. 위 보드 선도를 얻을 떄까지의 Python 코드는

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

# Band Pass Filter in Continuous Time
f_peak = 500
w0_peak = 2*np.pi*f_peak
bandWidth = 200
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])

w, h = sig.freqs(num, den, worN=np.logspace(0, 5, 1000))

plt.figure(figsize=(12,5))
plt.semilogx(toHz(w), 20 * np.log10(abs(h)))
plt.axvline(f_peak, color='k', lw=2)
plt.axvline(f_peak-bandWidth/2, color='k', lw=1)
plt.axvline(f_peak+bandWidth/2, color='k', lw=1)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude response [dB]')
tmpTitle = 'Band Pass Filter, Peak freq. at ' + str(f_peak) + 'Hz in continuous'
plt.title(tmpTitle)
plt.xlim((f_peak-bandWidth/2)*0.1, (f_peak+bandWidth/2)*10)
plt.grid()
plt.show()

입니다. num과 den을 계산한 후 freqs함수에서 주파수와 크기를 얻은 다음 그린거죠^^

연속시간 영역에서의 Band Pass Filter를 이산 시간 영역에서 구하기

이제 앞 선 장에서 이야기한 연속 시간 영역에서 보여준 Band Pass Filter를 이산 시간영역에서 표현하기 위해서는 z-domain으로 옮길 필요가 있습니다. 간편한 공식으로는 

를 대입하는 것입니다.

그러면 위 수식과 같은 형태로 표현이 될 겁니다. 여기서 각 계수는

위와 같습니다. 음~~ 연습장에서 좀 고생해야 합니다.^^ 이제 분모의 상수항을 1로 두는

형태로 표현하면, 각 계수는

a_0가 1인 형태로 표현되겠죠^^ 아까 했던 500Hz peak에 200Hz band의 필터를 적용하면....

요런 결과가 나타납니다.^^

그 결과에 대한 Bode 선도는 위와 같습니다. 고주파로 가면 20kHz로 샘플링을 설정했으므로 차이가 나지만, 관심 주파수 대역에서는 문제가 없습니다. discrete가 두개 인것은 위에 이야기한 손으로 푼 공식과 Python scipy가 제공하는 bilinear라는 함수와 결과를 비교한 겁니다. 제가 맞는지는 알아야죠^^ 위 보드 선도를 얻은 코드는...

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

numz, denz = sig.bilinear(num, den , Fs)

wz, hz = sig.freqz(numz, denz, worN=10000)

numz1 = np.array([2*H0*w0_peak**2/Ts, 0, 
                  -2*H0*w0_peak**2/Ts])
denz1 = np.array([4/Ts**2+2*w0_peak/Q/Ts+w0_peak**2, 
                  -8/Ts**2+2*w0_peak**2, 
                  4/Ts**2-2*w0_peak/Q/Ts+w0_peak**2])

numz1 = numz1/denz1[0]
denz1 = denz1/denz1[0]

wz1, hz1 = sig.freqz(numz1, denz1, worN=10000)

plt.figure(figsize=(12,5))
plt.semilogx(toHz(w), 20 * np.log10(abs(h)), label='continuous')
plt.semilogx(toHz(wz*Fs), 20 * np.log10(abs(hz)), 'c', label='discrete')
plt.semilogx(toHz(wz1*Fs), 20 * np.log10(abs(hz1)), 'r', label='discrete')
plt.axvline(f_peak, color='k', lw=2)
plt.axvline(f_peak-bandWidth/2, color='k', lw=1)
plt.axvline(f_peak+bandWidth/2, color='k', lw=1)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude response [dB]')
tmpTitle = 'Band Pass Filter, Peak freq. at ' + str(f_peak) + 'Hz in discrete'
plt.title(tmpTitle)
plt.xlim((f_peak-bandWidth/2)*0.02, (f_peak+bandWidth/2)*50)
plt.legend()
plt.grid()
plt.show()

입니다.^^

디지털 필터 구현하기

이제 필터를 적용해볼까요.. 그런데 이 부분에서는 쓸 이야기가 없습니다, 이미 [바로가기]와 [바로가기]에서 사용한 코드를 그대로 사용할 것이니까요.

def direct2FormModel(data, a1, a2, b0, b1, b2):
    from numpy import zeros, arange
    
    result = zeros((len(data),))
    timeZone = zeros((len(data),))
    
    for n in arange(2, len(data)):
        sum0 = -a1*timeZone[n-1] - a2*timeZone[n-2]
        timeZone[n] = data[n] + sum0
        result[n] = b0*timeZone[n] + b1*timeZone[n-1] + b2*timeZone[n-2]
        
    return result
    
def differentialEqForm(data, a1, a2, b0, b1, b2):
    from numpy import zeros, arange
    
    result = zeros((len(data),))
    
    for n in arange(2, len(data)):
        result[n] = b0*data[n] + b1*data[n-1] + b2*data[n-2] - a1*result[n-1] - a2*result[n-2]
        
    return result

위 함수 둘 중 하나를 사용하면 됩니다^^.

filteredOut = direct2FormModel(inputSig, denz[1], denz[2], numz[0], numz[1], numz[2])

적용은 위 코드처럼 한 줄이지요^^

아무튼 이런 시험 신호를 만들어서~

FFT를 해본 결과 500Hz에 약간 비중을 두고 만들었습니다.^^

필터를 적용하면 저주파와 고주파 대역은 살짝 사라지고, 제가 원한 부분만 이렇게 남게 되네요^^ 아~~ 전체 코드는 아래에 둡니다^^

# -*- coding: utf-8 -*-

def toHz(value):
    from numpy import pi
    return value/2/pi
    
def direct2FormModel(data, a1, a2, b0, b1, b2):
    from numpy import zeros, arange
    
    result = zeros((len(data),))
    timeZone = zeros((len(data),))
    
    for n in arange(2, len(data)):
        sum0 = -a1*timeZone[n-1] - a2*timeZone[n-2]
        timeZone[n] = data[n] + sum0
        result[n] = b0*timeZone[n] + b1*timeZone[n-1] + b2*timeZone[n-2]
        
    return result
    
def differentialEqForm(data, a1, a2, b0, b1, b2):
    from numpy import zeros, arange
    
    result = zeros((len(data),))
    
    for n in arange(2, len(data)):
        result[n] = b0*data[n] + b1*data[n-1] + b2*data[n-2] - a1*result[n-1] - a2*result[n-2]
        
    return result
    
def draw_FFT_Graph(data, fs, **kwargs):
    from numpy.fft import fft
    
    graphStyle = kwargs.get('style', 0)
    xlim = kwargs.get('xlim', 0)
    ylim = kwargs.get('ylim', 0)
    title = kwargs.get('title', 'FFT result')
    
    n = len(data)
    k = np.arange(n)
    T = n/fs
    freq = k/T 
    freq = freq[range(int(n/2))]
    FFT_data = fft(data)/n 
    FFT_data = FFT_data[range(int(n/2))]
    
    plt.figure(figsize=(12,5))
    if graphStyle == 0:
        plt.plot(freq, abs(FFT_data), 'r', linestyle=' ', marker='^') 
    else:
        plt.plot(freq,abs(FFT_data),'r')
    plt.xlabel('Freq (Hz)')
    plt.ylabel('|Y(freq)|')
    plt.vlines(freq, [0], abs(FFT_data))
    plt.title(title)
    plt.grid(True)
    plt.xlim(xlim)
    plt.ylim(ylim)
    plt.show()

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

# Create Test Signal
Fs = 20*10**3               # 20kHz
Ts = 1/Fs                   # sample Time
endTime = 2
t = np.arange(0.0, endTime, Ts)

inputSig = 3.*np.sin(2.*np.pi*500.*t)

sampleFreq = np.arange(50,1000,50)

for freq in sampleFreq:
    inputSig = inputSig + 2*np.sin(2*np.pi*freq*t)
    
plt.figure(figsize=(12,5))
plt.plot(t, inputSig)
plt.xlabel('Time(s)')
plt.title('Test Signal in Continuous')
plt.grid(True)
plt.show()

draw_FFT_Graph(inputSig, Fs, title='inputSig', xlim=(0, 1200))

# Band Pass Filter in Continuous Time
f_peak = 500
w0_peak = 2*np.pi*f_peak
bandWidth = 200
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])

w, h = sig.freqs(num, den, worN=np.logspace(0, 5, 1000))

plt.figure(figsize=(12,5))
plt.semilogx(toHz(w), 20 * np.log10(abs(h)))
plt.axvline(f_peak, color='k', lw=2)
plt.axvline(f_peak-bandWidth/2, color='k', lw=1)
plt.axvline(f_peak+bandWidth/2, color='k', lw=1)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude response [dB]')
tmpTitle = 'Band Pass Filter, Peak freq. at ' + str(f_peak) + 'Hz in continuous'
plt.title(tmpTitle)
plt.xlim((f_peak-bandWidth/2)*0.1, (f_peak+bandWidth/2)*10)
plt.grid()
plt.show()

numz, denz = sig.bilinear(num, den , Fs)

wz, hz = sig.freqz(numz, denz, worN=10000)

numz1 = np.array([2*H0*w0_peak**2/Ts, 0, 
                  -2*H0*w0_peak**2/Ts])
denz1 = np.array([4/Ts**2+2*w0_peak/Q/Ts+w0_peak**2, 
                  -8/Ts**2+2*w0_peak**2, 
                  4/Ts**2-2*w0_peak/Q/Ts+w0_peak**2])

numz1 = numz1/denz1[0]
denz1 = denz1/denz1[0]

wz1, hz1 = sig.freqz(numz1, denz1, worN=10000)

plt.figure(figsize=(12,5))
plt.semilogx(toHz(w), 20 * np.log10(abs(h)), label='continuous')
plt.semilogx(toHz(wz*Fs), 20 * np.log10(abs(hz)), 'c', label='discrete')
plt.semilogx(toHz(wz1*Fs), 20 * np.log10(abs(hz1)), 'r', label='discrete')
plt.axvline(f_peak, color='k', lw=2)
plt.axvline(f_peak-bandWidth/2, color='k', lw=1)
plt.axvline(f_peak+bandWidth/2, color='k', lw=1)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude response [dB]')
tmpTitle = 'Band Pass Filter, Peak freq. at ' + str(f_peak) + 'Hz in discrete'
plt.title(tmpTitle)
plt.xlim((f_peak-bandWidth/2)*0.02, (f_peak+bandWidth/2)*50)
plt.legend()
plt.grid()
plt.show()

filteredOut = direct2FormModel(inputSig, denz[1], denz[2], numz[0], numz[1], numz[2])

draw_FFT_Graph(filteredOut, Fs, xlim=(0, 1200))

plt.figure(figsize=(12,5))
plt.plot(t, inputSig, label='inputSig')
plt.plot(t, filteredOut, 'r', label='filteredOut')
plt.xlabel('Time(s)')
plt.legend()
plt.grid(True)
plt.show()