デジタルフィルタに関しての知識が何もない状態で、ローパス(ハイパス)フィルタを作ることになってしまって困ったのでメモです。
今回はscipyのバターワースフィルター(scipy.signal.butter)を使用しています。
フィルタの設計
まずはどういったフィルタにしたいかを考える必要があります。
ここで考えなくてはいけないのはフィルタの次数とカットオフ周波数です。
フィルタの次数
フィルタの次数は適用するフィルタがどのくらい急に変化するかを表しています。
このフィルタ次数が大きいほど急なフィルタになります。
大きいほど理想的なフィルタ形状となりますが、次数を大きくしすぎるとフィルタの周波数特性が悪くなるらしいので程よい値を使う必要があるみたいです。
以下に、次数を変えたときのハイパスフィルタのフィルタ形状を示します。
次数(order)を増やすとフィルタの傾きが大きくなっていることが確認できます。
カットオフ周波数
カットオフ周波数はどの周波数からフィルタを有効にするかを表しています。
(正確には元波形と比較して-3dB減衰する周波数のことですが。)
フィルタ次数の図でもわかる通り、カットオフ周波数できっかり0/1になるようなフィルタを作成することはできません。
そのため、除去したい周波数がわかっている場合はハイパスフィルタであれば少し上(ローパスフィルタの場合は少し下)の周波数を設定する必要があります。
コード例
フィルタの次数とカットオフ周波数が決まったので、実装していきます。
今回は10Hzと100Hzを合成した波形に対して、フィルタ次数4次・カットオフ周波数50Hzのローパスフィルタを適用してみます。
scipy.signal.butter
バターワースフィルタを作成します。
引数の説明は以下の通りです。
scipy.signal.butter(N, Wn, btype=’low’, analog=False, output=’ba’, fs=None)
N:フィルタ次数
Wn:カットオフ周波数(bandpass, bandstopの場合は2つ指定する)
btype:フィルタのタイプ(highpass, lowpass, bandpass, bandstop)
analog:アナログフィルタか
output:出力のタイプ(ba, sos, zpk)
fs:デジタルフィルタの場合、サンプリング周波数
使った感じanalogはFalseのままでよさそうでした。(というかアナログとデジタルの違いがわからないので)
また、ネット見るとカットオフ周波数をナイキスト周波数で割っている記述があったりしましたが、ナイキスト周波数で割る場合はデジタルフィルタでもfsの指定はいらないみたいでした。
カットオフ周波数をナイキスト周波数で割らないときはfsを指定する必要があります。
ドキュメント見ると出力はbaのままだと計算誤差があるようなので、sosを指定するほうが良いみたいです。
scipy.signal.sosfiltfilt
時間軸波形にフィルタを適用します。
scipy.signal.sosfiltfilt(sos, x, axis=- 1, padtype=’odd’, padlen=None)
sos:butterの戻り値
x:対象の波形
axis:どの軸に沿ってフィルタかけるか
padtype:パディングのタイプ?
padlen:パディングの長さ?
オプションの引数はよくわかんないですが、sosとxを指定すればなんとかなりました!
ソースコード
コード例です。
ローパスフィルタが無事に適用され、左の波形から高周波のノイズのようなものが除去できていることが確認できますね。
import matplotlib.pyplot as plt import numpy as np from scipy import signal # パラメータ設定 sample_freq = 32000 cutoff_freq = 50 filter_order = 4 sec = 1 f0 = 10 f1 = 100 # 合成派の作成 t = np.linspace(0, sec , sample_freq * sec) y = np.sin(2 * np.pi * f0 * t) + np.sin(2 * np.pi * f1 * t) # フィルタの作成・適用 sos = signal.butter(filter_order, cutoff_freq, 'lowpass', output='sos', fs=sample_freq) yf = signal.sosfiltfilt(sos, y) # グラフ描画 fig = plt.figure(figsize=(12,5), tight_layout=True) ax = plt.subplot(1, 2, 1) ax.plot(t, y) ax.set_xlabel('time') ax.set_ylabel('amplitude') ax.set_title('without Filter') ax = plt.subplot(1, 2, 2) ax.plot(t, yf) ax.set_xlabel('time') ax.set_ylabel('amplitude') ax.set_title('with Filter')
おしまい!
コメント