pythonでデジタルフィルタ FIRフィルタ

Pythonでデジタルフィルタを書きたい時が、まれにある。
でも自分で計算したフィルタの係数が、
本当にその周波数特性を持っているのか確認したくなる。

そこでpythonのscipyにあるfreqzを使用して確認してみる。

FIRフィルタとIIRフィルタ

宮崎のホームページ
デジタルフィルタはこのHPの図にあるように、入力信号や出力信号をクロックで遅延させたものに係数をかけたものの和をとっており
FIRでもIIRでも一般化して伝達関数は以下のように書ける。

$H(z)=\frac{B(z)}{A(z)} = \frac{b[0] + b[1]z^{-1} + ・・・+ b[M]z^{-M}}{a[0] + a[1]z^{-1} + ・・・+ a[N]z^{-N}}$

これの係数を得ることがデジタルフィルタを設計することとなる。
今回はこれがある前提で、周波数特性を得ることを目的にする。

FIRフィルタ設計

とりあえずフィルタ係数を得る。
設計の話はたくさん出てくる。
気が向いたらPythonでなんか書いてみる。
窓関数法の設計仕様

この右のグラフのようなものを得る予定。

scipy.signal.freqz

scipy.signal.freqz — SciPy v1.8.0 Manual
伝達関数を得られたら、freqzの式の係数に対応させる。

$z^{-M}$

と

$e^{-jwM}$

が対応する。

from signal import signal
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

a = [1]

b = [4.036409219416838e-03,
     -2.208366021329011e-18,
     -2.346162858786038e-02,
     -3.367761421639100e-02,
     7.203267485172524e-02,
     2.840739092673040e-01,
     4.000000000000000e-01,
     2.840739092673040e-01,
     7.203267485172524e-02,
     -3.367761421639100e-02,
     -2.346162858786038e-02,
     -2.208366021329011e-18,
     4.036409219416838e-03]

fs = 1e9  # 1GHz

f, h = signal.freqz(b, a, fs=fs)

plt.semilogx(f, 20*np.log10(abs(h)))
plt.vlines(2e8, -200, 10, color="red")
plt.hlines(-6, 1, 1e9, color="green")
plt.xlim(1e8, 5e8)
plt.ylim(-200, 4)
plt.grid(which="both")
plt.show()