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

전 아주 예전부터 제 블로그에 필터에 관한 글을 올리고 있었습니다. 물론 단편적인 글들이었지만요. 그 중에는 실제 자이로센서와 가속도센서를 융합하는 상보필터[바로가기]를 다루기 시작했고, 그 후 MATLAB이나 Python으로 구현하는 1차 혹은 2차 필터들에 대한 글들을 올렸습니다. 물론 이전에 올린 글들이 많지만, 문득 필터들에 대한 이야기를 1차 저역/고역 통과필터, 2차 필터, Band Pass, Band Stop 필터등등에 대해 한 번 쫘~악 정리하고 싶어지더라구요. 그래서 이런 것들을 정리할 겸 필터를 대상으로 연재를 시작할까 합니다.^^

Filter

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

Filter

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

Filter

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

오늘 다룰 이야기는 1차 저역 통과필터와 1차 고역 통과필터 흔히 LPF(Low Pass Filter), HPF(High Pass Filter)에 대한 이야기 입니다.

흔히 Hz 단위의 차단주파수(f_cut)가 결정되면 위의 간단한 변환으로 차단 각주파수가 계산됩니다.

이를 1차식으로 s-domain에서 표현하면 LPF는 위와 같이...

HPF는 위와 같이 표현됩니다.

만약 흔히 알려진 시정수(tau)로 LPF와 HPF를 표현하고 싶다면... 위 관계식을 이용해서

이렇게 표현할 수도 있습니다. 뭐 여하튼.. 1차의 저역통과필터식과 고역통과필터식은 합하면 '1'이 됩니다.^^.

1차 LPF와 HPF의 Bode선도를 그려보면 위와 같습니다.

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

def toHz(value):
    from numpy import pi
    return value/2/pi

f_cut = 1000
w_cut = 2*np.pi*f_cut

num_L1 = np.array([w_cut])
den_L1 = np.array([1., w_cut])
s_L1 = sig.lti(num_L1, den_L1)
w_L1, m_L1, P_L1 = sig.bode(s_L1)

num_H1 = np.array([1., 0.])
den_H1 = np.array([1., w_cut])
s_H1 = sig.lti(num_H1, den_H1)
w_H1, m_H1, P_H1 = sig.bode(s_H1)

plt.figure(figsize=(12,5))
plt.semilogx(toHz(w_L1), m_L1, lw=2, label='1st LPF')
plt.semilogx(toHz(w_H1), m_H1, lw=2, label='1st HPF')
plt.axvline(f_cut, color='k', lw=1)
plt.xlim(50, 15000)
plt.ylim(-25, 2)
plt.ylabel('Amplitude [dB]')
plt.xticks([100, 1000, 10000], ('','$f_{cut}$[Hz]',''), fontsize = 20)
plt.legend()
plt.grid()

plt.figure(figsize=(12,5))
plt.semilogx(toHz(w_L1), P_L1, lw=2, label='1st LPF')
plt.semilogx(toHz(w_H1), P_H1, lw=2, label='1st HPF')
plt.axvline(f_cut, color='k', lw=1)
plt.xlim(50, 15000)
plt.ylabel('Phase [degree]')
plt.legend()
plt.xticks([100, 1000, 10000], ('','$f_{cut}$[Hz]',''), fontsize = 20)
plt.grid()

plt.show()

f_cut 지점에서 -3dB가 떨어지는 것을 확인할 수 있습니다.

이제 이를 디지털 필터, Digital Filter로 변환하는 과정을 한 번 보도록 하겠습니다. 위 수식처럼 애초 s-domain에서의 LPF를 이산영역에서 표현할 수 있습니다. 위 과정의 상세한 전개는 [바로가기]에서 다루었습니다.

위 식은 HPF의 이산영역 수식입니다.

식 전개의 편의성을 위해 tau_ts를 위와 같이 두고, LPF와 HPF의 이산영역에서의 수식을 다시 편집하고 이 수식을 z-domain에서 표현하면...

와 같습니다. 즉, 설정하고 싶은 주파수 f_cut을 이용해서 알 수 있는 tau와 이산영역이 되면서 당연히 생겨난 ts(sample time)을 공식으로 위 수식을 완성할 수 있을 겁니다.

일반적으로 biquad라고 알려진 Direct2Form에 일반화해서 적용하기 위해 위의 수식을 사용하겠습니다.

필터는 위와 같은 모양이 될 거구요^^

그러면 Low Pass Filter의 계수는 위의 공식으로~~~

High Pass Filter의 계수는 위의 공식으로 결정됩니다. 이제 digital 영역에서 손쉽게 1차 저역/고역 통과필터의 계수를 찾을 수 있겠죠^^ 여기서... 

