본문 바로가기

Theory/ControlTheory

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

전 아주 예전부터 제 블로그에 필터에 관한 글을 올리고 있었습니다. 물론 단편적인 글들이었지만요. 그 중에는 실제 자이로센서와 가속도센서를 융합하는 상보필터[바로가기]를 다루기 시작했고, 그 후 MATLAB이나 Python으로 구현하는 1차 혹은 2차 필터들에 대한 글들을 올렸습니다. 물론 이전에 올린 글들이 많지만, 문득 필터들에 대한 이야기를 1차 저역/고역 통과필터, 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선도를 그려보면 위와 같습니다.

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는 위의 식입니다.

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

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

앞부분에서 사용한 대로 공식화하고 이 필터의 계수를 찾고 적용하는 부분에 대해 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차 LPFHPF의 계수를 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를 차단주파수로 가지게 해서 필터 계수를 얻고 이를 적용해서 데이터를 얻었네요.

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

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

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

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

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

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

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

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

반응형