n次のバターワースBiQuadフィルターを設計する

信号処理,数学

birdhouse.hateblo.jp
以前のBiQuadフィルターを改良して、n次のバターワース型に対応したデジタルフィルタを考えてみます。

バターワース型フィルターとは

バターワース型フィルタとは、カットオフ周波数で-3[dB]減衰しn×-6[dB/Octave]の減衰特性を持ったフィルターのことを指します。次数によってより急峻に減衰していくことになります。
バターワースフィルタ – Wikipedia

つまり、2次でカットオフが1000[Hz]なら周波数特性は
1000[Hz]で-3[dB]、2000[Hz]で-12[dB]するということ。

n次のバターワース型フィルタの伝達関数は以下の式で表されます。
\displaystyle H(s) = \frac{G_0}{\displaystyle \prod_{k=1}^{n}\frac{(s-s_k)}{\omega_c}}
なおs_k
\displaystyle s_k=\omega_c e^\frac{j(2k+n-1)\pi}{2n}
と、複素数形式での表現となってます。これをバターワース多項式といいます。

これをそのまま複素数形式でもいいけど、複素共益な極同士をかけることで実数形式で書くこともできます。
さらに\omega_c=1として正規化すると、正規化バターワース多項式とすることができます。
H(s)=\frac{G_0}{B_n(s)}とし、正規化バターワース多項式B_n(s)
B_n(s)=\displaystyle \prod_{k=1}^{\frac{n}{2}} \left(s^2-2s\cos(\frac{2k+n-1}{2n}\pi)+1\right)  ...n=even
B_n(s)=\displaystyle (s+1)\prod_{k=1}^{\frac{n}{2}} \left(s^2-2s\cos(\frac{2k+n-1}{2n}\pi)+1\right)  ...n=odd
と書き表すことができます。

つまり、次数が増えても2次または1次のバターワースフィルタの縦続接続ですべて表現することができるということです。
この基本形となる2次または1次のバターワースフィルタをBiQuadフィルタとし、これを縦続接続することでn次を実現することにします。

BiQuad係数の導出

1次BiQuad係数の導出

双1次変換とプリワーピング

