デジタルフィルタの実装

一般的なフルIIRフィルタ (IIRおよびFIRフィルタも含む) は時間ドメインでは入出力信号 $x, y$ にそれぞれフィルタ係数 $a, b$ をかけた畳み込みで実装できる: \begin{align} y[n] = \sum_{t=0}^{B} b_t x[n - t] - \sum_{t=1}^{A} a_t y[n - t] \end{align}

IIRフィルタが所望の周波数特性をもつように定義するにはZ変換・ラプラス変換・フーリエ変換を用いる。先の時間領域の等式の両辺を入力 $x[n]$ で割った入出力比を伝達関数(transfer function)と呼ぶ。さらに $x, y$ をZ変換とよばれる下記の置き換え \begin{align} X(z) = \sum_{t} x[t] z^{-t} \end{align} によりZ領域での伝達関数を考える: \begin{align} H(z) &= \frac{Y(z)}{X(z)} = \frac{X(z) \sum_{t=0}^{B} b_t z^{-t} - Y(z) \sum_{t=1}^{A} a_t z^{-t}}{X(z)} \nonumber \newline &= \sum_{t=0}^{B} b_t z^{-t} - H(z) \sum_{t=1}^{A} a_t z^{-t} \nonumber \newline &= \frac{\sum_{t=0}^{B} b_t z^{-t}}{1 + \sum_{t=1}^{A} a_t z^{-t}} \end{align} ここで次のような bilinear 変換によりラプラス領域の伝達関数 $H(s)$ にできる: \begin{align} \label{eq:laplace} s &= \frac{2}{T} \frac{z-1}{z+1}, \newline \end{align} ラプラス領域を介することで、さらにフーリエ変換として $s = j \omega$ を代入すれば周波数や位相の特性を得られる。絶対値をとれば $|H(j \omega)|$ が周波数ごとの音量 (db/octの傾きとかわかる), 位相をとれば $\angle{H(j \omega)}$ が周波数ごとの位相応答 (線形位相かなどわかる) を表す関数となる。

フィルタ設計の流れ

繰り返しになるが上記の各種変換により、デジタルの実装（離散時間領域）と、位相周波数応答の特性（フーリエ領域）を行き来することができる。ここでは触れないが、一時的にでてくるラプラス変換もアナログ回路との対応があり勉強すると楽しい。

ここまでの話を 4 から 1 に逆にたどると、デジタルフィルタの設計の流れとなる。つまり所望の位相・周波数特性を考え、逆変換により等価なフィルタ係数 $(a_t, b_t)$ を導出するという作業である。

Butterworth フィルタの理論

代表的なフィルタとしてButterworthフィルタが音響処理では普及している。他のフィルタと比べてカットオフ周波数のロールオフが緩やかでカットオフ後に荒ぶらない性質がある。アナログ回路に基づいているので多くの場合はラプラス領域での伝達関数が示される。デジタルフィルタとして実装する場合に必要である、Z変換した伝達関数は数式記号処理ソフトを使って式 \eqref{eq:laplace} を代入すると得られる。具体例は後述する実験結果か、以下の本の 8.10 節でも確認できる

より一般的な Butterworth フィルタについてはこの PDF が詳しい (余談だが、著者はNIのMassiveやReaktorの開発者らしい linkedin)

この辺の話は VA Filter Design の2章がよくまとまっている。それでは実際に使われるLPF/HPF/BPFを導出する。

1-pole Low Pass Filter (LPF) 6db/oct

天下り的に定義を覚えてもいいが、自分の中で納得する導出を考えてみる。

人間には音の高さや大きさが対数スケールで感じる(例えば倍の周波数が1オクターブ上として聞こえる)ので、指数的に周波数を倍にすると、音量が半分になるようなLPフィルタがほしい。そのためには周波数の逆数みたいな周波数応答がいい。 \begin{align} |H(j\omega)| &= \frac{1}{w}, \end{align} ただし周波数0のときに無限の音量になるのでやばい。そこで \begin{align} |H(j\omega)| &= \frac{1}{w + 1}, \end{align} とすれば周波数0では元の音量となるよう、先の関数を周波数 $w$ 軸に平行移動して音量が半減してゆく。ここまでは絶対値で議論してきたが、この性質をもつラプラス領域での伝達関数には以下のものが考えられる \begin{align} H(s) &= \frac{1}{s + 1}. \end{align} さらに $s = s’ / \omega_0$ を代入することで、今度は対数スケールで周波数軸を左右に平行移動できる。これは単なる変数変換なので、同じ代入式で任意のフィルタ伝達関数にカットオフ周波数 $\omega_0$ を導入できる。 \begin{align} H(s) = \frac{1}{\frac{s}{\omega_0} + 1} \end{align}

