前回に引き続き、今回は周波数応答のサンプルプログラムを公開します🐜
例としてローパスフィルタを扱います。
ローパスフィルタとハイパスフィルタ
以下にざっくり説明です。
| フィルタ | 特徴 |
|---|---|
| ローパスフィルタ | ・低い周波数帯域を通し、高い周波数帯域を遮断する ・ハイカットフィルタともいわれる ・FIRだと移動平均フィルタがよく使われる |
| ハイパスフィルタ | ・高い周波数帯域を通し、低い周波数帯域を遮断する ・ローカットフィルタともいわれる ・生データからローパス後のデータを減算することでも実現可能 |
意外と大事なのは、最後に書いてある
ハイパスフィルタ = （フィルタリング前のデータ） - （ローパスをかけたデータ）
です。
考えてみれば当たり前ですが、見落としがちな内容だと思います。ハイパスフィルタをさくっと作りたいときは、引き算するのが手っ取り早いですよ。
以降は、ローパスフィルタの代表である移動平均フィルタを題材に、実際にPythonで書いたサンプルプログラムを見ながら整理していきます。
開発環境
後述するサンプルプログラムを実行する開発環境は、次の通りです。
| 項目 | 内容 |
|---|---|
| OS | Windows 10 |
| Python | 3.7.6(Anaconda 4.8.1) |
| Numpy | 1.17.4 |
| Scipy | 1.3.2 |
| matplotlib | 3.1.1 |
サンプルプログラム
Pythonを使ったFFTのサンプルスクリプト(sample_fvtool.py)は以下の通りです。自由にコピペして、実際に動かしてみてください。
#!/usr/bin/python # -*- coding: utf-8 -*- import time import numpy as np from scipy import signal import matplotlib.pyplot as plt class TestFFT(): def __init__(self): self.msg = 'Test script for fvtool' # 設定値 self.fs = 250 # サンプリング周波数 def process(self): print(self.msg) lowpass = make_MA_filter(5) fvtool(lowpass, self.fs) plt.show() def fvtool(b, fs, a=1, worN=8192): # ナイキスト周波数計算 fn = fs / 2 # 周波数応答計算 w, h = signal.freqz(b, a, worN) x_freq_rspns = w / np.pi * fn y_freq_rspns = db(abs(h)) # 複素数→デシベル変換 # 群遅延計算 w, gd = signal.group_delay([b, a], worN) x_gd = w / np.pi * fn y_gd = gd # グラフ表示 plt.figure(figsize=(8, 6)) # 周波数応答プロット plt.subplot(2, 1, 1) plt.plot(x_freq_rspns, y_freq_rspns) plt.ylim(-70, 2) # MATLAB fvtool仕様 plt.ylabel('Amplitude [dB]') plt.grid() # 群遅延プロット plt.subplot(2, 1, 2) plt.plot(x_gd, y_gd) plt.ylim(0, len(b)) # MATLAB fvtool仕様 plt.xlabel('Frequency [Hz]') plt.ylabel('Group delay [samples]') plt.grid() def db(x, dBref=1): # デシベル変換 y = 20 * np.log10(x / dBref) return y def make_MA_filter(N): b = np.ones(N) / N a = [1 for i in range(N) ] return b/a if __name__ == '__main__': proc = TestFFT() proc.process()
実行結果は以下の通りです。
解説｜フィルタの要求仕様
下記のような移動平均フィルタを設計します。
| 項目 | 内容 |
|---|---|
| サンプリング周波数[Hz] | 250 |
| 除去したい周波数[Hz] | 50 |
電気・電子回路をいじっていると、よく信号に電源ノイズが重畳しますよね1。今回はそれを除去します。
解説｜周波数応答
基本的には下記ブログを参考に実装していけばよいので、軽くコメントする程度にします。
ディジタル信号処理｜周波数応答とフィルタ設計 - ari23の研究ノート
sample_fvtool.pyのmake_MA_filter関数
引数に応じた移動平均フィルタの係数をリストで戻します。遅延の都合上、奇数であることが望ましいです。
sample_fvtool.pyのfvtool関数
前回と同様に上から順に解説していきます。
def fvtool(b, fs, a=1, worN=8192): # ナイキスト周波数計算 fn = fs / 2 # 周波数応答計算 w, h = signal.freqz(b, a, worN) x_freq_rspns = w / np.pi * fn y_freq_rspns = db(abs(h)) # 複素数→デシベル変換
ナイキスト周波数はサンプリング周波数の半分で、周波数応答ではこのナイキスト周波数までが評価できます。詳細については、他の文献を参考にしてください2。
scipy.signal.freqzの第3引数のworNはデフォルト値でも問題ありませんが、ここではMATLABのfvtoolと同じ結果が得られるようにしています。
x軸は正規化角周波数だとわかりにくいので、単位変換します。
また、デシベル変換はFFTで用いたものと同じ関数です。
# 群遅延計算
w, gd = signal.group_delay([b, a], worN)
x_gd = w / np.pi * fn
y_gd = gd
こちらもx軸が正規化角周波数だとわかりにくいので、単位変換します。
最後に、周波数応答と群遅延の結果をグラフプロットします。
sample_fvtool.pyのグラフプロットの結果
以下に再掲し、上から順に解説します。
まず、1つ目は設計した周波数応答を示しています。50Hzと100Hzに谷があり、目的の50Hzを遮断できることが確認できます。
2つ目は設計したフィルタの群遅延を示しています。設計したフィルタは、FIRフィルタであり次数（タップ数）が5つなので、(5 - 1) /2 = 2サンプル遅れるのは正しいとわかります。
これで、周波数応答のサンプルプログラムを作ることができました。
おわりに
以前の記事を参考に、サンプルプログラムを作成しました。
設計したローパスフィルタはノッチフィルタに似ていますが、もちろん移動平均の係数の数を増やせば、より強めのフィルタになります。
ぜひ、試してみてくださいね（たくさんの山ができるような周波数応答になるはずです）。
また、公開したスクリプトは、FIRフィルタしか想定していないので、IIRフィルタを設計する際は適宜修正してください。
参考になれば幸いです(^^)
