Translate

scipyで手っ取り早くバダーワースフィルタを使えるようにする【butterworth filter】

最近プログラミングの環境をMATLABから徐々にpythonに移行することができてきて、MATLABと同じくらいに使えるようになってきております。

やっぱり、pythonは無料でいろんなパッケージを使えるのがかなりのメリットですね。あとバージョンが古いMATLABだと最新のコードがうまく走らないことが多くて、それも解消できたため、非常にストレスフリーでプログラミングができております。

そんな中で、最近頻回に使っているscipyの中のバダーワースフィルタ(butterworth filter)について備忘録として残しておこうと思います。


バダーワースフィルタは一度関数さえ書いておけば簡単。MATLABっぽく使える。



波形解析ではほぼ必須ともいえるバターワースフィルター(butterworth filter )ですが、MATLABだとあらかじめ関数がパッケージ化されているため、比較的簡単に使うことができます。

MATLABからのpython移民であった甘味としては、これを何とかまずはうまく実装できないかなと思いいろいろと調べてみました。

なんかいろいろ見ていると、いろんな設定とか(通過域、阻止域など)とっつきにくいフィルタ設計などがあり、敷居が高そうに見える。

(;甘)。。。MATLABの時には、カットしたい周波数領域と次数の設定くらいしか入れていなかったんだけどな。なんかMATLABっぽく簡単にかけるscipyのpythonコードはないのだろうか。


大まかに説明(私の拙い理解の範囲です。間違っていたらすみません)すると、はじめにフィルタの設計を行う関数を定義して、そのあとに設計された値を用いてフィルタを実際に走らせるという二つの工程に分かれた関数になっています。

上記のサイトを参考にローパスフィルタを書いてみたらこんな感じになります。4次のローパスフィルタで、引数に、フィルタをかけたいデータ(x)カットオフ周波数(cutoff)と、サンプリング周波数(fs→例えば100Hzとか)になります。
from scipy.signal import butter, filtfilt,lfilter
from scipy import signal
def butter_lowpass(cutoff, fs, order=4):
    nyq = 0.5 * fs
    cutoff = cutoff / nyq
    b, a = butter(order, cutoff, btype='low', analog=False, fs=None)
    return b, a
def butter_lowpass_filter(x, cutoff, fs, order=4):
    b, a = butter_lowpass(cutoff, fs, order=order)
    return filtfilt(b, a, x,axis=0)
    
この関数をスクリプトの最初の行に入れておいて必要な時に使うといった感じですね。ハイパスフィルタの場合には、上記の関数のlowの部分をhighに書き換えるだけでほぼ行けるはず。

実際に使う時にはこんな感じ。

data = butter_lowpass_filter(data,10,100)
    
簡単なプロセスの説明をすると、100Hzで計測されたデータ(data)に対して、10Hzのローパスフィルタをかける処理を行うという流れになります。先ほどの関数では二つありましたが、二つ目の関数を呼び出すとのその関数内で、一つ目に定義した関数(butter_lowpass)が使われて、フィルタ設計がされて、そこから戻される戻り値b,aを使ってフィルタ処理(butter_lowpass_filter)が行われる感じです。

たぶん本当は通過域、阻止域とかを厳密に使ってやるのだろうけど、さっくりとフィルタ処理を行いたいだけであれば、こんな感じでいけると思います。ちなみに、投入するデータはしっかりと一列毎に切り出さないとエラーを吐きます。

複数列にフィルタ処理をかけたい場合にはループ回路をかませるのが手っ取り速いでしょう(あんまりループ回路が多くなると、処理が重たくなりますが(笑))

 for k in range (0,data.shape[1]):
        data[:,k] = butter_lowpass_filter(data[:,k],10,100)
    
ループ回路使って複数列のフィルタ処理をするのだとこんな感じですかね。上記の二つの関数(butter_lowpass、butter_lowpass_filter)をあらかじめスクリプトに書いておきましょう。ループ回路の説明としては、m行k列の二次元データ配列dataを一列目からk列までローパスフィルタ処理を行うといった内容になります。基本numpyを使った二次元配列の場合にはこのコードで動くはす。。。(エラー出ていたらごめんなさい)