この式は 1 pole (分母 = 0の方程式が1つの解をもつ) であり、倍の周波数(1 octave上) で0.5倍の音量つまりデシベル(db)でいうと $20 \log_{10} 0.5 \approx -6$ db/octave で高域が減衰するという。

1-pole High Pass Filter (HPF) 6db/oct

LPFと同じような議論で、$s / \omega_0$ の代わりに二倍低い周波数が二倍小さいレベルとなるよう $\omega_0 / s$が使えることがわかる。代入するとこのような伝達関数が得られる: \begin{align} H(s) = \frac{\frac{s}{\omega_0}}{\frac{s}{\omega_0} + 1} \end{align}

2-pole Low Pass Filter (LPF) 12db/oct

単純に1-pole LPFを二回かけると 2 pole (分母 = 0の方程式が2つの解をもつ) で 12db/octave で高域が減衰するLPFになる。フィルタは畳み込みなので、そのラプラス変換は積つまり、もとのLPFの二乗になる。

単純に二回かけるかわりに、2-poleではレゾナンス $Q$ により周波数ピークをたたせることができる。 \begin{align} H(s) = \frac{1}{\frac{s^2}{\omega_0^2} + \frac{s}{\omega_0 Q} + 1} \end{align} TODO: なぜ Q でピークがたつかの説明 (天下り的に周波数応答みればよいが、直感的な説明があるとよい)。気になる人は VA Filter Designの 4.2 Resonanceを読むとよい。

2-pole High Pass Filter (HPF) 12db/oct

\begin{align} H(s) = \frac{\frac{s^2}{\omega_0^2}}{\frac{s^2}{\omega_0^2} + \frac{s}{\omega_0 Q} + 1} \end{align}

2-pole Band Pass Filter (BPF) 12db/oct

\begin{align} H(s) = \frac{\frac{s}{\omega_0 Q}}{\frac{s^2}{\omega_0^2} + \frac{s}{\omega_0 Q} + 1} \end{align}

LPF, HPFとBPFの伝達関数を全部足したら1になるな、と思った人は筋が良い。ラプラス変換には線形性があるので、実際の各フィルタの出力を足し合わせてもとに戻るように、何もしない伝達関数 $H(s)=1$ が現れる。

実装におけるフィルタ係数 (pythonによる自動導出)

ここまで得られたラプラス領域の伝達関数を式 \eqref{eq:laplace} を用いたZ変換することで冒頭の離散時間領域の多項式が求まり、実装に必要なデジタルフィルタの係数が得られる。1-poleくらいなら手計算してもいいが、2-poleになるとしんどいので面倒なことはPython (の数式処理ライブラリであるSymPy) にやらせよう。

# requires version '1.7.1
from sympy import *

s = Symbol('s')
z = Symbol('z')
Q = Symbol('Q')    # resonance
T = Symbol('T')    # sampling interval
w0 = Symbol('w0')  # cutoff freq

# z2s = 2 / T * (z - 1) / (z + 1)
s2z = 2 / T * (z - 1) / (z + 1)

def print_coeff(hs):
    hz = simplify(hs.subs(s, s2z))  # Z transform
    npole = degree(denom(hs), s)
    print("=== Transfer function ===")
    print("H(s) =", hs)  # transfer function in Laplace domain
    print("H(z) =", hz)  # transfer function in Z domain
    print("#pole =", npole)
    print("=== Filter coeffients ===")
    # FIR coeff
    dhz = collect(expand(denom(hz) * z ** -npole), z)
    nhz = collect(expand(numer(hz) * z ** -npole), z)
    a0 = dhz.coeff(z, 0)  # to normalize a0 = 1
    for i in range(npole + 1):
        print(f"b{i} =", nhz.coeff(z, -i) / a0)
    # IIR coeff
    for i in range(1, npole + 1):
        print(f"a{i} =", dhz.coeff(z, -i) / a0)


