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

　今回は、Pythonを用いて信号データに対してローパスフィルタやハイパスフィルタを適用する方法について紹介する。

1. フィルタ関数の定義

フィルタを適用する方法としてはScipyの関数を利用する。Scipyの信号処理関連のモジュールをロードするにはimport scipy import signalとする。

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

フィルタを適用する流れとしては、フィルタの次数と正規化周波数N,Wnを計算し、そのN,Wnからフィルタ伝達関数の分子、分母b,aを計算し最後にb,aを信号xに適用するという形になる。

今回はバタワースフィルタを設計、適用する例を紹介する。まずナイキスト周波数で正規化された周波数を入力パラメータとして使用するため、信号のサンプリングレート(サンプル数/秒）を用いてナイキスト種波数fnを算出する。次にローパスフィルタの通過周波数wp：指定した周波数以下は減衰を受けない、と阻止周波数ws：指定した周波数以上は減衰を受けるをナイキスト周波数で正規化された値を指定する。またwpでの減衰量gpassをdBで指定するが通常は3dBにすることが多い。あとgstopは阻止帯域での減衰量で今回は40dBとする。 4つの設計パラメータから次数N、正規化周波数Wnを求めるには signal.buttord(wp, ws, gpass, gstop)とする。得られたN,Wnからsignal.butter(N, Wn, "low")とすることでフィルタ伝達関数の係数b,aを取得し、最後にsignal.filtfilt(b, a, x)とすることで信号xにフィルタ伝達関数を適用した信号を取得できる。以上の処理をフィルタを適用する関数として定義する。

def lowpass(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, "low")            #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x)                  #信号に対してフィルタをかける
    return y

上記コードや細かい説明は以下を参照

watlab-blog.com

2. フィルタを適用するデータの生成

　フィルタを適用するテスト用の信号として正規分布乱数をノイズとするが、途中でベースレベルがステップ状に変わるものを用意する。

samplerate = 25600
x = np.arange(0, 12800) / samplerate                    # 波形生成のための時間軸の作成
data = np.random.normal(loc=0, scale=1, size=len(x))    # ガウシアンノイズを生成
data2=np.random.normal(loc=0, scale=1, size=len(x))+3
data3=np.concatenate([data[0:int(len(x)/2)],data2[0:int(len(x)/2)]])
plt.plot(x,data3)

上記コードで生成されるテスト信号は以下のようになる。

3. ローパスフィルタの生成と適用

　今回は通過周波数を300Hz,阻止周波数を600Hzとする。理想的には通過周波数と阻止周波数はなるべく近づけた方がよさそうだが、あまり両者を近づけると計算はできるが適切なフィルタ伝達関数の係数が見つからずフィルタがかからない。同様にgpass=0やgstopに大きな数値を入れても係数が見つからない場合がある。上で定義した関数lowpass()をフィルタ設計用の4つのパラメータ,fp,fs,gpass,gstopとテスト信号data3を用いてローパスフィルタ処理を実行する

fp = 300 # 通過域端周波数[Hz]
fs = 600 # 阻止域端周波数[Hz]
gpass = 3 # 通過域端最大損失[dB]
gstop = 40 # 阻止域端最小損失[dB]
 
# ローパスをする関数を実行
data_lofilt = lowpass(data3, samplerate, fp, fs, gpass, gstop)
plt.plot(x,data3)
plt.plot(x,data_lofilt,'r')
plt.legend(['Raw','Lowpass'])

フィルタ適用前のデータをRaw, フィルタ適用後をLowpassとして以下のように、ノイズが除去されており効果が確認できる。

3. ハイパスフィルタの生成と適用

　今回作成した関数Lowpassの一部をアレンジしてハイパスフィルタを生成・適用することもできる。途中まではすべてローパスと同じでフィルタ伝達関数の係数を算出するsignal.butter(N, Wn, "high")とするだけでハイパスフィルタを生成できる。

def hipass(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, "high")   #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x)    #信号に対してフィルタをかける
    return y

通過周波数と阻止周波数の数値はローパスの逆に設定して、ローパスフィルタと同じ手順でテスト信号にハイパスフィルタを適用する。

fp=600 
fs=300
data_hifilt = highpass(data3, samplerate, fp, fs, gpass, gstop)
plt.plot(x,data3)
plt.plot(x,data_hifilt,'r')
plt.legend(['Raw','Highpass'])

以下、ハイパスフィルタなのでノイズはそのままで、ステップ変化する部分は除去されていることが確認できる。