Python으로 구현해 보는 디지털 저역통과필터 (1차 Low Pass Filter)

1차 필터는 생각보다 블로그에서 많이 다루었더라구요^^. 처음 1차 저역/고역 통과필터를 C로 구현하는 방법에 대한 이야기[바로가기]때 부터 MATLAB[바로가기]뿐만 아나라 Python에서도 어떻게 구현할 것인지 이야기[바로가기]했지요. 심지어 전 엑셀에서 저역통과필터를 구현하는 것도 이야기[바로가기]를 했던 적이 있습니다.^^. 오늘은 그 대상이 Python입니다만, 실제로는 디지털 필터를 어떻게 구현할 것인지를 한 번 정리하는 것을 목적으로 합니다. 일단 1차 저역통과필터를 대상으로 차단주파수를 결정했을 때, 어떻게 디지털 필터로 변환하며 또 어떻게 실제 코드로 구현할 것인지를 보는 것이 목적입니다.^^.

연속시간 영역에서의 1차 저역통과필터

앞 선 필터 관련 글들에서도 참 자주 나온 수식이지만, 위 수식은 1차 저역통과필터(LPF)의 s-domain 표현식입니다. 위 식의 bode 선도도 이야기[바로가기]한 적이 있습니다. 뭐 아무튼... 만약 차단 주파수(f_cut)를 결정했다면~ 차단 각주파수는...

이 됩니다.~

그러면 tau는 위 식으로 결정됩니다.^^ 이제 차단 주파수를 정하면 tau를 결정할 수 있겠죠^^ 만약 차단 주파수를 100Hz에서 결정했다면~

tau는 위 수식으로 결정이 됩니다.^^

그러면 애초 1차 LPF는 위와 같이 결정이 됩니다.

f_cut = 100
w_cut = 2*np.pi*f_cut
tau = 1/w_cut

num = np.array([1.])
den = np.array([tau, 1.])
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_cut, color='k', lw=1)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude response [dB]')
tmpTitle = 'Low Pass Filter, cutoff freq. at ' + str(f_cut) + 'Hz in continuous'
plt.title(tmpTitle)
plt.xlim(1, Fs/2)
plt.grid()
plt.show()

그러면 위 Python 코드처럼 실행해주면, 보드선도를 그려서 확인할 수 있습니다. 아~ freqs는 scipy.signal 모듈안에 있습니다.

그렇게 해서 위 그림처럼 100Hz가 차단 주파수로 설정된 1차 저역통과필터를 확인할 수 있습니다. 100Hz 지점에서 3dB 떨어지고, -20dB의 기울기로 하강하는 것이 보이네요^^ 이렇게 해서 차단 주파수가 결정되었을 때 tau를 계산하는 법을 보았습니다. 이제 이를 디지털 필터로 한 번 변환해 보도록 하겠습니다.^^

연속시간 영역에서의 1차 저역통과필터를 이산 시간 영역에서 구하기

이제 다시 위 식에서~

위와 같이 정리하고~

라플라스 역변환으로 위와 같이 표현가능합니다. 여기서 미분항을 차분으로 표현하면~

연속시간의 수식이 이산시간에서의 수식으로 편하게 변환됩니다.^^

그러면 이렇게 표현되지요. 여기서 C, MATLAB, Python등등 글 초반에 이야기한 예제들을 보여드렸었는데요^^ 이 상태에서 약간 정형화된 모습을 보기 위해 z-변환을 하도록 하겠습니다.^^

그러면 위 결과를 얻지요. 이 상태에서~ 아까 100Hz에서 구한 tau와 현재 샘플 주파수를 20kHz라고 가정한 후,

num_z = np.array([Ts/(tau+Ts)])
den_z = np.array([1., -tau/(tau+Ts)])

wz, hz = sig.freqz(num_z, den_z, worN=10000)

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

위 코드를 실행하면 얻는 그림이

