Python으로 구현해 보는 Band Stop Filter

어쩌다 보니 Python으로 필터를 구현하는 법에 대해 꽤 많은 글을 적고 있네요^^. 그 시작은 아마, 1차 저역 통과 필터(1st order Low Pass Filter)[바로가기]일겁니다. 그리고 1차 고역통과필터(1st order High Pass Filter)[바로가기]를 이야기했구요. 그리고 나서 2차 Band Pass Filter[바로가기]를 이야기했습니다. 그러나 이제 Band Stop Filter로 오는건 당연한 일이겠죠^^. 특정한 Band만 저지(stop)하는 필터로 Notch라고 불리기도 합니다.

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

먼저

연속시간 영역에서 s-domain으로 표현된 Band Stop Filter는 위 식과 같이 표현됩니다.

stop시킬 구간의 peak를 f_peak로 하면 w0는 위 식으로 구해집니다.

그리고, Q도 역시 계산이 되구요. 여기서 bandWidth는 3dB 떨어지는 곳 사이의 구간의 길이가 됩니다.

이전의 예제들과 비스한 설정을 사용하면 위와 같이 결정됩니다. H0는 1로 그냥 두겠습니다.

그러면 위 그림처럼 아~주 깔끔한 Band Stop 필터가 만들어 졌네요^^ 아무튼 위 그림의 Band Pass Filter는 

입니다. 그리고, 위 그림을 얻기 위한 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
H0 = 1

num = np.array([H0,0,w0_peak**2])
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()

입니다.^^

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

앞선 예제에서도 줄기차게 이야기했지만, s-domain의 필터를 디지털로 뭐 어디 MCU라도 넣어볼려면 digital 필터로 변환되어야 하는데 그 때 사용되는 z-domin으로 변환하는 식은

입니다.

결국 대입해서 정리한 후에 위와 같은 형태로 정리하면 각 계수는...

와 같이 정리가 됩니다.

그걸~~ 위와 같이 정리한다고 보면...

이 될 것입니다.

그러면.. z-domain에서는 위와 같이 그려질 것이고, 다시 Bode Plot을 그려보면,

이 됩니다. 잘 변환되었음을 알 수 있죠^^ 위 그림을 그린 Python 코드는 아래와 같습니다.

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

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

numz1 = np.array([H0*4/Ts**2 + H0*w0_peak**2,
                  -2*H0*4/Ts**2 + 2*w0_peak**2,
                  H0*4/Ts**2 + H0*w0_peak**2])
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(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()

디지털 필터 구현하기

이제 다시 본래의 목적(글 처음에는 말하지 않았지만, 제 목적인^^)대로 디지털로 구현해보도록 하겠습니다.

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

위 코드를 이용해서.. 시험 신호를 만들도록 하겠습니다.

사실 이 신호는 글 처음에 언급한 예전 글에서도 사용한 형태인데요. 이번에느 Band Pass Filter와 동일한 시험신호를 사용하겠습니다. 이 신호는 코드로 만들었기 때문에 정확하지만, 혹시하는 마음에 FFT를 수행해보면...

이렇습니다.^^.

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


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

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

위 코드에 나타나듯이 DirectIIForm으로 직접 필터링을 해보면... 그 결과는

위와 같이 딱 500Hz 중심으로 얼마의 구간이 빠진 형태가 됩니다. 모양이 Bode 선도와 흡사하죠^^ 시간영역에서의 결과는

위 그림과 같이 원 신호에 대해 특정 주파수만 빠진 형태로 나타나게 됩니다.^^ 아래에는 참조로 본 글에 사용된 전체 코드입니다.

# -*- 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
H0 = 1

num = np.array([H0,0,w0_peak**2])
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([H0*4/Ts**2 + H0*w0_peak**2,
                  -2*H0*4/Ts**2 + 2*w0_peak**2,
                  H0*4/Ts**2 + H0*w0_peak**2])
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(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()