1次の場合、バターワース多項式B_n(s)(s+1)となります。
なのでLPFの場合、伝達関数H(s)=\frac{1}{s+1}です。
これを双1次変換していきます。
双一次変換の式は
\displaystyle S=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}
でした。正規化バターワース多項式では\omega_c=1としていますが、実際にはカットオフ周波数が入るのでプリワーピングした周波数で割ります。
プリワーピングの式は
\displaystyle \omega_a=\frac{2}{T}\tan(\frac{\omega_0}{2})
よって
\displaystyle s=\frac{1}{\tan(\frac{\omega_0}{2})}\frac{1-z^{-1}}{1+z^{-1}}
半角の公式から*1
\displaystyle s=\frac{1+\cos(\omega_0)}{\sin(\omega_0)}\frac{1-z^{-1}}{1+z^{-1}}
この式をもとに1を
\displaystyle 1=\frac{1+\cos(\omega_0)}{1+\cos(\omega_0)}\frac{1+z^{-1}}{1+z^{-1}}
として、s1の共通因数である\displaystyle \frac{1+\cos(\omega_0)}{1+z^{-1}}を除く
\displaystyle
\left[
\begin{split}
s&=&\frac{1-z^{-1}}{\sin(\omega_0)}\\
1&=&\frac{1+z^{-1}}{1+\cos(\omega_0)}
\end{split}\right.
それぞれに\sin^2(\omega_0)をかける*2
\displaystyle
\left[
\begin{split}
s&=&(1-z^{-1})\sin(\omega_0)\\
1&=&(1+z^{-1})(1-\cos(\omega_0))
\end{split}\right.
展開して整理すると
\displaystyle
\left[
\begin{split}
s&=\sin(\omega_0)&-\sin(\omega_0)z^{-1}\\
1&=(1-\cos(\omega_0))&+(1-\cos(\omega_0))z^{-1}
\end{split}\right.
となりました。

1次LPF係数の導出

1次バタワースのLPF伝達関数\displaystyle H(s)=\frac{1}{s+1}に先程の双一次変換でもとめた式を代入していく
\displaystyle
\begin{split}
H(z)&=\frac{(1-\cos(\omega_0))+(1-\cos(\omega_0))z^{-1}}{\sin(\omega_0)-\sin(\omega_0)z^{-1}+(1-\cos(\omega_0))+(1-\cos(\omega_0))z^{-1}}\\
&=\frac{(1-\cos(\omega_0))+(1-\cos(\omega_0))z^{-1}}{(1-\cos(\omega_0)+\sin(\omega_0))+(1-\cos(\omega_0)-\sin(\omega_0))z^{-1}}
\end{split}
よって
\displaystyle
\left[
\begin{split}
a_0&=1-\cos(\omega_0)+\sin(\omega_0)\\
a_1&=1-\cos(\omega_0)-\sin(\omega_0)\\
a_2&=0\\
b_0&=1-\cos(\omega_0)\\
b_1&=1-\cos(\omega_0)\\
b_2&=0
\end{split}
\right.

1次HPF係数の導出

1次バタワースのHPF伝達関数\displaystyle H(s)=\frac{s}{s+1}に先程の双一次変換でもとめた式を代入していく
\displaystyle
\begin{split}
H(z)&=\frac{\sin(\omega_0)-\sin(\omega_0)z^{-1}}{(1-\cos(\omega_0)+\sin(\omega_0))+(1-\cos(\omega_0)-\sin(\omega_0))z^{-1}}
\end{split}
よって
\displaystyle
\left[
\begin{split}
a_0&=1-\cos(\omega_0)+\sin(\omega_0)\\
a_1&=1-\cos(\omega_0)-\sin(\omega_0)\\
a_2&=0\\
b_0&=\sin(\omega_0)\\
b_1&=-\sin(\omega_0)\\
b_2&=0
\end{split}
\right.

2次BiQuad係数の導出

双1次変換とプリワーピング

途中までは1次のときと同じです。
\displaystyle s=\frac{1}{\tan(\frac{\omega_0}{2})}\frac{1-z^{-1}}{1+z^{-1}}
として
\displaystyle
\left[
\begin{split}
s&=\frac{1}{\tan(\frac{\omega_0}{2})}\frac{1-z^{-1}}{1+z^{-1}}\\
s^2&=\frac{1}{\tan^2(\frac{\omega_0}{2})}\frac{(1-z^{-1})^2}{(1+z^{-1})^2}
\end{split}
\right.
半角の公式より*3
\displaystyle
\left[
\begin{split}
s&=\frac{1+\cos(\omega_0)}{\sin(\omega_0)}\frac{1-z^{-1}}{1+z^{-1}}\\
s^2&=\frac{1+\cos(\omega_0)}{1-\cos(\omega_0)}\frac{(1-z^{-1})^2}{(1+z^{-1})^2}
\end{split}
\right.
これを変形して
\displaystyle
\left[
\begin{split}
1&=\frac{1+\cos(\omega_0)}{1+\cos(\omega_0)}\frac{(1+z^{-1})^2}{(1+z^{-1})^2}\\
s&=\frac{1+\cos(\omega_0)}{\sin(\omega_0)}\frac{(1-z^{-1})(1+z^{-1})}{(1+z^{-1})^2}\\
s^2&=\frac{1+\cos(\omega_0)}{1-\cos(\omega_0)}\frac{(1-z^{-1})^2}{(1+z^{-1})^2}
\end{split}
\right.
共通因数\displaystyle \frac{1+\cos(\omega_0)}{(1+z^{-1})^2}を除き
\displaystyle
\left[
\begin{split}
1&=\frac{(1+z^{-1})^2}{1+\cos(\omega_0)}\\
s&=\frac{(1-z^{-1})(1+z^{-1})}{\sin(\omega_0)}\\
s^2&=\frac{(1-z^{-1})^2}{1-\cos(\omega_0)}
\end{split}
\right.
すべての項にsin^2(\omega_0)をかけて*4
\displaystyle
\left[
\begin{split}
1&=(1+z^{-1})^2(1-\cos(\omega_0)\\
s&=(1-z^{-1})(1+z^{-1})\sin(\omega_0)\\
s^2&=(1-z^{-1})^2(1+\cos(\omega_0))
\end{split}
\right.
展開して整理すると
\displaystyle
\left[
\begin{split}
1&=(1-\cos(\omega_0))&+(2-2\cos(\omega_0))z^{-1}&+(1-\cos(\omega_0))z^{-2}\\
s&=\sin(\omega_0)&&-\sin(\omega_0)z^{-2}\\
s^2&=(1+\cos(\omega_0))&-(2+2\cos(\omega_0))z^{-1}&+(1+\cos(\omega_0))z^{-2}
\end{split}
\right.

2次LPF係数の導出

バターワース多項式B_n(s)=\displaystyle \prod_{k=1}^{\frac{n}{2}} \left(s^2-2s\cos(\frac{2k+n-1}{2n}\pi)+1\right)  ...n=even
ここで-2s\cos(\frac{2k+n-1}{2n}\pi)に注目し
-2cos(\frac{2k+n-1}{2n}\pi) = A(n,k)と置くと
B_n(s)=\displaystyle \prod_{k=1}^{\frac{n}{2}} \left(s^2-A(n,k)s+1\right)  ...n=even
と表せる。つまり、2次LPFの伝達関数H(s)
\displaystyle H(s)=\frac{1}{s^2-A(n,k)s+1}となる
よって、先程の双一次変換した式を代入していくと
\displaystyle
\begin{split}
H(z)&=\frac{(1-\cos(\omega_0))+(2-2\cos(\omega_0))z^{-1}+(1-\cos(\omega_0))z^{-2}}{(1+\cos(\omega_0))-(2+2\cos(\omega_0))z^{-1}+(1+\cos(\omega_0))z^{-2}+A(n,k)(\sin(\omega_0)-\sin(\omega_0)z^{-2})+(1-\cos(\omega_0))+(2-2\cos(\omega_0))z^{-1}+(1-\cos(\omega_0))z^{-2}}\\
&=\frac{(1-\cos(\omega_0))+(2-2\cos(\omega_0))z^{-1}+(1-\cos(\omega_0))z^{-2}}{(2+A(n,k)\sin(\omega_0))-4\cos(\omega_0)z^{-1}+(2-A(n,k)\sin(\omega_0))z^{-2}}
\end{split}
つまり
\displaystyle
\left[
\begin{split}
a_0&=2+A(n,k)\sin(\omega_0)\\
a_1&=-4\cos(\omega_0)\\
a_2&=2+A(n,k)\sin(\omega_0)\\
b_0&=1-\cos(\omega_0)\\
b_1&=2-2\cos(\omega_0)\\
b_2&=1-\cos(\omega_0)
\end{split}
\right.
すべてを2で割る
\displaystyle
\left[
\begin{split}
a_0&=1+\frac{A(n,k)}{2}\sin(\omega_0)\\
a_1&=-2\cos(\omega_0)\\
a_2&=1+\frac{A(n,k)}{2}\sin(\omega_0)\\
b_0&=\frac{1-\cos(\omega_0)}{2}\\
b_1&=1-\cos(\omega_0)\\
b_2&=\frac{1-\cos(\omega_0)}{2}
\end{split}
\right.
2次ではA(n,k)n=2k=1を代入してA(n,k)\fallingdotseq1.41421...

2次HPF係数の導出

LPFと同様に考えると2次HPFの伝達関数H(s)
\displaystyle H(s)=\frac{s^2}{s^2-A(n,k)s+1}となる
よって、先程の双一次変換した式を代入していくと
\displaystyle
\begin{split}
H(z)&=\frac{(1+\cos(\omega_0))-(2+2\cos(\omega_0))z^{-1}+(1+\cos(\omega_0))z^{-2}}{(1+\cos(\omega_0))-(2+2\cos(\omega_0))z^{-1}+(1+\cos(\omega_0))z^{-2}+A(n,k)(\sin(\omega_0)-\sin(\omega_0)z^{-2})+(1-\cos(\omega_0))+(2-2\cos(\omega_0))z^{-1}+(1-\cos(\omega_0))z^{-2}}\\
&=\frac{(1+\cos(\omega_0))-(2+2\cos(\omega_0))z^{-1}+(1+\cos(\omega_0))z^{-2}}{(2+A(n,k)\sin(\omega_0))-4\cos(\omega_0)z^{-1}+(2-A(n,k)\sin(\omega_0))z^{-2}}
\end{split}
つまり
\displaystyle
\left[
\begin{split}
a_0&=2+A(n,k)\sin(\omega_0)\\
a_1&=-4\cos(\omega_0)\\
a_2&=2+A(n,k)\sin(\omega_0)\\
b_0&=1+\cos(\omega_0)\\
b_1&=-2-2\cos(\omega_0)\\
b_2&=1+\cos(\omega_0)
\end{split}
\right.
すべてを2で割る
\displaystyle
\left[
\begin{split}
a_0&=1+\frac{A(n,k)}{2}\sin(\omega_0)\\
a_1&=-2\cos(\omega_0)\\
a_2&=1+\frac{A(n,k)}{2}\sin(\omega_0)\\
b_0&=\frac{1+\cos(\omega_0)}{2}\\
b_1&=-1-\cos(\omega_0)\\
b_2&=\frac{1+\cos(\omega_0)}{2}
\end{split}
\right.
2次ではA(n,k)n=2k=1を代入してA(n,k)\fallingdotseq1.41421...

周波数特性を算出する

導いた伝達関数を持つフィルターの周波数特性を調べる場合にはH(z)の式にz=e^{j\omega}=e^{2j\pi fT}を代入すれば計算できる。
fの範囲を0 < f < \frac{F_s}{2}の範囲で変化させればよい。
\displaystyle
\begin{split}
H(e^{j\omega})&=\frac{b_0+b_1e^{-j\omega}+b_2e^{-2j\omega}}{a_0+a_1e^{-j\omega}+a_2e^{-2j\omega}}\\
\end{split}
オイラーの公式から*5
\displaystyle
\begin{split}
H(e^{j\omega})&=\frac{b_0+b_1(\cos\omega-j\sin\omega)+b_2(\cos2\omega-j\sin2\omega)}{a_0+a_1(\cos\omega-j\sin\omega)+a_2(\cos2\omega-j\sin2\omega)}
\end{split}
複素数なのでH(e^{j\omega})の絶対値を取り、dBにする。
20\log_{10}|H(e^{j\omega})|

で、実装したフィルターの周波数特性を出してみました。
カットオフ周波数は1[kHz]でサンプリング周波数は48[kHz]です。
次数はnを1~4まで確認しました。

f:id:tsubakurame-1913:20210615174112p:plain
LPF特性
f:id:tsubakurame-1913:20210615174116p:plain
HPF特性

カットオフの-3[dB]とn×-6[dB/Octave]の減衰が実現されてます。

余談

多分、これを応用すればn次のチェビシェフとかも作れるはず。多分。
チェビシェフ多項式の方が面倒そうだし、仕事でとりあえずバターワースを使う感じだったので作ってはないけどね。
そのうち作ってはみたいなとは思うけど、めんどそう。