본문으로 바로가기

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

신고