입니다. 샘플주파수의 절반(Nyquist Frequency)에 가까워질 수록 연속시간과는 차이가 나지만, 일단 나머지 부분에서는 비슷한 결과네요. 아 Python에서 연속시간의 주파수 응답을 얻기 위한 함수가 freqs이고, 이산시간에서는 freqz입니다 그런데.... freqs는 그냥 주파수를 domain으로 출력해서 plot할 때 그냥 사용하면 되는데, freqz느 normalize된 주파수를 내놓기 때문에 plot할 때, Fs(sample freqeuncy)를 곱해야 합니다. 뭐 아무튼 여기까지 해서 라플라스 표현으로 나타내는 1차 저역통과필터를 z-변환의 형태로 표현하도록 공식화 할 수 있네요.

디지털 필터 구현하기

이제 실제 필터의 구현에 대해 이야기를 해보겠습니다. 필터의 형태에서 아마도~ 가장 많이 사용하는 형태가

2차의 Direct II Form

이런 Direct II Form일텐데요. 이 형태에 맞춰서 계수를 정리할 수 있습니다.

이렇게 말이죠^^ 위 그림의 2차 Direc II Form의 경우는 def로 만들어 두면 편할 겁니다.^^.

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 filterSimpleLPF(data, tau, Ts):
    from numpy import zeros, arange
    
    result = zeros((len(data),))
    
    for n in arange(1, len(data)):
        result[n] = (tau*result[n-1] + Ts*data[n])/(tau + Ts)
        
    return result

그게 위 함수인 direc2FormModel입니다. 블럭을 그대로 그냥 뭐 아무 생각없이 작성했다고 생각하시면 됩니다^^. 이왕하는 김에 이전 글에서 tau로 표현하던 간편한 형태를 또 filterSimpleLPF로 작성해 두었습니다. 또 direct2form이 아니라 바로 차분 방정식으로 표현하는 것도 나중을 위해 2차로 구현을 해 두었습니다.(differentialEqForm)

이제 제가 맨날 하듯이~~~ 시험신호 하나 만들고 가죠^^

# 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*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()

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

2Hz의 주 주파수를 두고, 10부터 500까지 50간격으로 주파수를 만들어 더해둡니다. 아 그리고 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()

입니다. 이제~ 시험신호는

이렇습니다. 흠~~ 이쁘네요^^

주파수 응답은 또 위와 같습니다.^^. 원했던 주파수가 다 표현되어 있네요^^

filteredSig1 = filterSimpleLPF(inputSig, tau, Ts)
filteredSig2 = differentialEqForm(inputSig, a1, a2, b0, b1, b2)
filteredSig3 = direct2FormModel(inputSig, a1, a2, b0, b1, b2)

draw_FFT_Graph(filteredSig1, Fs, title='filterSimpleLPF', xlim=(0, 500))
draw_FFT_Graph(filteredSig2, Fs, title='differentialEqForm', xlim=(0, 500))
draw_FFT_Graph(filteredSig3, Fs, title='direct2FormModel', xlim=(0, 500))

이제 간편한 제가 자주 블로그에서 언급하던 방법(filterSimpleLPF)과 일반적으로 필터를 구현할 때 많이 사용하는 Direct II Form, 그리고 차분방정식으로 그대로 표현하는 differentialEqForm을 모두 적용한 후 FFT로 확인했더니....

모두 위와 같이 나왔습니다 큰 차이가 없다는 거겠죠^^

그 결과도 위와 같이 나타납니다. LPF가 잘 적용된 듯 보입니다.^^

각 방법별 에러를 확인했는데 아주 작거나 혹은 첫 샘플에서만 약간의 에러를 가지는 것으로 확인됩니다. 괜찮네요^^

오늘은

차단 주파수를 결정했을 때, 연속시간 영역에서 1차 저역통과필터 결정하기
연속시간 영영에서 결정된 1차 저역통과필터를 이산 시간영역에서 표현하기
결정된 필터를 알려진 몇몇 방법으로 직접 구현해보기

를 보여드리고 싶었습니다.^^ 마지막으로 전체 Python 코드는 아래에 있습니다.^^

# -*- 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 filterSimpleLPF(data, tau, Ts):
    from numpy import zeros, arange
    
    result = zeros((len(data),))
    
    for n in arange(1, len(data)):
        result[n] = (tau*result[n-1] + Ts*data[n])/(tau + Ts)
        
    return result
    
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()

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*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()

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

