Pythonによるデータ処理5 ～フィルタ処理2

　今回は、Pythonを用いて信号データに対してローパスフィルタ、ハイパスフィルタに続いてバンドパス/ストップフィルタを適用する方法について紹介する。

atatat.hatenablog.com

1. バンドパスフィルタ

バンドパスフィルタについても前回のローパス・ハイパスと同様でフィルタ伝達関数の係数を算出するsignal.butter(N, Wn, "band")と設定する。

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

#バターワースフィルタ（バンドパス）
def bandpass(x, samplerate, fp, fs, gpass, gstop):
    fn = samplerate / 2 #ナイキスト周波数
    wp = fp / fn  #ナイキスト周波数で通過域端周波数を正規化
    ws = fs / fn  #ナイキスト周波数で阻止域端周波数を正規化
    N, Wn = signal.buttord(wp, ws, gpass, gstop)  #オーダーとバターワースの正規化周波数を計算
    b, a = signal.butter(N, Wn, "band") #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける
    return y

バンドパスフィルタは、ある周波数からある周波数まで通過させる特性上、通過周波数は値を2つ有する。同様に阻止周波数の値も2つ有することになる。 2つの値を指定するためには、通過域端周波数と阻止域端周波数はarrayを用いてsignal.buttordに入力する。今回の例は通過周波数を1000Hzから3000Hz, 阻止周波数は500Hz以下と6000Hz以上に設定するのでfp = np.array([1000,3000]),fs = np.array([500,6000])とする。

テスト用の信号はガウシアンノイズに加えて、阻止帯域下側の信号として200Hzの正弦波、阻止帯域上側の信号として8000Hzの正弦波、あと通過帯域の信号として2000Hzの正弦波を加算した信号を用意する。通過帯域と阻止帯域の損失は前回と同じくgpass=3,gstop=40とする。

samplerate = 25600
fl=200 # 阻止帯域下側の信号
fh=8000 # 阻止帯域下側の信号
fb=2000 # 通過帯域の信号

x = np.arange(0, 12800) / samplerate  # 波形生成のための時間軸の作成
data = np.random.normal(loc=0, scale=1, size=len(x))    # ガウシアンノイズを生成
data2= np.sin(2*np.pi*fl*x)+np.sin(2*np.pi*fh*x)+np.sin(2*np.pi*fb*x) # 阻止帯域と通過帯域の信号生成
data3=data+data2
 
fp = np.array([1000,3000])   # 通過域端周波数[Hz]
fs = np.array([500,6000]) # 阻止域端周波数[Hz]
gpass = 3 # 通過域端最大損失[dB]
gstop = 40# 阻止域端最小損失[dB]
 
# バンドパスをする関数を実行
data_bandfilt = bandpass(data3, samplerate, fp, fs, gpass, gstop)
plt.plot(x,data3)
plt.plot(x,data_bandfilt,'r')
plt.legend(['Raw','Bandpass'])

実行すると上ののように全体的には信号幅が狭まったような波形になるが、詳細がわかりづらいので一部分を拡大して再プロットする。

plt.plot(x,data3)
plt.plot(x,data_bandfilt,'r')
plt.xlim([0,0.01])
plt.legend(['Raw','Bandpass'])

高周波のノイズ成分と、200Hzのうねりは除去できているが、通過帯に設定した2000Hzの信号は通過帯域中のガウシアンノイズにより歪んではいるが残っている。

バンドパスフィルタの効果を明確に確認するため、ガウシアンノイズなしにしたテスト用の信号を用いる。

data_bandfilt2 = bandpass(data2, samplerate, fp, fs, gpass, gstop)
plt.plot(x,data2)
plt.plot(x,data_bandfilt2,'r')
plt.xlim([0,0.01])
plt.legend(['Raw','Bandpass'])

以下のように低域の200Hzと高域の8000Hzの正弦波が除去され、2000Hzの信号だけが残っていることがはっきりと確認できる。

"バンドパスフィルタ適用(ノイズなし)" — バンドパスフィルタ適用(ノイズなし)

2. バンドストップフィルタ

バンドパスと同じ要領で、signal.butter(N, Wn, "bandstop")とすればバンドストップフィルタを適用することができる。

#バターワースフィルタ（バンドストップ）
def bandstop(x, samplerate, fp, fs, gpass, gstop):
    fn = samplerate / 2 #ナイキスト周波数
    wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化
    ws = fs / fn   #ナイキスト周波数で阻止域端周波数を正規化
    N, Wn = signal.buttord(wp, ws, gpass, gstop)  #オーダーとバターワースの正規化周波数を計算
    b, a = signal.butter(N, Wn, "bandstop")            #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x)                  #信号に対してフィルタをかける
    return y

テスト用の信号は上記のガウシアンノイズなしの信号を用いる。通過帯域と阻止帯域の関係が逆になるだけであとはバンドパスと同じになる。

fp = np.array([500, 6000]) #通過域端周波数[Hz]
fs = np.array([1000, 3000]) #阻止域端周波数[Hz]
gpass = 3 #通過域端最大損失[dB]
gstop = 40 #阻止域端最小損失[dB]
data_bandstopfilt = bandstop(data2, samplerate, fp, fs, gpass, gstop)
plt.plot(x,data2)
plt.plot(x,data_bandstopfilt,'r')
plt.xlim([0,0.01])
plt.legend(['Raw','Bandstop'])