tau_ts는 위의 식입니다.

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

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

sampleFreq = np.arange(10,500,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()

이제 시험신호를 만들어서 한 번 적용해 봐야죠^^& 위 시험신호는 10kHz 샘플 주파수에서 생성된 것으로 가정했구요. 기본적으로 2Hz 주 성분에 10Hz부터 500Hz까지 50Hz 간격으로 신호를 생성해서 합했습니다. 이 신호의 FFT결과를 보면

위와 같이 의도한 대로 깔끔하게 나타나네요^^

def draw_FFT_Graph(data, fs, **kwargs):
    from numpy.fft import fft
    import matplotlib.pyplot as plt
    
    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()

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

앞부분에서 사용한 대로 공식화하고 이 필터의 계수를 찾고 적용하는 부분에 대해 Python으로 코드를 확인해 보겠습니다.

def get1stFilterCoeffi(f_cut, ts, isLPF):
    from numpy import pi
    
    w_cut = 2*pi*f_cut
    tau = 1/w_cut
    
    tau_ts = 1/(tau+ts)
    if isLPF=='LPF':
        return -tau*tau_ts, ts*tau_ts, 0
    elif isLPF=='HPF':
        return -tau*tau_ts, tau*tau_ts, -tau*tau_ts

먼저 위 코드는 1차 LPF나 HPF의 계수를 direct2form에 맞춰 계산하는 코드입니다. 뭐 크게 어렵지 않습니다.^^.

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 d2f_1st(data, a1, b0, b1):
    return direct2FormModel(data, a1, 0, b0, b1, 0)

먼저 direct2FormModel 함수는 2차 필터까지 고려한 것으로 d2f_1st라는 함수로 2차에서만 사용하는 계수는 0을 넣은 것입니다. 아무튼... 위 함수는 direct2form에 맞춰 필털를 동작시키는 함수입니다. 계수는 아까 get1stFilterCoeffi함수로 찾았으니까 적용할 수 있는거죠^^...

이제 이 함수들을 이용해서 LPF를 적용해보면

이런 결과가 나타납니다. 

a1, b0, b1 = get1stFilterCoeffi(200, Ts, 'LPF')
dataLPF = d2f_1st(inputSig, a1, b0, b1)

위 저역통과필터의 결과는 이미 만들어둔 함수로 인해 위 코드 딱 두줄이네요^^ 200Hz를 차단주파수로 가지게 해서 필터 계수를 얻고 이를 적용해서 데이터를 얻었네요.

plt.figure(figsize=(12,5))
plt.plot(t, inputSig, 'r', label='inputSig')
plt.plot(t, dataLPF, lw=2, label='filtered')
plt.xlabel('Time(s)')
plt.title('Test Signal in Continuous')
plt.legend()
plt.grid(True)
plt.show()

실제 FFT 결과도 LPF에 어울리게 나타납니다.

어때요? 그렇죠?^^ 실제 200Hz 다음에 나오는 값이 대략 0.7쯤인데 이를 20log10을 적용시켜보면 -3dB가 되거든요^^. 

draw_FFT_Graph(dataLPF, Fs, title='dataLPF', xlim=(0, 600))

이제 HPF를 적용한 경우를 볼까요?

a1, b0, b1 = get1stFilterCoeffi(200, Ts, 'HPF')
dataHPF = d2f_1st(inputSig, a1, b0, b1)

위와 같이 HPF라고 알려주고 계수를 받아 적용해보면~

요렇게 나타납니다. 저주파 영역이 사라졌기 때문에 필터의 결과가 센터(0)에 딱 붙어 있네요^^

plt.figure(figsize=(12,5))
plt.plot(t, inputSig, 'r', label='inputSig')
plt.plot(t, dataHPF, lw=2, label='filtered')
plt.xlabel('Time(s)')
plt.title('Test Signal in Continuous')
plt.legend()
plt.grid(True)
plt.show()

그리고... FFT의 결과도 예상대로 입니다.^^

draw_FFT_Graph(dataHPF, Fs, title='dataLPF', xlim=(0, 600))

이렇게 해서 시간영역에서 1차 저역통과필터와 1차 고역통과필터의 개념을 이야기했고, 이를 digital로 구현하기 위해 각 계수를 찾는 공식을 도출하고 이를 code로 표현하는 것 까지 다루어 보았습니다. 실제 예시로 Python을 사용했구요^^ 오늘 글 전체에서 다룬 Python의 전체 예시는 Github[바로가기]에 두었습니다....