# Design 1st Order Low Pass Filter
f_cut = 100
w_cut = 2*np.pi*f_cut
tau = 1/w_cut

num = np.array([1.])
den = np.array([tau, 1.])
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_cut, color='k', lw=1)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude response [dB]')
tmpTitle = 'Low Pass Filter, cutoff freq. at ' + str(f_cut) + 'Hz in continuous'
plt.title(tmpTitle)
plt.xlim(1, Fs/2)
plt.grid()
plt.show()

# Design 1st Order Low Pass Filter Z-Transform
num_z = np.array([Ts/(tau+Ts)])
den_z = np.array([1., -tau/(tau+Ts)])

wz, hz = sig.freqz(num_z, den_z, worN=10000)

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

# Implementation Signal
a1 = den_z[1]
a2 = 0
b0 = num_z[0]
b1 = 0.
b2 = 0.

filteredSig1 = filterSimpleLPF(inputSig, tau, Ts)
filteredSig2 = differentialEqForm(inputSig, a1, a2, b0, b1, b2)
filteredSig3 = direct2FormModel(inputSig, a1, a2, b0, b1, b2)

draw_FFT_Graph(filteredSig1, Fs, title='filterSimpleLPF', xlim=(0, 500))
draw_FFT_Graph(filteredSig2, Fs, title='differentialEqForm', xlim=(0, 500))
draw_FFT_Graph(filteredSig3, Fs, title='direct2FormModel', xlim=(0, 500))

plt.figure(figsize=(12,5))
plt.plot(t, inputSig)
plt.plot(t, filteredSig1, 'r', label='SimpleLPF')
plt.plot(t, filteredSig2, 'c', label='DE')
plt.plot(t, filteredSig3, 'k', label='D2F')
plt.xlabel('Time(s)')
plt.title('Filtered Signal in Discrete')
plt.legend()
plt.grid(True)
plt.show()

# Err
err = filteredSig2 - filteredSig1
plt.figure(figsize=(12,5))
plt.plot(t, err)
plt.xlabel('Time(s)')
plt.title('Error between differentialEqForm and filterSimpleLPF')
plt.grid(True)
plt.show()

err = filteredSig2 - filteredSig3
plt.figure(figsize=(12,5))
plt.plot(t, err)
plt.xlabel('Time(s)')
plt.title('Error between differentialEqForm and direct2FormModel')
plt.grid(True)
plt.show()

err = filteredSig1 - filteredSig3
plt.figure(figsize=(12,5))
plt.plot(t, err)
plt.xlabel('Time(s)')
plt.title('Error between filterSimpleLPF and direct2FormModel')
plt.grid(True)
plt.show()

'Software > Python' 카테고리의 다른 글

Python으로 구현해 보는 Band Stop Filter (6)	2016.09.22
Python으로 구현해보는 Band Pass Filter (18)	2016.09.08
Python으로 구현해보는 1차 고역통과필터 (8)	2016.09.01
Python에서 보드 선도 Bode Plot 그려보고 그래프 있는 척 치장하기^^ (6)	2016.08.24
IPython Notebook에서 Markdown 사용하기... (6)	2016.08.19
Anaconda를 이용한 Python 설치와 Spyder의 편리함 (12)	2016.08.05
Python에서 그래프를 애니메이션으로 표현하고 GIF로 저장하기 (18)	2016.04.01

PinkWink

Python으로 구현해 보는 디지털 저역통과필터 (1차 Low Pass Filter)

연속시간 영역에서의 1차 저역통과필터

연속시간 영역에서의 1차 저역통과필터를 이산 시간 영역에서 구하기

디지털 필터 구현하기

'Software > Python' 카테고리의 다른 글

티스토리툴바

Python으로 구현해 보는 디지털 저역통과필터 (1차 Low Pass Filter)

연속시간 영역에서의 1차 저역통과필터

연속시간 영역에서의 1차 저역통과필터를 이산 시간 영역에서 구하기

디지털 필터 구현하기

'Software > Python' 카테고리의 다른 글

'Software/Python' Related Articles

티스토리툴바