print("Filter: 1-pole LPF")
print_coeff(hs = 1 / (s / w0 + 1))
print()
print("Filter: 1-pole HPF")
print_coeff(hs = s / (s + w0))
print()
print("Filter: 2-pole LPF")
print_coeff(hs = 1 / (s**2 / w0**2 + s / w0 / Q + 1))
print()
print("Filter: 2-pole HPF")
print_coeff(hs = (s**2 / w0**2) / (s**2 / w0**2 + s / w0 / Q + 1))
print()
print("Filter: 2-pole BPF")
print_coeff(hs = (s / w0 / Q) / (s**2 / w0**2 + s / w0 / Q + 1))

Filter: 1-pole LPF
=== Transfer function ===
H(s) = 1/(s/w0 + 1)
H(z) = T*w0*(z + 1)/(T*w0*(z + 1) + 2*z - 2)
#pole = 1
=== Filter coeffients ===
b0 = T*w0/(T*w0 + 2)
b1 = T*w0/(T*w0 + 2)
a1 = (T*w0 - 2)/(T*w0 + 2)

Filter: 1-pole HPF
=== Transfer function ===
H(s) = s/(s + w0)
H(z) = 2*(z - 1)/(T*w0*(z + 1) + 2*z - 2)
#pole = 1
=== Filter coeffients ===
b0 = 2/(T*w0 + 2)
b1 = -2/(T*w0 + 2)
a1 = (T*w0 - 2)/(T*w0 + 2)

Filter: 2-pole LPF
=== Transfer function ===
H(s) = 1/(s**2/w0**2 + 1 + s/(Q*w0))
H(z) = Q*T**2*w0**2*(z + 1)**2/(Q*T**2*w0**2*(z + 1)**2 + 4*Q*(z - 1)**2 + 2*T*w0*(z - 1)*(z + 1))
#pole = 2
=== Filter coeffients ===
b0 = Q*T**2*w0**2/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
b1 = 2*Q*T**2*w0**2/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
b2 = Q*T**2*w0**2/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
a1 = (2*Q*T**2*w0**2 - 8*Q)/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
a2 = (Q*T**2*w0**2 + 4*Q - 2*T*w0)/(Q*T**2*w0**2 + 4*Q + 2*T*w0)

Filter: 2-pole HPF
=== Transfer function ===
H(s) = s**2/(w0**2*(s**2/w0**2 + 1 + s/(Q*w0)))
H(z) = 4*Q*(z - 1)**2/(Q*T**2*w0**2*(z + 1)**2 + 4*Q*(z - 1)**2 + 2*T*w0*(z - 1)*(z + 1))
#pole = 2
=== Filter coeffients ===
b0 = 4*Q/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
b1 = -8*Q/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
b2 = 4*Q/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
a1 = (2*Q*T**2*w0**2 - 8*Q)/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
a2 = (Q*T**2*w0**2 + 4*Q - 2*T*w0)/(Q*T**2*w0**2 + 4*Q + 2*T*w0)

Filter: 2-pole BPF
=== Transfer function ===
H(s) = s/(Q*w0*(s**2/w0**2 + 1 + s/(Q*w0)))
H(z) = 2*T*w0*(z - 1)*(z + 1)/(Q*T**2*w0**2*(z + 1)**2 + 4*Q*(z - 1)**2 + 2*T*w0*(z - 1)*(z + 1))
#pole = 2
=== Filter coeffients ===
b0 = 2*T*w0/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
b1 = 0
b2 = -2*T*w0/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
a1 = (2*Q*T**2*w0**2 - 8*Q)/(Q*T**2*w0**2 + 4*Q + 2*T*w0)
a2 = (Q*T**2*w0**2 + 4*Q - 2*T*w0)/(Q*T**2*w0**2 + 4*Q + 2*T*w0)

ある程度、展開できればいいやと思ってたが、完全に自動化できると思ってなかった…。ちなみに古いバージョンのsympyでやるとうまくいかないので注意。