PythonのSciPyでローパスフィルタをかける!

  • このエントリーをはてなブックマークに追加

信号にノイズが含まれる場合、フィルタ処理をして波形を整形することがあります。フィルタ処理は難しい理論がありますが、PythonのSciPyであれば使うだけなら簡単にできます。

こんにちは。wat(@watlablog)です。
今回はPythonにおけるローパスフィルタの処理について学びます

ローパスフィルタについて

ローパスフィルタとは?

フィルタとは、日本語で言うと濾過器のことをさします。信号処理の世界では、ある信号から必要な成分のみを抽出する時に使います。

Advertisements

ローパスフィルタ(Low-pass filter)とは、信号の低周波数帯域の成分のみを通過(パス)させ、高周波帯域の成分は阻止(カット)するフィルタのことを指します。

発祥は電気分野とのことですが、電気の世界ではとにかく信号にノイズが乗るような対象を取り扱うので、特に高周波のノイズ成分を除去できるローパスフィルタが多用されています。

ローパスフィルタで重要な因子と用語(図解)

ローパスフィルタは低周波数帯域の信号を通し、高周波帯域の信号を阻止するので、周波数軸で考えると下図のように3つの領域(「通過域」「遷移域」「阻止域」)に分けられます。

ローパスフィルタの図

フィルタは通過域の信号はほとんど減衰させないものを指しますが、阻止域で直角に減衰するのではなく、遷移域と呼ばれる徐々に減衰量が変わっていく範囲を持ちます。

さらに、通過域の端部である通過域端周波数でもある程度の減衰量を持ちます。

各周波数でどのような減衰特性を持たせるか、どこまで許容させるのかといった要素を決めて行くのが、フィルタ設計の主目的※です。

※フィルタ処理も元の波形を変化させてしまう信号処理技術であることには変わりがないので、減衰量の設計以外にも信号の位相変化を許容するかどうか、安定性はあるかどうか…といった要素もフィルタ設計にはありますが、今回はそこまで細かい話には踏み込まず、工学の現場でよく使われているバターワースフィルタを扱います。

PythonのSciPyローパスフィルタ

ソースコード

メインファイル

ソースコードを以下に示します。フィルタ処理を実際に行う部分はdef文を使って関数化しました。関数の使い方については「Pythonの関数def文の使い方!引数や別ファイル式も解説」の記事に詳細を示しています。

まずはメイン(main_filter.py)のコードです。

サンプルとしてガウシアンノイズを生成していますが、こちらについては「Pythonでガウス分布を持つノイズの作り方と調整方法」の記事を参照下さい。

中身ですが、上の図で説明した通過域端周波数fp[Hz]と阻止域端周波数fs[Hz]を設定します。またそれぞれの周波数における減衰量もgpassgstopで設定します。

関数ファイル

以下の関数を使ってローパスフィルタを通した波形を取得します。ここで、設定した周波数をナイキスト周波数で除していますが、これはデジタルフィルタを使う時に周波数を正規化する必要があるからです。

以下に引用する公式ページに使い方が書いてあります(日本語検索では意味まで説明していた記事は見当たらず…)。

For digital filters, these are normalized from 0 to 1, where 1 is the Nyquist frequency.

SciPy.org:scipy.signal.buttord

僕はここの正規化の使い方、考え方のエビデンスをとるのに時間がかかってしまいました。初心者の壁かも知れません!

1がナイキスト周波数である必要があるので、サンプリングレートの半分であるナイキスト周波数で割っています。公式ページでは周波数が[Hz]ではありませんが、比をとっているだけなので、設定した周波数とナイキスト周波数がともに[Hz]で書いてあれば問題ありませんね。

Wnも同様です。正規化する必要があります。そしてこのWnには-3[dB]ポイントという特別な意味合いを持ちます。バターワースフィルタのWn周波数では減衰量が\(1/\sqrt{2}\)(半値、-3[dB])になります。

For a Butterworth filter, this is the point at which the gain drops to 1/sqrt(2) that of the passband (the “-3 dB point”).

SciPy.org:scipy.signal.butter

signal.butterで計算しているbとaは伝達関数の分母分子です。そもそもローパスフィルタを始めとした各種フィルタはラプラス変換を利用して伝達関数の形式で表すことができます。伝達関数の分母と分子はそれぞれ多項式になっているため、このb, aはベクトル形式になっています。

