#1. はじめに
MOOGをはじめとしたシンセサイザーに使われているラダーフィルタは，下図(Fig. 1.1)のようなシステムとして表される[1, 2]．
Fig.1.1: ラダーフィルタのアナログブロック線図．
ここで，$G(s, \omega_c)$は，カットオフ角周波数$\omega_c$の1次系のローパスフィルタの伝達関数である．上図からわかるように，ラダーフィルタは，1次系フィルタが4つ連結したフィードバックシステムとなっている．また，フィードバックゲイン$k$を調整することで，レゾナンスを発生させることが可能である．
本稿では，Fig. 1.1のシステムをコンピュータプログラムとして実装可能なデジタルシステムへと変換するまでの道筋を記述した．
#2. 1次系のアナログローパスフィルタ
まず，以下のような，1次系のローパスフィルタ特性を示すRC回路(Fig. 2.1)を考える．
Fig. 2.1: RCローパスフィルタ．
ここで，入出力電圧をそれぞれ$x$，$y$とすると，その関係は次の微分方程式によって記述することができる．
\dot{y}(t) = -\frac{1}{RC}y(t) + \frac{1}{RC}x(t) \tag{2.1}
ここで，$\frac{1}{RC} = \omega_c$，$y(0) = 0$とし，上式をラプラス変換すると次のようになる．
sY(s) = -\omega_cY(s) + \omega_cX(s) \tag{2.2}\\
従って，システムの伝達関数は次のように表現される．
G(s, \omega_c) = \frac{\omega_c}{s + \omega_c} \tag{2.3}\\
この伝達関数の周波数特性を以下(Fig. 2.2)に示す．
Fig. 2.2: 1次系アナログローパスフィルタの周波数特性．$f_c$はカットオフ周波数．
上記の伝達関数(Eq. 2.3)は，確かに，ローパスフィルターとしての特性を有することが確認できる．
#3. 1次系のデジタルローパスフィルタ
続けて，先の1次系のアナログローパスフィルタの伝達関数$G(s, \omega_c)$(Eq. 2.3)をサンプリング周波数$Fs$の下で，双一次Z変換(Eq. 3.1)を用いてデジタル伝達関数$H(z, \omega_c)$(Eq. 3.2)に変換すると，次のような形になる．
s\rightarrow2Fs\frac{1-z^{-1}}{1+z^{-1}} \tag{3.1}\\
H(z, \omega_c) = \frac{b_0 + b_1 z^{-1}}{1 + a_1 z^{-1}} \tag{3.2}\\
ここで，上記の$a_n$，$b_n$はフィルタ係数と呼ばれる値であり，今回の場合は下記のようにな値をとる．
\begin{align}
b_0 &= \frac{\omega_c}{\omega_c + 2Fs} \tag{3.3.1}\\
b_1 &= b_0 \tag{3.3.2}\\
a_1 &= \frac{\omega_c - 2Fs}{\omega_c + 2Fs} \tag{3.3.3}\\
\end{align}
この伝達関数(Eq. 3.2)の周波数特性を以下(Fig. 3.1)に示す．
Fig. 3.1: 1次系デジタルローパスフィルタの周波数特性($Fs = 48000\ Hz$)．
さらに，$z^{-n}$を掛けることは，$n$サンプル時間遅延させるということに等しいので，求めた$H(z, \omega_c)$(Eq. 3.2)に対する入出力それぞれ$x$，$y$として，次のような差分方程式の形で入出力関係を記述することができる(Eq. 3.4)．
y[n] = b_0x[n] + b_1x[n-1] - a_1y[n-1] \tag{3.4}\\
ここまでの操作を整理すると，まず，連続時間軸におけるアナログローパスフィルタの入出力関係を示した微分方程式に対し，ラプラス変換を施すことでその伝達関数を求めた．続いて，Z変換によりアナログ伝達関数を離散化した後に，再び時間軸での表現に戻した．この時，システムの入出力関係は差分方程式，すなわち，離散時間軸での表現として記述されているため，コンピュータプログラムとして実装が可能となる．
#4. ラダーフィルタの伝達関数を求める
ここで，改めてラダーフィルタのブロック線図を下記に示す(Fig. 4.1)．
Fig. 4.1: デジタルラダーフィルタのブロック線図．
ブロック線図の等価交換を用いると，上記のブロック線図全体の伝達関数$H'(z)$は1次のデジタルローパスフィルタの伝達関数$H(z, \omega_c)$(Eq. 3.2)を用いて，次のような形として表すことができる(Eq. 4.1)．
H'(z, \omega_c) = \frac{\{H(z, \omega_c)\}^4}{1 + k\{H(z, \omega_c)\}^4} \tag{4.1}
上式を整理すると，次のような形にすることができる(Eq. 4.2)．
H'(z, \omega_c) = \frac{\displaystyle \sum_{n = 0}^{4}b_n'z^{-n}}{1 + \displaystyle \sum_{n = 1}^{4}a_n'z^{-n}} \tag{4.2}
ここで，各フィルタ係数$b_n'$，$a_n'$は，1次のデジタルローパスフィルタの伝達関数$H(z, \omega_c)$(Eq. 3.2)のフィルタ係数$b_n$，$a_n$(Eq. 3.3)を用いて，次のように表すことができる(Eq. 4.3)．
\begin{align}
b_n' &= \frac{{}_4 C_n b_0^{(4-n)} b_1^n}{1 + kb_0^4} \tag{4.3.1}\\
a_n' &= \frac{{}_4 C_n (a_0^n + k b_0^{(4-n)} b_1^n)}{1 + kb_0^4} \tag{4.2.2}\\
\end{align}
この伝達関数の周波数特性(Fig. 4.2)と極配置(Fig. 4.3)は次のようになる．
Fig. 4.2: 伝達関数$H'(z)$(Eq. 4.2)の周波数特性($Fs = 48000\ Hz, k = 3$)．
Fig. 4.3: 伝達関数$H'(z)$(Eq. 4.2)の極配置($Fs = 48000\ Hz, k = 3$)
$k = contant$あれば，レゾナンスの強度は一定かつ，極もz平面の単位円内に収まるためシステムは安定であることがわかる．
#5. リアルタイム処理での問題点
先の伝達関数$H'(z, \omega_c)$を用いたモデルには，1つ問題点がある．それは，カットオフ$\omega_c$や，フィードバックゲイン$k$が変化するたびに，全てのフィルタ係数(Eq. 4.3)を再計算しなければならないという点である．LFOやフィルタENVを用いて常にパラメータ変化させることを考えると，その計算コストは無視できない[3]．
ここで，ブロック線図の等価交換を適用せずにラダーフィルタの出力を求めることを考えてみる．システムの入出力をそれぞれ$u$，$y$とすると，等価交換適用前のブロック線図は，Fig. 4.1と同様に，下記のようになる．
Fig. 5.1: デジタルラダーフィルタのブロック線図．
Fig. 5.1中の$H(z, \omega_c)$は，Eq. 3.4の差分方程式として表現できるため，出力$y$は次のように計算できる．
\begin{align}
\begin{cases}
x_0[n] &= u[n] - ky[n] \\
x_1[n] &= b_0x_0[n] + b_1x_0[n-1] - a_1x_1[n-1]\\
x_2[n] &= b_0x_1[n] + b_1x_1[n-1] - a_1x_2[n-1]\\
x_3[n] &= b_0x_2[n] + b_1x_2[n-1] - a_1x_3[n-1]\\
y[n] &= b_0x_3[n] + b_1x_3[n-1] - a_1y[n-1]\\
\end{cases}
\tag{5.1}
\end{align}
このモデルに従えば，カットオフ$\omega_c$が変化しても，1次系ローパスのフィルタ係数$b$, $a$(Eq. 3.3)を再計算するだけで済む．しかし，この場合フィードバックループに遅延器を含んでいない状態（ディレイフリーループ）を残すため，$x_0[n]$を求めるために$y[n]$の値が必要となり，また，$y[n]$を求めるためには$x_0[n]$の値が必要となるため，事実上その実装は不可能である．
#6. ディレイフリーループの対処
先の場合，システムのフィードバックループ経路が遅延を含んでいないことが問題となった．従って，一番単純な対処方法は，ループ経路に1つ遅延器を挿入することであるである．すなわち，次のようなブロック線図を考える(Fig. 6.1)．
Fig. 6.1: ループに遅延を挿入したデジタルラダーフィルタのブロック線図．
このブロック線図からなるシステムの周波数特性(Fig. 6.2)と極配置(Fig. 6.3)は次のようになる．
Fig. 6.2: Fig. 6.1システムの周波数特性($Fs = 48000\ Hz, k = 3$)．
Fig. 6.3: Fig. 6.1システムの極配置($Fs = 48000\ Hz, k = 3$)．
遅延器を挿入することで，実装上の問題は回避することができたが，上のボード線図と極配置からわかるように，今度は，レゾナンスを調整するフィードバックゲイン$k$に対する応答性が安定しないという問題が生じる．それ以前に，カットオフ周波数によっては一部の極がz平面上の単位円外に位置してしまう．この時，システムは不安定になるため，Fig. 6.1システムは可変フィルタとして用いるのには適さない．
#7. フィードバックゲインに対する応答性を安定させる
ラダーフィルタのフィードバックループ経路に遅延器を挿入すると，フィードバックゲイン$k$に対する応答性が安定しないという問題に対する対処法は，T. Stilsonらの論文[1]が参考になる．結論だけを示すと，1次系ローパスのアナログ伝達関数をZ変換によって離散化する際に，双一次Z変換ではなく下記のような変換(Eq. 7.1)を用いることで，フィードバックゲイン$k$に対する応答性を改善できる．
s\rightarrow 1.3Fs\frac{1-z^{-1}}{1+0.3z^{-1}} \tag{7.1}\\
この時の1次系ローパスのデジタル伝達関数とフィルタ係数をそれぞれ$H''(z, \omega_c)$，$b_n'', a_n''$として，次のように表すことができる．
H''(z, \omega_c) = \frac{b_0'' + b_1'' z^{-1}}{1 + a_1'' z^{-1}} \tag{7.2}\\
\begin{align}
b_0'' &= \frac{\omega_c}{\omega_c + 1.3Fs} \tag{7.3.1}\\
b_1'' &= \frac{0.3\omega_c}{\omega_c + 1.3Fs} \tag{7.3.2}\\
a_1'' &= \frac{0.3\omega_c - 1.3Fs}{\omega_c + 1.3Fs} \tag{7.3.3}\\
\end{align}
ここで，上記の1次系ローパスのデジタル伝達関数$H''(z)$を用いてFig. 6.1のループ経路に遅延を含むデジタルラダーフィルタのブロック線図を書き直すと次のようになる(Fig. 7.1)．
Fig. 7.1: Eq. 7.1によってデジタル化したラダーフィルタ（遅延あり）
このシステムのシステムの周波数特性(Fig. 7.2)と極配置(Fig. 7.3)は次のようになる．
Fig. 7.2: Fig 7.1システムの周波数特性($Fs = 48000\ Hz, k = 3$)．
Fig. 7.3: Fig 7.1システムの極配置($Fs = 48000\ Hz, k = 3$)．
Eq. 7.1のZ変換を用いることで，フィードバックゲイン$k = constant$下でのレゾナンス応答をほぼ一定に保ちつつ，システムも安定な状態に保つことが可能となった．
#8. おわりに
今回は，ディレイフリーループに遅延器を挿入することでその問題に対処したが，遅延器の挿入なしにディレイフリーループを近似する手法も存在する[2]．
記事を書き上げて思ったのですが，文中に具体的なプログラミング要素が皆無になってしまったので，そのうちJUCEでの具体的な実装について書くので今回のところは許してていただきたく．
#9. 参考文献
- T. Stilson, J. O. Smith: Analyzing the Moog VCF with Considerations for Digital Implementation, in Proc. Int. Computer Music Conf., Hong Kong, China, 1996, pp. 398−401. [Link]
- V. Zavalishin: THE ART OF VA FILTER DESIGN, 2015. [Link]
- 越田 俊介, 阿部 正英, 川又 政征: 可変ディジタルフィルタ -設計から応用まで-, IEICE Fundamentals Review, vol.10, no.1, 2016, pp.34-45. [Link]
- 三谷政昭: 信号解析のための数学, 森北出版, 1990.
- 大類重範: ディジタル信号処理, 日本理工出版会, 2001.
- 佐藤和也, 平本和彦，平田研二: はじめての制御工学, 講談社, 2010.