信号にノイズが含まれる場合、フィルタ処理をして波形を整形することがあります。フィルタ処理は難しい理論がありますが、PythonのSciPyであれば使うだけなら簡単にできます。
こんにちは。wat(@watlablog)です。
今回はPythonにおけるローパスフィルタの処理について学びます！
目次（項目クリックでジャンプできます）
ローパスフィルタについて
ローパスフィルタとは？
フィルタとは、日本語で言うと濾過器のことをさします。信号処理の世界では、ある信号から必要な成分のみを抽出する時に使います。
ローパスフィルタ(Low-pass filter)とは、信号の低周波数帯域の成分のみを通過（パス）させ、高周波帯域の成分は阻止（カット）するフィルタのことを指します。
発祥は電気分野とのことですが、電気の世界ではとにかく信号にノイズが乗るような対象を取り扱うので、特に高周波のノイズ成分を除去できるローパスフィルタが多用されています。
ローパスフィルタで重要な因子と用語（図解）
ローパスフィルタは低周波数帯域の信号を通し、高周波帯域の信号を阻止するので、周波数軸で考えると下図のように3つの領域（「通過域」「遷移域」「阻止域」）に分けられます。
フィルタは通過域の信号はほとんど減衰させないものを指しますが、阻止域で直角に減衰するのではなく、遷移域と呼ばれる徐々に減衰量が変わっていく範囲を持ちます。
さらに、通過域の端部である通過域端周波数でもある程度の減衰量を持ちます。
各周波数でどのような減衰特性を持たせるか、どこまで許容させるのかといった要素を決めて行くのが、フィルタ設計の主目的※です。
※フィルタ処理も元の波形を変化させてしまう信号処理技術であることには変わりがないので、減衰量の設計以外にも信号の位相変化を許容するかどうか、安定性はあるかどうか…といった要素もフィルタ設計にはありますが、今回はそこまで細かい話には踏み込まず、工学の現場でよく使われているバターワースフィルタを扱います。
PythonのSciPyローパスフィルタ
ソースコード
メインファイル
ソースコードを以下に示します。フィルタ処理を実際に行う部分はdef文を使って関数化しました。関数の使い方については「Pythonの関数def文の使い方！引数や別ファイル式も解説」の記事に詳細を示しています。
まずはメイン(main_filter.py)のコードです。
サンプルとしてガウシアンノイズを生成していますが、こちらについては「Pythonでガウス分布を持つノイズの作り方と調整方法」の記事を参照下さい。
中身ですが、上の図で説明した通過域端周波数fp[Hz]と阻止域端周波数fs[Hz]を設定します。またそれぞれの周波数における減衰量もgpassとgstopで設定します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
import filter_function import numpy as np from matplotlib import pyplot as plt samplerate = 25600 #波形のサンプリングレート x = np.arange(0, 12800) / samplerate #波形生成のための時間軸の作成 data = np.random.normal(loc=0, scale=1, size=len(x)) #ガウシアンノイズを生成 fp = 3000 #通過域端周波数[Hz] fs = 6000 #阻止域端周波数[Hz] gpass = 3 #通過域端最大損失[dB] gstop = 40 #阻止域端最小損失[dB] data_filt = filter_function.lowpass(data, samplerate, fp, fs, gpass, gstop) |
関数ファイル
以下の関数を使ってローパスフィルタを通した波形を取得します。ここで、設定した周波数をナイキスト周波数で除していますが、これはデジタルフィルタを使う時に周波数を正規化する必要があるからです。
以下に引用する公式ページに使い方が書いてあります（日本語検索では意味まで説明していた記事は見当たらず…）。
For digital filters, these are normalized from 0 to 1, where 1 is the Nyquist frequency.
SciPy.org:scipy.signal.buttord
僕はここの正規化の使い方、考え方のエビデンスをとるのに時間がかかってしまいました。初心者の壁かも知れません！
1がナイキスト周波数である必要があるので、サンプリングレートの半分であるナイキスト周波数で割っています。公式ページでは周波数が[Hz]ではありませんが、比をとっているだけなので、設定した周波数とナイキスト周波数がともに[Hz]で書いてあれば問題ありませんね。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
import numpy as np from scipy import signal #バターワースフィルタ（ローパス） 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 #フィルタ後の信号を返す |
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はベクトル形式になっています。
実行結果
下図上段が時間波形、下段が上段の波形をFFTした後の周波数波形です。それぞれフィルタ前（original）とフィルタ後（filtered）を比較しています。
ここではFFTやデシベル変換も使っていますが、それぞれの詳細は以下の記事にまとめていますので、是非ご覧ください。
PythonでFFT！SciPyのFFTまとめ
Pythonで音圧のデシベル(dB)変換式と逆変換式
元の時間波形はガウシアンノイズなので周波数全体でフラットな特性を持っていますが、フィルタ後は設定した周波数(fp=3000[Hz]）から綺麗に減衰していく結果を得ました。
ノイズ波形だけではつまらないので、サイン波のノイズ除去を行ってみましょう。上のメインコード（main_filter.py）の波形生成部分とフィルタ周波数を以下のように変更してみます。
|
1 2 3 4 5 6 7 8 |
samplerate = 25600 x = np.arange(0, 12800) / samplerate # 波形生成のための間軸の作成 data = np.sin(2.0 * np.pi * 50 * x) + 0.3 * np.random.normal(loc=0, scale=1, size=len(x)) #サイン波にガウシアンノイズを重畳 fp = 200 #通過域端周波数[Hz] fs = 400 #阻止域端周波数[Hz] gpass = 3 #通過域端最大損失[dB] gstop = 40 #阻止域端最小損失[dB] |
すると元の時間波形はノイジーなサイン波、フィルタ後はノイズが軽減している結果を得ます。このページでは極端な例を紹介していますが、これがフィルタ処理の効果です。
ソースコードの周波数や減衰量を変えると、減衰が始まる位置、遷移域の傾きを自在に変化させることができますので、是非トライしてみて下さい。
まとめ
フィルタ処理の一種であるローパスフィルタをPythonで使えるようになりました。
ローパスフィルタは周波数軸で通過域、遷移域、阻止域における周波数と減衰量を決める必要があります。
フィルタをかける時に使用する伝達関数の入力値はSciPyの関数を使うことで簡単に求めることができました。
理解するのに時間はかかったけど、ようやくフィルタ処理ができるようになってきたぞ！ここから各種フィルタに発展させていこう！
また、英語の公式ページを見ると理解が深まることがわかったので、これからは積極的に海外Webページも活用していこう！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！
コメント