参考:フーリエ変換まで含めた全コード

先ほどまではフィルタ処理の部分のみを解説していましたが、ファイルを2つ作ったり面倒だと思いますので、コピペしてすぐ動かせる全コードを以下に示します。

このコードはサンプル波形を生成して、フーリエ変換とデシベル変換を行っており、フィルタ前後の波形を比較する事が可能です。是非動かして遊んでみて下さい。

実行結果

下図上段が時間波形、下段が上段の波形をFFTした後の周波数波形です。それぞれフィルタ前(original)とフィルタ後(filtered)を比較しています。

ここではFFTやデシベル変換も使っていますが、それぞれの詳細は以下の記事にまとめていますので、是非ご覧ください。

PythonでFFT!SciPyのFFTまとめ
Pythonで音圧のデシベル(dB)変換式と逆変換式

実行結果の図

元の時間波形はガウシアンノイズなので周波数全体でフラットな特性を持っていますが、フィルタ後は設定した周波数(fp=3000[Hz])から綺麗に減衰していく結果を得ました。

ノイズ波形だけではつまらないので、サイン波のノイズ除去を行ってみましょう。上のメインコード(main_filter.py)の波形生成部分とフィルタ周波数を以下のように変更してみます。

すると元の時間波形はノイジーなサイン波、フィルタ後はノイズが軽減している結果を得ます。このページでは極端な例を紹介していますが、これがフィルタ処理の効果です。

ソースコードの周波数や減衰量を変えると、減衰が始まる位置、遷移域の傾きを自在に変化させることができますので、是非トライしてみて下さい。

まとめ

フィルタ処理の一種であるローパスフィルタをPythonで使えるようになりました。

ローパスフィルタは周波数軸で通過域、遷移域、阻止域における周波数と減衰量を決める必要があります。

フィルタをかける時に使用する伝達関数の入力値はSciPyの関数を使うことで簡単に求めることができました。

理解するのに時間はかかったけど、ようやくフィルタ処理ができるようになってきたぞ!ここから各種フィルタに発展させていこう!
また、英語の公式ページを見ると理解が深まることがわかったので、これからは積極的に海外Webページも活用していこう!

Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。

コメント

  1. Tuk より:

    いつも参考にさせていただいています。
    時系列データの全体に載っているノイズを除去したく、ローパスフィルターの上記コードを利用しているのですが、gstopの値の設定方法がよくわかりません。
    ①まずFFTによって周波数データに変換し、縦軸の振幅値をdB変換するのですが、0Hzにおける振幅値をA0として、dB値は20*log10(A/A0)で計算は合っていますでしょうか?
    ②gstopの値はどのようにして設定するのが妥当でしょうか?
    バタワースフィルターでは減衰傾度がある以上、一定周波数以上のデータを完全にカットする、というのは難しいのですが、何とかきれいにノイズ除去できる方法がないかと模索中です。。
    よろしくお願いします。

    1. wat より:

      いつもありがとうございます。
      gpassとgstopの使い方は、まずは正式に1次情報である以下の公式サイトをご参照ください。
      https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.buttord.html

      この記事の図の「通過域端減衰量」がgpassで、「阻止域端減衰量」がgstopとなります。
      0[Hz]における振幅値は関係なく、単純に「信号をどれだけ減衰させたいか」で設定します。

      ただ、仰るとおりの理由で完全にノイズをカットすることは難しいです。
      完全に周波数でノイズをカットする1手法として、フーリエ変換して関連帯域を0にしてから(複素数のフーリエスペクトル、負の周波数帯域も0にする)逆フーリエ変換することが考えられます。
      但し、フーリエ変換を使う方法はある程度の区間を平均化するため、時間軸で変動のある信号に対しては万能ではありません。

      フィルタを綺麗にかけるというのは「フィルタ設計」と呼ばれるほど奥が深い分野ですので、より専門的にはこのキーワードで勉強する必要があると思います。
      ざっとGoogle検索するとやはりこの通過域端、阻止域端の工夫が研究されているようです。
      (僕はこの分野のプロではないので、この領域はまだ未熟。。。あまり深い回答はできず申し訳ないです。このページに疑問を投げかけることで、もしかしたらプロから回答があるかも…?)

      是非色々調べてください。ただ不明点があればこちらも勉強できるので、気軽にお声がけください!

コメントを残す

*