PythonのSciPyでバンドパスフィルタをかける!

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

ローパス、ハイパスフィルタの応用となるバンドパスフィルタは信号のノイズ除去以外にも、特徴量の抽出に役立ちます。フィルタには実際難しい理論がありますが、PythonのSciPyなら簡単に使うことが可能です。

こんにちは。wat(@watlablog)です。
これまで学んだフィルタ処理の応用編として、このページではバンドパスフィルタを紹介します

バンドパスフィルタについて

バンドパスフィルタとは?

WATLABブログでは前回までに、高周波帯域ノイズを除去したい時に使うローパスフィルタ、低周波帯域の信号をキャンセルして高周波帯域成分を抽出するハイパスフィルタをPythonでコーディングして来ました。

今回はその両方の特性を持つバンドパスフィルタband-pass filter)を紹介します。

バンドパスフィルタとは、ある周波数幅(バンド)のみの信号を通過(パス)させるフィルタのことで、ローパスフィルタとハイパスフィルタの組み合わせによって構成されます。

バンドパスフィルタは光学の分野でよく使われており、特定波長の光のみを抽出する目的で利用されます。

しかし、「特定の周波数帯域だけを通す」という機能は、工学的にも幅広い用途があるためあらゆる場面で活用されています。

バンドパスフィルタで重要な因子と用語(図解)

バンドパスフィルタはある周波数幅の帯域だけを通過させるため、以下の図のように1つの通過域がありますが、阻止域遷移域が2つあるのが特徴です。

そのため、通過域端と阻止域端が2つずつあることになり、これまでのローパスフィルタやハイパスフィルタと比べるとプログラムする上で設定項目が増えることになります。

今回もバターワースフィルタを例にコーディングしていきます。

PythonのSciPyバンドパスフィルタ

Advertisements

ソースコード

メインファイル

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

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

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

中身ですが、ローパスやハイパスの時とパラメータの種類は同じです。しかし通過域端周波数fpと阻止域端周波数fsは上図の理由で2つあるので、以下のコードのようにベクトル形式で入力する必要があります。そしてこのベクトルはNumPy形式で記述する必要があるため、np.arrayを使っています。

関数ファイル

関数ファイル(filter_function.py)もローパスやハイパスとほぼ同じです。signal.butterでbandを指定する所が変化点です。

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

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

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

実行結果

バンドパスフィルタをガウシアンノイズに対してかけた結果です。
下図上段が時間波形、下段が上段の波形をFFTした後の周波数波形です。それぞれフィルタ前(original)とフィルタ後(filtered)を比較しています。

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

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

これまでと同じように、ノイズ波形だけではつまらないので、サイン波を二つ合成させた波形を作って効果を確認していきます。

上のメインコード(main_filter.py)の波形生成部分とフィルタ周波数を以下のように変更してみます。

元波形は1000[Hz]の信号が200[Hz]でうねり成分を持ち、かつ全体的に高周波ノイズが含まれていましたが、バンドパスフィルタをかけたことで1000[Hz]の信号を比較的綺麗に抽出することができました。

まとめ

Advertisements

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

バンドパスフィルタはローパスフィルタやハイパスフィルタと比べ阻止域が2つあるため、周波数の設定はベクトル形式で用意する必要がありました。

SciPyを使えばバンドパスフィルタを容易に設計できることを学びました。

バンドパスは特定の信号を抽出する用途で使うため、何か特徴量を抽出する分野では万能に使えそうだね!フィルタ処理の基礎もあと少し!

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

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

SNSでもご購読できます。

