双二次フィルタの係数計算方法
双二次フィルタの係数計算方法を勉強したので備忘録として残しておきます。下記のサイトが大変参考になりました。他のIIRフィルタに関する記事もぜひご覧ください。
まずは簡単に今回対象とする３つのフィルタの特徴をまとめます。チェビシェフには遮断領域にリプルを持たせるチェビシェフⅡ型がありますが今回は触れません。
バタワース特性：最も代表的なフィルタ特性で、他の型よりも簡単な式で表される。振幅特性が優れており（最大平坦特性）、次数を上げていくと理想フィルタ特性に近づいていく。通過域の特性が良いため音声処理ではよく使用される。
チェビシェフ特性：通過域に一定振幅のリプルを持つのが最大の特徴である。バタワースフィルタが周波数の増加に従ってなだらかに減衰していくのに対して、チェビシェフフィルタはカットオフ周波数付近で急峻な特性を持っている。
ベッセル特性：群遅延特性に優れており、通過域の位相特性が周波数に対して直線的に変化する（位相直線）。単位ステップ応答の立ち上がりおよびリンギングが最も良好である。その代わりカットオフ特性がなだらかとなる。
設計には大きく分けて4つのStepがあります。Step1：フィルタの原型となる基準LPF（低域原型フィルタ）の作成、Step2：双一次変換の際の周波数軸のゆがみに対応するのための周波数変換（予補償、プリワーピング）、Step3：目的のカットオフ周波数とフィルタ形式に変換する周波数・型変換、Step4：アナログ→デジタル空間に変換する双一次s-z変換で成り立っています。
Step1 基準ＬＰＦ設計
まずはアナログの基準LPFを設計します。説明の単純化のために双二次フィルタ（２次IIRフィルタ）の場合に限った伝達関数を示します。すべての次数に一般化した式は先ほどの参考リンクに載っています。各フィルタ型の基準LPFの伝達関数は下記の通りです。
・バタワース型
・チェビシェフ型
A：通過域リプル＝遮断周波数におけるゲイン
・ベッセル型
Step2 プリワーピング（予補償）
後段の双一次s-z変換は周波数領域で非線形変換が発生するため、アナログフィルタとデジタルフィルタの間で周波数がイコールになりません。設計しようとしているデジタルフィルタのカットオフ周波数に対して、前もってその逆変換をアナログフィルタに対して実施します。
T：サンプリング周期
Ωc：カットオフ角周波数（デジタル領域）
Fs：サンプリング周波数
Fc：カットオフ周波数
ωc：カットオフ角周波数（アナログ領域）
Step3 周波数・型変換
作成した基準LPFに対して周波数とフィルタタイプの変換を行います。ここでの角周波数ωcはStep2 プリワーピングで求めた角周波数です。作成したいフィルタタイプに合わせて、Step1で求めた基本LPF伝達関数のsを置き換えます。これを実施する事で基準LPFのカットオフ周波数がプリワーピングで求めた角周波数にシフトします。
・LPF
・HPF
・BPF
・BRF
ωc1：低域側カットオフ角周波数
ωc2：高域側カットオフ角周波数
BPFとBRFの場合は次数が4次になってしまうため、双二次フィルタとして設計するには因数分解して二つの二次式の積の形で表現する必要があります。真正面から計算すると大変なので、最初の参考リンクの表2-12を参考にするといいと思います。
Step4 双一次s-z変換
最後にアナログフィルタをデジタルフィルタに変換する双一次s-z変換を行います。下記の式をStep3で求めた伝達関数に代入します。式を整理して分母の定数項が1になるように分子分母を割ってやれば（正規化）双二次フィルタの伝達関数の形になるはずです。
ちなみに過去記事でも紹介していますが、ピーキングフィルタとかシェルビングフィルタとか他のフィルタの伝達関数はこちらのpdfが見やすいです。
https://wsignal.sakura.ne.jp/onsei2007/kousaku/column.pdf
おまけ フィルタ特性の確認
求めた伝達関数からフィルタ特性を求めるpythonプログラムを作成しました。設計した所望のフィルタになりましたか？私は計算ミスしまくってなかなか思い通りのフィルタが出来ませんでした。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 48000
a2 = 0.1716
a1 = 0
a0 = 1
b2 = 0.2929
b1 = 0.5858
b0 = 0.2929
a = [a0, a1, a2]
b = [b0, b1, b2]
w, h = signal.freqz(b, a, 65536)
f = w / (2 * np.pi) * fs
amp = 20 * np.log10(abs(h))
plt.figure(figsize=(10, 4))
plt.xlabel('Freq [Hz]')
plt.ylabel('Amp [dB]')
plt.xlim(20, fs)
plt.ylim(-120, 20)
plt.grid(which='both', axis='both')
plt.plot(f, amp, color="tab:blue")
plt.show()