コメント

  1. TKI より:

    質問失礼いたします。Python初心者です。
    ソースコードを実行すると”name ‘signal’ is not defined”とエラーが吐き出されてしまいます。
    signalは別で定義が必要ということでしょうか?
    又は、scipyからimportすれば使用出来るのでしょうか?

    おかしな質問でしたら申し訳ありません。
    ご回答いただけると嬉しいです。
    宜しくお願い致します。

    1. wat より:

      ご訪問ありがとうございます。
      signalはfrom scipy import signalとインポートします。
      こちらのページでは省略していましたので、わかりにくかったかも知れません。
      新しく「参考:フーリエ変換まで含めた全コード」という項目を書いてみました。
      こちらのコードを貼り付けてそのまま動きますでしょうか?
      よろしくお願い致します。

  2. hiroko より:

    初心者ですが、ブログを楽しみにしています。ありがとうございます。

    任意のcsvファイルのデータに対してバンドパスフィルタ(0.5Hz~30Hz)を使いたいと思っています。
    この例はガウシアンノイズに対してバンドパスフィルタを用いたと思い、の部分を書き換えてみました。
    用いたCSVファイルは、「ただPythonでcsvから離散フーリエ変換をするだけのコード」に載っていたcsvファイルを参考にしています。

    samplerate = 1024
    x = np.arange(0, 512) / samplerate
    df_master = pd.read_csv(‘signals.csv’, index_col=0)
    data = df_master.iloc[:,[2]]

    fp = np.array([1, 30])
    fs = np.array([0.5, 60])
    gpass = 3
    gstop = 40


    すると、以下のエラーが出たのでてしまいましたがよく理解できません。gpass/gstopの値を変更すると、”45″の値は変化します。
    The length of the input vector x must be greater than padlen, which is 45.

    設定方法について、教えて頂けますか。(参考になるページでも嬉しいです。)

    1. wat より:

      ご質問ありがとうございます。
      こちらでエラーが再現せず回答が遅れてしまって申し訳ございません。
      (解決されていましたら無視して下さい…)

      そのエラーはベクトルの長さに関するもののようで、データ長に不整合があるようです。

      xは新たに定義されていますが、csvはそのまま使っているのでしょうか?

      ただPythonでcsvから離散フーリエ変換をするだけのコード」を参考にして頂いているとのことで、
      僕の方でこの記事で紹介しているコードをベースにバンドパスフィルタを追加して動作させてみたところ、特にエラーなく動作しました。

      コード変更部分は以下です。
      先頭にfrom scipy import signalを書き、バンドパスの関数を置いておいた上で、フーリエ変換前にフィルタ処理を追加しただけです。

      # 列方向に順次フーリエ変換(DFT)をするコード
      for i in range(len(df.T)-1):
      print(i)
      data = df.T.iloc[i+1] # フーリエ変換するデータ列を抽出

      # フィルタ処理
      fp = np.array([1, 30]) # 通過域端周波数[Hz]※ベクトル
      fs = np.array([0.5, 60]) # 阻止域端周波数[Hz]※ベクトル
      gpass = 3 # 通過域端最大損失[dB]
      gstop = 40 # 阻止域端最小損失[dB]
      data = bandpass(x=data, samplerate=1/dt,
      fp=fp, fs=fs,
      gpass=gpass, gstop=gstop)
      data = pd.Series(data)

      spectrum, amp, phase, freq = calc_fft(data.values, 1/dt) # フーリエ変換をする関数を実行

      まずは上記記事のコードをベースに動かしてみて、必要であればカスタマイズをお試し頂ければと思います。

      フィルタに関する質問は多いのと、どうやら他の質問サイト様でも需要がある様子がわかりましたので、別途まとめ記事を書こうと思います。
      近日公開しますので、もしよろしければご確認下さい。

      以上、よろしくお願い致します。

      1. hiroko より:

        ありがとうございます。出来ました!
        先日更新されていたブログも拝見しました。ありがとうございます。

        ブログで紹介されていたPyQも始めてみました(^^ゞ。

        1. wat より:

          ご連絡ありがとうございました。
          問題ないようで安心しました!
          PyQ色々なコンテンツあって面白いですよね。
          是非今後もブログでわからない事やリクエストがありましたらご連絡ください!

TKI へ返信する コメントをキャンセル

*