DTMerのためのデジタル信号処理 ～IIRなEQを作ろう～

Published: 2025-03-16

はじめに

こんにちは！Kashiwadeです。普段はオーケストラ系の曲を作っています。

最近円安で辛いですね…。全然プラグイン買えてないです。プラグインを買えないなら…自作してしまいましょう！！というわけで過去何回かプラグイン自作記事を書いてきました。JUCEを使うと簡単にプラグインを作ることができます。

もちろんEQも簡単に作れます。JUCEにはめっちゃ便利なクラスがあるので、これにサンプルを渡すだけでEQをかけれます。

でも、折角プラグインを自作するならこだわりたいですよね…？肝心要の信号処理部分は自分で1から練り上げたものを使いたいですよね…？

というわけでこの記事では、デジタル信号処理の基礎と、IIR型のEQの原理と実装方法を説明します！ライブラリに頼らずIIRなEQを自作プラグインに組み込めるようになりましょう！

目次

• DTMerのためのデジタル信号処理 ～IIRなEQを作ろう～
  • はじめに
  • 目次
  • Chapter 1. 前提知識
    • オイラーの公式
      • 直感的な証明
      • オイラーの公式の嬉しさ
    • フーリエ変換
      • フーリエ級数
    • 離散フーリエ変換
    • ラプラス変換
    • Z変換
    • 離散時間システムの表現
    • 離散時間システムの解析
    • 標本化定理(サンプリング定理)
  • Chapter 2. IIRデジタルフィルタの設計
    • 電気回路からの直接設計
      • 周波数特性の確認
      • コンピュータ実装のためのブロックダイアグラム
    • 双一次変換法
    • その他の著名な方法(Orfanidis)
    • その他の著名な方法(MZTi)
    • その他の著名な方法(Balance Mastering Magpha EQ / DDMF GrandEQ / Kush Audio q.632)
    • ちょっと乱暴な方法(双一次変換+オーバーサンプリング)
  • Chapter 3. IIRデジタルフィルタの実装
    • コンピュータの世界の数値表現（浮動小数点数）と演算誤差
      • 丸め誤差
      • 情報落ち
      • 桁落ち
    • ブロック法
      • 縦続接続
      • 並列接続
      • 帰還接続
    • 直接形I (Direct Form I)
    • 直接形II (Direct Form II)
    • 縦続形構成(Biquad Direct Form I)
    • 縦続形構成(Biquad Direct Form II)
    • 縦続形構成(Biquad Transposed Direct Form I)
    • 縦続形構成(Biquad Transposed Direct Form II)
    • 実験
  • Chapter 4. アナログフィルタの伝達関数
    • アナログ1次フィルタ
      • 低域通過フィルタ
      • 高域通過フィルタ
      • 低域シェルビングフィルタ
      • 高域シェルビングフィルタ
    • アナログ2次フィルタ
      • 低域通過フィルタ
      • 高域通過フィルタ
      • 帯域通過フィルタ
      • Notchフィルタ
      • 全域通過フィルタ
      • ピーキングフィルタ
      • 低域シェルビングフィルタ
      • 高域シェルビングフィルタ
    • アナログ高次フィルタ
  • Chapter 5. まとめ
  • おすすめ文献

Chapter 1. 前提知識

本文を理解するうえで必要な前提知識を浚う。

オイラーの公式

信号解析で頻出の「複素数乗」とは一体なんなのかを理解する必要がある。

$$e^{jx} = \cos x + j \sin x$$

直感的な証明

$e^x$をマクローリン展開すると、

$$e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!}$$

実は$x$に複素数を代入することができるので、

$$\begin{aligned} e^{jx} &= \sum_{n=0}^{\infty} \frac{(jx)^n}{n!} = \sum_{n=0}^{\infty} \frac{j^n}{n!}x^n \\ &= \sum_{k=0}^{\infty} \frac{j^{2k}}{(2k)!}x^{2k} + \sum_{k=0}^{\infty} \frac{j^{2k+1}}{(2k+1)!}x^{2k+1} \\ &= \sum_{k=0}^{\infty} \frac{(-1)^{k}}{(2k)!}x^{2k} + j \sum_{k=0}^{\infty} \frac{(-1)^{k}}{(2k+1)!}x^{2k+1} \\ \end{aligned}$$

ちなみに、

$$\begin{aligned} \cos x &= \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n)!}x^{2n} \\ \sin x &= \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)!}x^{2n+1} \end{aligned}$$

なので、

$$e^{jx} = \cos x + j \sin x$$

を得る。

オイラーの公式の嬉しさ

オイラーの公式の嬉しさを説明する為に、複素数平面の嬉しさを説明する。 複素数平面とは、横軸に実軸(Reと書くことが多い)、縦軸に虚軸(Imと書くことが多い)をとった平面である。通常の二次元平面では点Aの座標を$(1, 2)$と表現するが、複素数平面ではこの点を$1+2j$と表現する。ここで$j$とは虚数単位であり、$j^2=-1$であるような$j$のことを指す。

従来であれば1つの点の位置を表現するのに2つの値が必要だったが、複素数平面を用いることで、1つの点の位置を表現するのに1つの値で十分になる。

続いて、オイラーの公式の両辺に$A$を掛けた$Ae^{jx} = A \cos x + j A \sin x$を複素数平面上に描画するとこうなる。

複素数平面なら1つの点の位置を表現するのに1つの値で十分になるという意味が明白になったと思う。

フーリエ変換

下記の動画が分かりやすい。信号を単位円上に巻きつけて、それの重心を考えるというアイデアである。

https://www.youtube.com/watch?v=fGos3wrKeHY

この結果、連続的な信号$f(t)$に対して、下記のようなフーリエ変換の式を得る。

$$\mathcal{F}[f(t)](\omega) = \int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt$$

また、関数$f(t)$をフーリエ変換した物を$\mathcal{F}[f(t)]$や$\mathcal{F}[f(t)](\omega)$や$F(\omega)$と書くこともある。

$\omega$は角周波数であり、1秒間に何ラジアン進むかを意味する。1回転は$2 \pi \, [\mathit{\mathrm{rad}}]$なので、周期$T=\frac{2 \pi}{\omega}$であり、周波数$f = \frac{1}{T} = \frac{\omega}{2 \pi}\, [\mathrm{Hz}]$である。

性質	数式	説明
線形性	$\mathcal{F}[a f(t) + b g(t)] = a F(\omega) + b G(\omega)$	フーリエ変換は線形演算を保つ
時間シフト	$\mathcal{F}[f(t - t_0)] = e^{-j\omega t_0} F(\omega)$	関数の遅延は周波数領域で位相変化を引き起こす
周波数シフト	$\mathcal{F}[e^{j \omega_0 t} f(t)] = F(\omega - \omega_0)$	時間領域での指数関数の乗算は周波数領域での平行移動に対応
時間スケーリング	$\mathcal{F}[f(a t)] = \frac{1}{\lvert a \rvert} F ( \frac{\omega}{a} )$	時間領域でのスケーリングは周波数領域で逆スケーリングと振幅の変化を引き起こす
時間反転	$\mathcal{F}[f(-t)] = F(-\omega)$	時間反転は周波数領域での反転に対応する
畳み込み	$\mathcal{F}[f(t) * g(t)] = F(\omega) G(\omega)$	時間領域での畳み込みは周波数領域での積に対応する
時間領域での積	$\mathcal{F}[f(t) g(t)] = \frac{1}{2\pi} (F * G)(\omega)$	時間領域での積は周波数領域での畳み込みに対応する
積分定理	$\mathcal{F} \left[ \int_{-\infty}^{t} f(\tau) d\tau \right] = \frac{1}{j\omega}F(\omega) + \pi F(0) \delta(\omega)$	時間領域での積分は周波数領域での分数形式に変換される
微分定理	$\mathcal{F}[f'(t)] = j\omega F(\omega)$	時間領域の微分は周波数領域での乗算に対応する
フーリエ変換のフーリエ変換	$\mathcal{F}[\mathcal{F}[f(t)]] = 2\pi f(-t)$	時間反転と振幅の変化を伴う

フーリエ級数

前述のフーリエ変換を得るために、本来はフーリエ級数→フーリエ積分→フーリエ変換と導出していくのが普通だが、本記事ではフーリエ変換は道具の1つなので、天下り的にフーリエ変換と級数を紹介する。フーリエ級数は、周期信号を三角関数の和に展開する。

区間$-\frac{T}{2} \leq t \leq \frac{T}{2}$で定義された信号$f(t)$が区間外で周期$T$の周期信号であるなら、その信号は、

$$\begin{aligned} f(t) &= \sum_{n=-\infty}^{\infty} C_n e^{jn\frac{2 \pi}{T}t} \\ C_n&=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) e^{-jn\frac{2 \pi}{T}t} \end{aligned}$$

というフーリエ級数に展開できる。

離散フーリエ変換

前に示した動画のように、フーリエ変換のことを、ある信号を単位円上に巻きつけたその重心と考えるなら、離散的な信号を単位円上に巻きつけた時の重心は

$$\sum_{n=-\infty}^{\infty} f(nT)e^{-j \omega nT}$$

と書ける。これが離散フーリエ変換の式である。

ラプラス変換

フーリエ変換には便利な性質がある。例えば、微分記号を外すことができる。

$$\mathcal{F}[f^{(n)}](\omega) = (j \omega)^{n} \mathcal{F}[f](\omega)$$

現実世界の解析では微分方程式を解くことが多く、微分を外せるというフーリエ変換の性質は非常に便利である。

しかし、フーリエ変換の値が存在するためには、$f(t)$が$(- \infty, \infty)$で絶対可積分であること。つまり

$$\int_{-\infty}^{\infty} |f(t)| dt < \infty$$

であることが必要である。これだと$f(t)$の制約が強く、少し使い勝手が悪い。

そこで$f(t)$を $t \rightarrow \infty$で有界にすべく、$t$を$(0, \infty)$とし、新たなパラメータ$\sigma \ (\sigma > 0)$を用いて、$e^{-\sigma t}f(t)$と考える。これでフーリエ変換の式を書き直すと、

$$\begin{aligned} \mathcal{F}[e^{-\sigma t}f(t)](\omega) &= \int_{0}^{\infty} e^{-\sigma t}f(t) e^{-j\omega t} dt \\ &= \int_{0}^{\infty} f(t) e^{- (\sigma + j\omega) t} dt \end{aligned}$$

ここで、$\sigma + j\omega$は単なる複素数であるから$s$と置くことで、ラプラス変換の式

$$\mathcal{L}[f(t)](s) = \int_{0}^{\infty} f(t) e^{- st} dt$$

を得る。 関数$f(t)$をラプラス変換した物を$\mathcal{L}[f(t)]$や$\mathcal{L}[f(t)](s)$や$F(s)$と書くこともある。

性質	数式	説明
線形性	$\mathcal{L}[a f(t) + b g(t)] = a F(s) + b G(s)$	ラプラス変換は線形演算を保つ
時間シフト	$\mathcal{L}[f(t - t_0) u(t - t_0)] = e^{-st_0} F(s)$	関数の遅延が指数関数の乗算に変換
指数関数乗算の性質	$\mathcal{L}[e^{s_0t} f(t)] = F(s - s_0)$	時間領域での指数関数の乗算は変数の平行移動に対応
畳み込みのラプラス変換	$\mathcal{L}[f(t) * g(t)]= F(s) G(s)$	畳み込みはラプラス変換での積に対応
時間領域での積	$\mathcal{L}[f(t) g(t)] = \frac{1}{2\pi j} \int_{c - j\infty}^{c + j\infty} F(\sigma) G(s - \sigma) d\sigma$	積はラプラス変換での複素畳み込みに対応
積分のラプラス変換	$\mathcal{L} [\int_0^t f(\tau) d\tau] = \frac{F(s)}{s}$	時間領域での積分はラプラス変換での除算に対応
微分のラプラス変換	$\mathcal{L} [f'(t)] = s F(s) - f(0)$	関数の微分がラプラス変換では乗算と初期値に対応
2階微分のラプラス変換	$\mathcal{L}[f''(t)] = s^2 F(s) - s f(0) - f'(0)$	高階微分も同様に変換できる

Z変換

離散版のラプラス変換のような物である。離散フーリエ変換の式に対して、ラプラス変換の導出と同じように求める。$z$は複素数である。

$$\begin{aligned} \sum_{n=0}^{\infty} e^{-\sigma nT}f(nT)e^{-j \omega nT} &= \sum_{k=0}^{\infty} f(nT)e^{- (\sigma + j \omega) nT} \\ &= \sum_{n=0}^{\infty} f(nT)z^{-n} \end{aligned}$$

関数$f(t)$をZ変換した物を$\mathcal{Z}[f(t)]$や$\mathcal{Z}[f(t)](z)$や$F(z)$と書くこともある。

性質	数式	説明
線形性	$\mathcal{Z}[a f(t) + b g(t)] = a F(z) + b G(z)$	Z変換は線形演算を保つ
時間シフト	$\mathcal{Z}[f[n - k]] = z^{-k} F(z)$	時間領域のシフトは $z$ の乗算に対応
指数関数乗算の性質	$\mathcal{Z}[\alpha^n f(t)] = F(z/\alpha)$	時間領域での指数関数の乗算は変数のスケーリングに対応
畳み込みのZ変換	$\mathcal{Z}[f(t) * g(t)] = F(z) G(z)$	畳み込みはZ変換での積に対応
和のZ変換	$\mathcal{Z} \left[ \sum_{k=0}^{n} f[k] \right] = \frac{F(z)}{1 - z^{-1}}$	累積和はZ変換での除算に対応
微分のZ変換（離散時間）	$\mathcal{Z} [\Delta f(t)] = (z - 1) F(z)$	離散差分はZ変換では乗算に対応
2階差分のZ変換	$\mathcal{Z} [\Delta^2 f(t)] = (z - 1)^2 F(z)$	高階差分も同様に変換できる

離散時間システムの表現

デジタルフィルタは離散時間システムである。離散時間システムへの入力信号と出力信号は整数$n$とサンプリング周期$T$及び連続時間信号$f(t)$を用いて$f(nT)$と表現される。離散時間信号を、サンプリング周期を考えずに値の配列として考えて、$f[n]$と表現することもある。

離散時間信号を表現する重要なツールとして単位インパルス信号が挙げられる。

$$\delta (nT)= \begin{cases} 1, & n=0 \\ 0, & n \neq 0 \end{cases}$$

これは次のような波形を持つ。

インパルス信号の頂点は並行移動することができ、整数$k$を用いて

$$\delta (nT-kT)= \begin{cases} 1, & n=k \\ 0, & n \neq k \end{cases}$$

と書ける。これを用いれば任意の信号$x(nT)$は

$$x(nT) = \sum_{k=-\infty}^{\infty} x(kT) \delta (nT - kT)$$

と書ける。これはかなり当たり前な式だが、離散時間システムを表現するうえで役立つ。

インパルス信号を用いて離散時間システムを表現することを考えよう。離散時間システムは、離散時間信号$x(nT)$を入力して処理を行い離散時間信号$y(nT)$を得ることと言えるため、離散時間システムを$f$という関数と見れば

$$y(nT) = f[x(nT)]$$

と書ける。図で表すなら下記の通りである。

システムにインパルス信号を入力したときの出力、インパルス応答を$h(nT)=f[\delta(nT)]$とする。

前に述べたように$x(nT) = \sum_{k=-\infty}^{\infty} x(kT) \delta (nT - kT)$であるから、

$$y(nT) = f[x(nT)] = f \left[ \sum_{k=-\infty}^{\infty} x(kT) \delta (nT - kT) \right]$$

と書ける。離散時間システムは線形性(EQ掛けてから信号を足しても、信号を足してからEQを掛けても結果は同じ)を持つから、

$$f[a_1x_1(nT)+a_2x_2(nT)] = a_1f[x_1(nT)]+a_2f[x_2(nT)]$$

が成り立つので、

$$\begin{aligned} y(nT) = f[x(nT)] &= f \left[ \sum_{k=-\infty}^{\infty} x(kT) \delta (nT - kT) \right] \\ &= \sum_{k=-\infty}^{\infty} f \left[x(kT) \delta (nT - kT) \right] \\ &= \sum_{k=-\infty}^{\infty} x(kT) f \left[ \delta (nT - kT) \right] \end{aligned}$$

ここで、離散時間システムは時不変性(いつ入力しても結果は同じ)を持つので、インパルス応答の式から$h(nT-kT)=f[\delta(nT-kT)]$であるから、

$$\begin{aligned} y(nT) = f[x(nT)] &= \sum_{k=-\infty}^{\infty} x(kT) f \left[ \delta (nT - kT) \right] \\ &= \sum_{k=-\infty}^{\infty} x(kT) h(nT-kT) \end{aligned}$$

この式の右辺の形を畳み込み和（コンボリューション）と呼ぶ。つまり、任意の入力信号に対する出力信号は、システムのインパルス応答を入力信号に畳み込むと得られるということが分かる。

また、当たり前な式だと紹介した$x(nT) = \sum_{k=-\infty}^{\infty} x(kT) \delta (nT - kT)$は、実は畳み込み和の計算であり、インパルス信号$\delta(nT)$は畳み込み和における単位元ということも分かる。

システムのインパルス応答が有限長であればそのシステムはFIRシステム(Finite Impulse Response)であり、無限長であればIIRシステム(Infinite Impulse Response)である。 現実のシステムを考える際、FIRシステムであれば有限長のインパルス応答を畳み込めば良いので実現可能である。ただし、無限長のインパルス応答を畳み込むのは不可能であるためIIRシステムは実現できないように見える。しかし、差分方程式を使うことでIIRシステムは上手いこと実現可能である。

具体的に見てみる。インパルス応答が下記式で表せるようなIIRシステムを考える。

$$h(nT) = \begin{cases} \alpha ^n, & 0 \leq n \\ 0, & \text{otherwise} \end{cases}$$

インパルス応答$h(nT)$は$n<0$では常にゼロである。もし$h(nT)$が$n<0$でゼロでないなら、入力のインパルスが来る前に何等かの応答が表れていることになる。未来の入力を先取りして応答するシステムは現実では存在しない。ちなみに、インパルス応答$h(nT)$が$n<0$で常にゼロであるシステムを因果的なシステムと呼ぶ。現実に存在する全てのシステムは因果的である。このインパルス応答を用いれば入力信号$x(nT)$に対する出力信号$y(nT)$は畳み込みを利用して

$$\begin{aligned} y(nT) = f[x(nT)] &= \sum_{k=-\infty}^{\infty} x(kT) h(nT-kT) \\ &= \sum_{k=-\infty}^{\infty} h(kT) x(nT-kT) \quad \scriptsize{※畳み込みは可換}\\ &= \sum_{k=0}^{\infty} \alpha^k x(nT-kT) \end{aligned}$$

具体的な$n$と$k$を書き出してみると、

$$\begin{aligned} y(0 \times T) &= \sum_{k=0}^{\infty} \alpha^k x(-kT) \\ &= \alpha^0 x(0) + \alpha^1 x(-T) + \alpha^2 x(-2T) + \alpha^3 x(-3T) + \alpha^4 x(-4T) + \cdots \\ y(1 \times T) &= \sum_{k=0}^{\infty} \alpha^k x(T-kT) \\ &= \alpha^0 x(T) + \alpha^1 x(0) + \alpha^2 x(-T) + \alpha^3 x(-2T) + \alpha^4 x(-3T) + \cdots \\ &= \alpha^0 x(T) + \alpha ( \alpha^0 x(0) + \alpha^1 x(-T) + \alpha^2 x(-2T) + \alpha^3 x(-3T) + \cdots )\\ &= \alpha^0 x(T) + \alpha y(0 \times T) \\ y(2 \times T) &= \sum_{k=0}^{\infty} \alpha^k x(2T-kT) \\ &= \alpha^0 x(2T) + \alpha^1 x(T) + \alpha^2 x(0) + \alpha^3 x(-T) + \alpha^4 x(-2T) + \cdots \\ &= \alpha^0 x(2T) + \alpha (\alpha^0 x(T) + \alpha^1 x(0) + \alpha^2 x(-T) + \alpha^3 x(-2T) + \cdots ) \\ &= \alpha^0 x(2T) + \alpha y (1 \times T) \end{aligned}$$

段々法則も分かって来たので、

$$\begin{aligned} y(nT) = f[x(nT)] &= \sum_{k=0}^{\infty} \alpha^k x(nT-kT) \\ &= x(nT) + \sum_{k=1}^{\infty} \alpha^{k} x(nT-kT) \\ &= x(nT) + \sum_{k=0}^{\infty} \alpha^{k+1} x(nT-(k+1)T) \\ &= x(nT) + \alpha \sum_{k=0}^{\infty} \alpha^{k} x((n-1)T-kT) \\ &= x(nT) + \alpha y((n-1)T) \end{aligned}$$

以上より、前の結果を定数倍して現在の信号に加算すれば実現可能ということが分かった。

これを一般化すれば、離散時間システムの入出力関係の式は

$$\begin{aligned} y(nT) &= a_0x(nT) + a_1x(nT-T) + \cdots + a_Mx(nT-MT) \\ & \qquad -(b_0y(nT) + b_1y(nT-T) + \cdots + b_My(nT-NT)) \\ &= \sum_{k=0}^{M} a_k x(nT-kT) - \sum_{k=1}^{N}b_ky(nT-kT) \\ \end{aligned}$$

となる。これを差分方程式と呼ぶ。$M$と$N$のうち大きい方をシステムの次数と呼ぶ。

離散時間システムの解析

離散時間システムの解析にはZ変換を用いる。これは離散時間におけるラプラス変換のような物である。電気回路などの連続時間システムの解析にラプラス変換を用いることと対応づけると良い。

直前の項で、離散時間システムの入出力関係は差分方程式で表せると述べた。この差分方程式の両辺をZ変換してみる。

$$\begin{aligned} Y(z) &= \sum_{k=0}^{M} a_k z^{-k} X(Z) - \sum_{k=1}^{N} b_k z^{-k}Y(z) \\ \frac{Y(z)}{X(z)} &= \frac{\sum_{k=0}^{M} a_k z^{-k}}{1+\sum_{k=1}^{N} b_k z^{-k}} \end{aligned}$$

入出力の信号のZ変換の関係式$\frac{Y(z)}{X(z)}$を伝達関数$H(z)$と言う。Z変換では$z=e^{j\omega T}$を代入すれば離散フーリエ変換の結果を得ることができるため、この伝達関数に$z=e^{j\omega T}$を代入すればこのシステムの周波数特性$H(\omega)$を得ることができる。

周波数特性$H(\omega)$の振幅をとった$|H(\omega)|$を周波数振幅特性、あるいは単に振幅特性と言い、偏角をとった$\arg(H(\omega))$を周波数位相特性、あるいは単に位相特性と言う。

プラグインでよく見るスペクトラムアナライザは振幅特性の2乗に対数をとったもの、つまり、

$$\begin{aligned} \log_{10}(|H(\omega)|^2) &= 2 \log_{10}(|H(\omega)|)\ [\mathrm{B}] \quad \footnotesize{※単位はベル} \\ &= 20 \log_{10}(|H(\omega)|)\ [\mathrm{dB}] \quad \footnotesize{※単位はデシベル} \end{aligned}$$

の式で得られたものの内、正の周波数を表示したものである。

標本化定理(サンプリング定理)

アナログ世界の信号は連続信号であることが多い。これをデジタル世界で処理する為には、連続信号を離散信号に変換する必要がある。時間方向に離散化することを「標本化」や「サンプリング」と呼ぶ。振幅方向に離散化することを「量子化」と呼ぶ。

連続信号$f(t)$を周期$T$でサンプリングした信号の周波数特性を考えてみる。単純に考えれば連続信号$f(t)$を周期$T$でサンプリングした信号は

$$f_{\mathit{sampled}}(t) = \begin{cases} f(nT), & t = nT\quad (n は整数) \\ 0, & \text{otherwise} \end{cases}$$

と書けるため、フーリエ変換の式は下記のようになりそう。

$$\mathcal{F}=\int_{-\infty}^{\infty} f_{\mathit{sampled}}(t) e^{-i\omega t} dt$$

しかし、$f_{\mathit{sampled}}(t)$は面積がゼロであるため、上式の結果はゼロになる。これは困るので、$f_{\mathit{sampled}}(t)$の定義を考え直す必要がある。（ちなみに、これを短冊の面積と考えて計算すると、離散フーリエ変換の式が出てくる。それは今回欲しい物ではない。）

そこでディラックのデルタ関数を用意する。この関数は任意の連続な関数$g(x)$を用いて下記のように定義される。

$$\int_{-\infty}^{\infty} g(x) \delta(x) \, dx = g(0)$$

ディラックのデルタ関数は「関数(function)」という名前がついているが正しくは「超関数(distribution)」であり、関数ではない。なので$\delta(x)=〇〇$という形では書けない。ただ、なんとなくのイメージとして表現はできる。例えば$g(x)=1$の時、ディラックのデルタ関数の定義式から

$$\int_{-\infty}^{\infty} \delta(x) \, dx = 1$$

を得る。つまり、ディラックのデルタ関数の面積は1である。また定義式から、$g(x) \delta(x)$の演算を通して、どんな関数に対しても$x=0$の時の値$g(0)$しか拾わず、$x=0$以外の値は打ち消している様子から

$$\delta(x) = \begin{cases} \infty, & x = 0 \\ 0, & x \neq 0 \end{cases}$$

のようなイメージができる。なぜ$1$ではなく$\infty$なのか。前述の通りディラックのデルタ関数の面積は1であるので、面積が1の長方形を考える。ただし、この長方形の幅は0である。そのような長方形を考えると高さを$\infty$と表現するしかない。

注意しなければならないのが、これはあくまでもイメージということである。信号処理の分野では、$\delta (0)$の値を便宜上$1$と解釈してデルタ関数を扱うと直感的に分かりやすい場面が多い。

このディラックのデルタ関数を間隔$T$で並べたデルタ関数列は、

$$\delta_T(t) = \sum_{n=-\infty}^{\infty} \delta(t-nT)$$

と書ける。デルタ関数列と連続信号$f(t)$の積を取ることで、連続信号であることを保ったまま$f_{\mathit{sampled}}(t)$を得ることができる。

$$f_{\mathit{sampled}}(t)=f(t)\delta_T(t)$$

いよいよ$f_{\mathit{sampled}}(t)=f(t)\delta_T(t)$をフーリエ変換する。時間領域での積は周波数領域での畳み込みに対応するため、

$$\begin{aligned} \mathcal{F}[f(t)\delta_T(t)](\omega) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} \mathcal{F}[f(t)](u) \cdot \mathcal{F}[\delta_T(t)](\omega - u) \, du \end{aligned}$$

を計算すれば良い。まず、デルタ関数列のフーリエ変換から行う。 デルタ関数列は周期$T$の周期関数なのでフーリエ級数展開して、

$$\begin{aligned} \delta_T(t) &= \sum_{n=-\infty}^{\infty} \left( \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} \delta(t) e^{-jn\frac{2 \pi}{T}t} \right) e^{jn\frac{2 \pi}{T}t} \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} e^{jn\frac{2 \pi}{T}t} \end{aligned}$$

これをフーリエ変換すると

$$\begin{aligned} \mathcal{F}[\delta_T(t)](\omega) &= \int_{-\infty}^{\infty} \delta_T(t) e^{-j \omega t} \, dt \\ &= \int_{-\infty}^{\infty} \left( \frac{1}{T} \sum_{n=-\infty}^{\infty} e^{jn\frac{2 \pi}{T}t} \right) e^{-j \omega t} \, dt \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} e^{jn\frac{2 \pi}{T}t} e^{-j \omega t} \, dt \\ \end{aligned}$$

ここで、ディラックのデルタ関数のフーリエ変換$\mathcal{F}[\delta(t)](\omega)=1$と、フーリエ変換のフーリエ変換を用いれば$\mathcal{F}[1](\omega)=2 \pi \delta(\omega)$であり、更に周波数シフトを考えれば$\mathcal{F}[1 \cdot e^{j \omega_0 t}](\omega)=2 \pi \delta(\omega - \omega_0)$なので、

$$\begin{aligned} \mathcal{F}[\delta_T(t)](\omega) &= \int_{-\infty}^{\infty} \delta_T(t) e^{-j \omega t} \, dt \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} 2 \pi \delta(\omega - n\frac{2 \pi}{T}) \\ &= \frac{2 \pi}{T} \sum_{n=-\infty}^{\infty} \delta(\omega - n\frac{2 \pi}{T}) \\ \end{aligned}$$

これを使えば$f_{\mathit{sampled}}(t)=f(t)\delta_T(t)$のフーリエ変換は

$$\begin{aligned} \mathcal{F}[f(t)\delta_T(t)](\omega) &= \frac{1}{2 \pi} \int_{-\infty}^{\infty} \mathcal{F}[f(t)](u) \cdot \mathcal{F}[\delta_T(t)](\omega - u) \, du \\ &= \frac{1}{2 \pi} \int_{-\infty}^{\infty} \mathcal{F}[f(t)](u) \cdot \frac{2 \pi}{T} \sum_{n=-\infty}^{\infty} \delta(\omega - n\frac{2 \pi}{T} - u) \, du \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} \mathcal{F}[f(t)](u) \cdot \delta(\omega - n\frac{2 \pi}{T} - u) \, du \\ \end{aligned}$$

ディラックのデルタ関数は畳み込み演算に置ける単位元なので、

$$\begin{aligned} \mathcal{F}[f(t)\delta_T(t)](\omega) &= \frac{1}{T} \sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} \mathcal{F}[f(t)](u) \cdot \delta(\omega - n\frac{2 \pi}{T} - u) \, du \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} \mathcal{F}[f(t)](\omega - n\frac{2 \pi}{T}) \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} F(\omega - n\frac{2 \pi}{T}) \end{aligned}$$

この式を解釈すると、連続時間信号をサンプリングした離散信号のスペクトルは、連続時間信号のスペクトルを間隔$\frac{2 \pi}{T}$で並べた物と言える。間隔$\frac{2 \pi}{T}$とはサンプリング周波数である。これを図で描くと次のようになる。

連続時間信号のスペクトルの最大値が$\sigma$以下、つまり、連続時間信号が$\sigma$以下で帯域制限されている場合を図示している。これをサンプリングした信号のスペクトルは、連続時間信号のスペクトルをサンプリング周波数間隔で並べた物になる。もし$\sigma \geq \frac{1}{2} \times \frac{2 \pi}{T}$なら、サンプリング後のスペクトルが一部重複し信号が歪んでしまう。この歪みのことをエイリアシング歪みと呼ぶ。元信号をサンプリングする際にエイリアシング歪みが起こらないサンプリング周波数の下限のことをナイキストレートと呼び、あるサンプリング周波数でサンプリングしたときに表現できる周波数の上限のことをナイキスト周波数と呼ぶ。

折角なのでサンプリングした信号から連続時間信号を復元する処理も行ってみよう。ここではエイリアシング歪みが無かった時を考える。

連続時間信号をサンプリングした離散信号のスペクトルは、連続時間信号のスペクトルを間隔$\frac{2 \pi}{T}$で並べた物である。離散信号から連続信号を復元するには、離散信号のスペクトルのうち中心以外にあるものを打ち消すような処理を行えばいい。つまり、低域通過フィルタを適応させれば良い。

カットオフ周波数が$\omega_c$である理想低域通過フィルタの周波数特性は次のようにあらわせる。

$$H(\omega) = \begin{cases} 1, & |\omega| \leq \omega_C \\ 0, & \omega_C < |\omega| \end{cases}$$

理想低域通過フィルタを時間領域で表現するには$H(\omega)$を逆フーリエ変換すればよく、

$$h(t) = \frac{1}{2 \pi} \int_{-\omega_c}^{\omega_c} e^{j \omega t} = \frac{\omega_c}{\pi} \frac{\sin \omega_ct}{\omega_ct}$$

である。これはsinc関数である。

周波数領域にて、サンプリングした離散信号のスペクトルに$H(\omega)$を掛ければ、中心以外にあるものを打ち消すことができ、連続時間信号を復元できる。周波数領域に置ける積は時間領域に置ける畳み込みに対応する。つまり、離散信号にsinc関数を畳み込めば連続時間信号を復元できるのである。

Chapter 2. IIRデジタルフィルタの設計

いよいよIIRデジタルフィルタの設計に移る。実は、IIRデジタルフィルタを無から作り出すのは難しい。そこで、アナログフィルタの理論は昔からよく研究されていることから、アナログフィルタを何等かの形で変換してデジタルフィルタを得ることを考えよう。ここではいくつかの方法を紹介する。

電気回路からの直接設計

下記はRCローパスフィルタの回路である。これをデジタルフィルタとして設計し、プログラムで実装可能なブロックダイアグラムを得ることを目的とする。

周波数特性の確認

まずはこの回路がどのような周波数特性を持つのかを確認する。 この回路について成り立つ式は次の通り。

$$\begin{aligned} V_x &= V_r + V_c \\ V_y &= V_c \\ V_r &= IR \\ I &= \frac{d Q_c}{dt} \\ Q_c &= C V_c \end{aligned}$$

これを整理すると次のようになる。

$$\begin{aligned} V_x &= V_r + V_c \\ &= IR + V_y \\ &= R \frac{d Q_c}{dt} + V_y \\ &= R \frac{d C V_c}{dt} + V_y \\ &= R \frac{d C V_y}{dt} + V_y \\ &= R C \frac{d V_y}{dt} + V_y \end{aligned}$$

微分方程式となり少々厄介なので、両辺をラプラス変換して微分を消す。

$$\begin{aligned} \mathcal{L}[V_x](s) &= s R C \mathcal{L}[V_y](s) + \mathcal{L}[V_y](s) \\ &= (1+sRC) \mathcal{L}[V_y](s) \end{aligned}$$

この回路の周波数特性を知るためには、入出力信号の周波数特性の比$\frac{\mathcal{F}[V_y](\omega)}{\mathcal{F}[V_x](\omega)}$を求めればよい。ラプラス変換を用いると、$\mathcal{F}[V_y](\omega)=\mathcal{L}[V_y](j \omega)$かつ$\mathcal{F}[V_x](\omega)=\mathcal{L}[V_x](j \omega)$と書けるから、伝達関数$H(s)$は

$$H(s) = \frac{\mathcal{L}[V_y](s)}{\mathcal{L}[V_x](s)} = \frac{1}{1+sRC}$$

となり、周波数特性$H(j \omega)$は

$$H(j \omega) = \frac{\mathcal{L}[V_y](j \omega)}{\mathcal{L}[V_x](j \omega)} = \frac{1}{1+j \omega RC}$$

となる。カットオフ周波数のパラメータを$\omega_c = \frac{1}{RC}$のように定めれば、

$$H(j \omega) = \frac{1}{1+j \omega / \omega_c}$$

振幅特性と位相特性はそれぞれ次のように書ける。

$$\begin{aligned} |H(j \omega)| &= \left| \frac{1}{1+j \omega / \omega_c} \right| \\ \arg H(j \omega) &= \arg \frac{1}{1+j \omega / \omega_c} \end{aligned}$$

カットオフ周波数を1000Hzにするために、$\omega_c = \frac{1}{RC} = 2\pi \times 1000$ とした時の振幅特性と位相特性は下記のようになる。

コンピュータ実装のためのブロックダイアグラム

この回路をデジタルフィルタとしてプログラムで実装するために、信号の処理過程を示すブロックダイアグラムを得る。 まず、実際の信号がどのように処理されるのかを知るため、微分方程式を実際に解く。 伝達関数を書き換える。

$$\begin{aligned} (1+sRC) \mathcal{L}[V_y](s) &= \mathcal{L}[V_x](s)\\ sRC \mathcal{L}[V_y](s) &= \mathcal{L}[V_x](s) - \mathcal{L}[V_y](s) \\ \mathcal{L}[V_y](s) &= \frac{1}{s} \cdot \frac{1}{RC} \cdot (\mathcal{L}[V_x](s) - \mathcal{L}[V_y](s)) \end{aligned}$$

逆ラプラス変換を行う。

$$\begin{aligned} \mathcal{L}[V_y](s) &= \frac{1}{s} \cdot \frac{1}{RC} \cdot (\mathcal{L}[V_x](s) - \mathcal{L}[V_y](s)) \\ &= \frac{1}{RC} \cdot \left\{\frac{1}{s} \cdot \mathcal{L}[V_x](s) - \frac{1}{s} \cdot \mathcal{L}[V_y](s) \right\} \\ V_y &= \frac{1}{RC} \cdot \left\{ \int_0^{t} V_x(\tau) d\tau - \int_0^{t} V_y(\tau) d\tau \right\} \\ &= \frac{1}{RC} \cdot \left\{ \int_0^{t} V_x(\tau) - V_y(\tau) d\tau \right\} \end{aligned}$$

$V_x$と$V_y$の式を得た。しかし、$V_x$と$V_y$は連続的な値である。プログラムで実装するするためには、離散的な信号の処理を知る必要がある。 そこで$V_x$と$V_y$を離散的な値として見なす。すると、$\int_0^{t} V_x(\tau) - V_y(\tau) d\tau$は幅$T$の短冊の面積の合計と考えることができ、積分記号を積和記号に置き換えることができる。

$$V_y(kT) = \frac{1}{RC} \cdot T \left\{ \sum_{n=0}^{k} V_x(nT) - V_y(nT) \right\}$$

これのブロックダイアグラムを考えればゴールだが、積和記号をどのようにブロックダイアグラムにすれば良いのかすぐには分からない。

問題を簡単にし、例として$V_y(t) = \int_0^{t} V_x(\tau) d\tau$という連続信号を離散化してそれのダイアグラムを考えてみる。 離散化すると、

$$V_y(kT) = T \sum_{n=0}^{k} V_x(nT)$$

を得る。分かりやすいように$k$に具体的な値を入れて書き出すと、

$$\begin{aligned} V_y(0) &= T \sum_{n=0}^{0} V_x(nT) = TV_x(0) \\ V_y(T) &= T \sum_{n=0}^{1} V_x(nT) = TV_x(0) + TV_x(T) = V_y(0) + TV_x(T) \\ V_y(2T) &= T \sum_{n=0}^{2} V_x(nT) = TV_x(0) + TV_x(T) + TV_x(2T)= V_y(T) + TV_x(2T) \\ V_y(3T) &= T \sum_{n=0}^{3} V_x(nT) = TV_x(0) + TV_x(T) + TV_x(2T) + TV_x(3T)= V_y(2T) + TV_x(3T) \\ \cdots \\ V_y(kT) &= V_y\left((k-1)T\right) + TV_x(kT) \end{aligned}$$

これからブロックダイアグラムを書くと、 である。ここで$z^{-1}$とは信号を1サンプル遅らせるという意味である。

これを使えば、先の式のダイアグラムを書くことができる。

$$\begin{aligned} V_y(kT) &= \frac{1}{RC} \cdot T \left\{ \sum_{n=0}^{k} V_x(nT) - V_y(nT) \right\} \\ \end{aligned}$$

これをプログラムで実装すれば、デジタルフィルタを実現できる。 カットオフ周波数が1000Hzになるようにパラメータを設定して実装し、振幅特性と位相特性を図示すると下記のようになる。サンプリング周波数は48000Hzである。青色で示されるAnalog RC LPFがアナログフィルタの特性を示し、橙色のDigital Naive RC LPFが上記ダイアグラムを実装して作ったデジタルフィルタの特性である。ナイキスト周波数に近づくにつれてアナログフィルタの特性と乖離していることが分かる。

双一次変換法

デジタルフィルタをアナログフィルタ電気回路から直接設計する手法はあまり効率的ではない。出来れば、アナログフィルタの伝達関数を良い感じに変形すればデジタルフィルタを得られるようにしたい。その方法が双一次変換法である。

標本化定理の所で扱ったように、連続時間信号の周波数領域は $-\infty \sim \infty$ の範囲に広がるのに対し、離散時間信号の周波数領域は $-\frac{\pi}{T} \sim \frac{\pi}{T}$ に限られる。双一次変換法は等角写像を利用し、この $-\infty \sim \infty$ の範囲を $-\frac{\pi}{T} \sim \frac{\pi}{T}$ に圧縮する手法である。

アナログフィルタの伝達関数から周波数特性を得るには$s=j\omega$を代入すれば良い。このことは、$s$平面の虚軸上で伝達関数の値を評価することと同じである。一方、デジタルフィルタの伝達関数から周波数特性を得るには$z=e^{j\omega}$を代入すればよい。このことは、$z$平面の単位円上で伝達関数の値を評価することと同じである。つまり、$s$平面上の虚軸を$z$平面上の単位円に対応させる変換を見つければ良い。結論を述べると、そのような変換は$s=\frac{1-z^{-1}}{1+z^{-1}}$で実現できる。

双一次変換法は、アナログフィルタの解析の結果得た伝達関数$H(s)$に$s=\frac{1-z^{-1}}{1+z^{-1}}$を代入することで、デジタルフィルタの伝達関数を得る方法である。双一次変換法はこのように、非常に簡単な処理でデジタルフィルタを得ることができる。しかし、問題がある。 アナログフィルタの周波数特性は $s = j\Omega$ の代入により得られ、デジタルフィルタの周波数特性は $z = e^{j\omega T}$ の代入により得られる。双一次変換を通して得たデジタルフィルタの周波数特性は、元のアナログフィルタの伝達関数$H(s)$を用いて表すと、$H\left( \frac{1-e^{-j \omega T}}{1+e^{-j \omega T}} \right)$となる。元のアナログフィルタの周波数特性$H(j\Omega)$と、双一次変換後の周波数特性$H\left( \frac{1-e^{-j \omega T}}{1+e^{-j \omega T}} \right)$について、アナログ領域の$\Omega$とデジタル領域の$\omega$の関係式は、

$$\begin{aligned} j \Omega &= \frac{1-e^{-j \omega T}}{1+e^{-j \omega T}} = \tanh \left( j\frac{\omega T}{2} \right) = j \tan \frac{\omega T}{2} \\ \Omega &= \tan \frac{\omega T}{2} \end{aligned}$$

となる。つまり対応関係が非線形であり、周波数尺度が変わってしまう。この対応関係を多少まともにする方法として、双一次変換の際にアナログフィルタの伝達関数$H(s)$に$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$を代入する方法が挙げられる。

これにより、アナログ領域の$\Omega$とデジタル領域の$\omega$の関係式は

$$\begin{aligned} \Omega &= \frac{2}{T} \tan \frac{\omega T}{2} \end{aligned}$$

となる。この関係を、角周波数を周波数に変更して対数軸にプロットした物が下図である。サンプリング周波数は48000Hzとしている。このグラフから、$\Omega$と$\omega$は低周波数領域では理想的な関係$\Omega=\omega$に沿って遷移するが、ナイキスト周波数に近づくにつれて段々ずれていくことが分かる。

実際に双一次変換法を用いて、先に説明したRC低域通過フィルタの伝達関数からデジタルフィルタを得てみよう。 RC低域通過フィルタの伝達関数は次の通り。

$$H(s) = \frac{1}{1+sRC}$$

$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$を代入して、

$$\begin{aligned} H(z) &= \frac{1}{1+RC \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}} \\ &= \frac{T + Tz^{-1}}{T + Tz^{-1} + 2RC- 2RCz^{-1}} \\ &= \frac{T + Tz^{-1}}{(T + 2RC) + (T - 2RC)z^{-1}} \\ &= \frac{\frac{T}{T + 2RC} + \frac{T}{T + 2RC}z^{-1}}{1 + \frac{T - 2RC}{T + 2RC}z^{-1}} \\ \end{aligned}$$

離散時間信号の差分方程式が

$$\begin{aligned} y(nT) &= \sum_{k=0}^{M} a_k x(nT-kT) - \sum_{k=1}^{N}b_ky(nT-kT) \\ \end{aligned}$$

である時、そのZ変換が

$$\begin{aligned} \frac{Y(z)}{X(z)} &= \frac{\sum_{k=0}^{M} a_k z^{-k}}{1+\sum_{k=1}^{N} b_k z^{-k}} \end{aligned}$$

なので、係数$a_0, a_1, b_1$は

$$a_0 = \frac{T}{T + 2RC}, \quad a_1 = \frac{T}{T + 2RC}, \quad b_1 = \frac{T - 2RC}{T + 2RC}$$

である。

RC低域通過フィルタは1次のフィルタなので、改めて差分方程式を書くと、

$$\begin{aligned} y(nT) &= a_0 x(nT) + a_1 x(nT-T) - b_1 y(nT-T) \end{aligned}$$

である。

カットオフ周波数が1000Hzになるようにパラメータを設定して実装し、振幅特性と位相特性を図示すると下記のようになる。サンプリング周波数は48000Hzである。青色で示されるAnalog RC LPFがアナログフィルタの特性を示し、橙色のDigital RC LPFが上記を実装して作ったデジタルフィルタの特性である。

概ね形がアナログフィルタに近いもののナイキスト周波数に近づくにつれてずれ始めることが分かる。

更にカットオフ周波数を高くすると、指定した周波数より早く減衰が始まってしまう。そこで$\Omega = \frac{2}{T} \tan \frac{\omega T}{2}$の式を用いて事前に補正を掛ける（この処理をプリワーピングと呼ぶ）方法もある。例えばカットオフ周波数を12000Hzにしたいなら、$\Omega$が12000Hzの値になるよう全体をスケールする。下図の緑色のグラフがプリワーピングの結果である。

その他の著名な方法(Orfanidis)

双一次変換法は非常に便利だが、オーディオ信号処理においては割と致命的な欠陥を持つ。そこで非常に多くの研究が行われてきた。ここでは比較的有名な研究を紹介する。

1つ目が1997年のS.J.Orfanidisによる研究である。

S. J. Orfanidis, “Digital parametric equalizer design with prescribed Nyquist-frequency gain,” Journal of the Audio Engineering Society, vol. 45, no. 6, pp. 444-455, Jun. 1997.

ナイキスト周波数に近づくにつれてどんどん減衰して行ってしまうので、その分ゲインして補正してあげようというアイデアである。

Orfanidisの手法は2次のピークフィルタでのみ利用できる。詳しくは後述するが、高次のフィルタも演算誤差を抑える為に2次以下のフィルタの連結で表現することが普通であるため、2次フィルタのみ可という制約はそれほど甚大ではない。ベル型以外のフィルターは別の方法を適用する。

まず、ゲインが DC($f=0\,[\mathrm{Hz}]$) と無限大($f=\infty\,[\mathrm{Hz}]$)で$G_0$である 2 次アナログフィルタの伝達関数は、ピーク中心周波数$\Omega_0$と係数$A,\ B$を用いて、

$$H(s) = \frac{G_0 s^2 + B s + G_0 \Omega_0^2}{s^2 + A s + \Omega_0^2}$$

のように表せる。なお$A,\ B$は、周波数$\Omega_0$でのゲイン$G$、フィルタの幅$\Delta \Omega$、その幅を持つ時のゲイン$G_B$を用いて下記である。

$$A = \sqrt{\frac{G_B^2 - G_0^2}{G^2 - G_B^2}} \Delta\Omega, \quad B = G A$$

ここから、ナイキスト周波数が$G_1$である時、$f=\infty\,[\mathrm{Hz}]$の地点のゲインも$G_1$であるような新たなフィルタ

$$H(s) = \frac{G_1 s^2 + B s + G_0 W^2}{s^2 + A s + W^2}$$

の係数を求める。詳細な導出は元論文を辿って欲しい。結論を述べると、

$$W^2 = \sqrt{\frac{G^2 - G_1^2}{G^2 - G_0^2}} \Omega_0^2, \quad A = \sqrt{\frac{C + D}{|G^2 - G_B^2|}}, \quad B = \sqrt{\frac{G^2 C + G_B^2 D}{|G^2 - G_B^2|}} \\ C = (\Delta\Omega)^2 |G_B^2 - G_1^2| - 2W^2 \left( |G_B^2 - G_0 G_1| - \sqrt{(G_B^2 - G_0^2)(G_B^2 - G_1^2)} \right) \\ D = 2W^2 \left( |G^2 - G_0 G_1| - \sqrt{(G^2 - G_0^2)(G^2 - G_1^2)} \right)$$

である。これに双一次変換を適応して、以下を得る。

$$H(z) = \frac{\left( \frac{G_1+G_0W^2+B}{1+W^2+A} \right) -2 \left( \frac{G_1-G_0W^2}{1+W^2+A} \right)z^{-1} + \left( \frac{G_1+G_0W^2-B}{1+W^2+A} \right)z^{-2}}{1-2\left( \frac{1-W^2}{1+W^2+A} \right)z^{-1} + \left( \frac{1+W^2-A}{1+W^2+A} \right)z^{-2}}$$

なお双一次変換は$s=\frac{1-z^{-1}}{1+z^{-1}}$を使っており、デジタル領域の周波数は正規化済みであるから、中心周波数$f_0\, \mathrm{[Hz]}$とサンプリング周波数$f_s\, \mathrm{[Hz]}$の時、$\omega_0 = \frac{2 \pi f_0}{f_s}, \Delta \omega = \frac{2 \pi \Delta f}{f_s}$として、

$$\Omega_0 = \tan \left( \frac{\omega_0}{2} \right), \quad \Delta \Omega = \left( 1 + \sqrt{\frac{G_B^2 - G_0^2}{G_B^2 - G_1^2}} \sqrt{\frac{G^2 - G_1^2}{G^2 - G_0^2}} \Omega_0^2 \right) \tan \left( \frac{\Delta \omega}{2} \right)$$

とする。

実際にこの手法を用いてみよう。アナログ2次ピークフィルタの一例として、次のような伝達関数を持つものを考える。$R$は抵抗値だがゲインのような働きをする。

$$H(s) = \frac{\Omega_0^2-s^2}{s^2+2R\Omega_0s+\Omega_0^2}$$

まず$W$を求める為には$G$及び$G_0,\ G_1$を求める必要がある。$G$は中心周波数におけるゲインであるから、$s=j\Omega_0$を代入して

$$\begin{aligned} G=|H(j\Omega_0)| &=\left| \frac{\Omega_0^2-(j\Omega_0)^2}{(j\Omega_0)^2+2R\Omega_0(j\Omega_0)+\Omega_0^2} \right| \\ &= \frac{1}{R} \end{aligned}$$

$G_0$はDC成分のゲインである。DC成分のゲインは$|H(0)|$なので、

$$G_0 = |H(0)| = 1$$

$G_1$はナイキスト周波数におけるゲインである。中心周波数$f_0\, \mathrm{[Hz]}$とサンプリング周波数$f_s\, \mathrm{[Hz]}$の時、$\omega_0 = \frac{2 \pi f_0}{f_s}, \Delta \omega = \frac{2 \pi \Delta f}{f_s}$として、

$$G_1^2 = \frac{ G_0^2 (\omega_0^2 - \pi^2)^2 + G^2 \pi^2 (\Delta \omega)^2 \left( G_B^2 - G_0^2 \right) / \left( G^2 - G_B^2 \right) }{ (\omega_0^2 - \pi^2)^2 + \pi^2 (\Delta \omega)^2 \left( G_B^2 - G_0^2 \right) / \left( G^2 - G_B^2 \right) }$$

続いて、$A$を求める為には$G_B$および$\Delta\Omega$が必要である。どちらかを決めればもう片方は自動で決まるので、今回は$G_B=\sqrt{G_0G}$とする。

サンプリング周波数を48000Hz、中心周波数を10000Hz、$R=0.2$にした時の振幅特性と位相特性を下図に示す。 青色がアナログフィルタの特性を示し、橙色が双一次変換を行って作ったデジタルフィルタ、緑色が双一次変換とプレワーピングを行って作ったデジタルフィルタ、赤色が上記のOrfanidisの手法を実装して作ったデジタルフィルタの特性である。 振幅特性のグラフを見ると、橙・緑・赤の中で赤色のグラフが最も青色に近いことが分かる。 ただし、位相特性のグラフを見ると赤色のグラフは他と傾向が違うことが分かる。これは、Orfanidisの方法は実質、元々のアナログフィルタにハイシェルフを適応させたものと同じであると見なせるからである。

その他の著名な方法(MZTi)

D. W. Gunness and O. S. Chauhan, “Optimizing the magnitude response of matched z-transform filters (“mzti”) for loudspeaker equalization”, journal of the audio engineering society, no. 2, Sep. 2007.

Orfanidisの手法は、通常の双一次変換と比較して大きな改善を達成している。しかし、いくつか問題がある。まず、位相特性がアナログフィルタと大きく異なっている。そして、ナイキスト周波数付近の振幅特性がアナログフィルタと一致するようにデジタルフィルタを設計したが、通常、ナイキスト周波数付近の音は聴覚できない。また、リサンプリングの際に用いるアンチエイリアシングフィルタはナイキスト周波数付近のゲインをゼロに向かわせるため、そこでの特性を一致させることには意味がない。更に、聴覚可能な周波数帯域において、アナログフィルタとの特性の乖離が生じており、問題である。MZTiは、この課題意識から生まれた手法である。

前々節では、アナログフィルタをデジタルフィルタに変換する方法として双一次変換を紹介した。前節のOrfanidisでも双一次変換を使っている。この変換は、$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$を代入することで実現される。つまり、何らかの方法で$s$を$z$に変換している。実は、この変換にはさまざまな選択肢が存在する。双一次変換は、周波数特性をなるべく保ったまま$s$を$z$に変換することを目的とした手法である。これから紹介するMZTiは、双一次変換ではなくMatched Z-transform Method (MZT)を用いて$s$を$z$に変換し、その後、多少の補正を加える(improved)ことで改善を図る手法である。

Matched Z-transform Method (MZT)は、アナログシステムの全ての極と零点を保ったままデジタルシステムに変換する方法である。具体的には、伝達関数が

$$\begin{aligned} H(s) &= \frac{a_Ms^M+\cdots+a_1s+a_0}{b_Ns^N+\cdots+b_1s+b_0} \\ &= \frac{a_M}{b_N} \frac{\prod_{i=1}^M(s-\xi_i)}{\prod_{i=1}^N(s-p_i)} \end{aligned}$$

と書かれている時、零点$\xi_i$と極$p_i$を、サンプリング周期を$T$として$z=e^{sT}$という対応関係を使って、

$$H(z)=g\frac{a_M}{b_N}\frac{\prod_{i=1}^M(z-e^{\xi_iT})}{\prod_{i=1}^N(z-e^{p_iT})}$$

とする。このように変形することで、元のアナログシステムが安定であれば、変換後のデジタルシステムも安定となる。なお、$g$は任意のスケール値であり、例えば$H(s)$と$H(z)$のDC成分のゲインを一致させるように設定する。

なぜMZTiでは双一次変換ではなくMZTを用いるのか。それは、アナログフィルタとデジタルフィルタの特性の差分を考えたとき、MZTの方が差分の形が単純であり、補正が容易だからである。実際に見てみよう。Orfanidisでも扱ったピークフィルタを使う。中心周波数は$\omega_c$である。

$$\begin{aligned} H(s) &= \frac{\omega_c^2-s^2}{s^2+2R\omega_cs+\omega_c^2} \\ &= -\frac{(s-\omega_c)(s+\omega_c)}{(s+R\omega_c-\omega_c\sqrt{R^2-1})(s+R\omega_c+\omega_c\sqrt{R^2-1})} \end{aligned}$$

アナログシステムにおける零点は$-\omega_c$と$\omega_c$で、極点は$-R\omega_c \plusmn \omega_c\sqrt{R^2-1}$。$z=e^{sT}$という対応関係を使って、$z$平面上での零点と極点を考えれば、

$$\begin{aligned} H(z) &= -g\frac{(z-e^{\omega_cT})(z-e^{-\omega_cT})}{(z-e^{(-R\omega_c-\omega_c\sqrt{R^2-1})T})(z-e^{(-R\omega_c+\omega_c\sqrt{R^2-1})T})}\\ &= g\frac{-z^2+(e^{\omega_cT}+e^{-\omega_cT})z-1}{z^2-(e^{(-R\omega_c+\omega_c\sqrt{R^2-1})T}+e^{(-R\omega_c-\omega_c\sqrt{R^2-1})T})z + e^{-2R\omega_cT}} \\ &= g\frac{-1+(e^{\omega_cT}+e^{-\omega_cT})z^{-1}-z^{-2}}{1 - (e^{(-R\omega_c+\omega_c\sqrt{R^2-1})T}+e^{(-R\omega_c-\omega_c\sqrt{R^2-1})T})z^{-1} + e^{-2 R\omega_c T }z^{-2} } \end{aligned}$$

ここで、$\cosh x = \frac{e^x+e^{-x}}{2}$なので

$$\begin{aligned} H(z) &= g\frac{-1+(e^{\omega_cT}+e^{-\omega_cT})z^{-1}-z^{-2}}{1 - (e^{(-R\omega_c+\omega_c\sqrt{R^2-1})T}+e^{(-R\omega_c-\omega_c\sqrt{R^2-1})T})z^{-1} + e^{-2 R\omega_c T }z^{-2}} \\ &= g\frac{-1+(2 \cosh(ω_cT))z^{-1}-z^{-2}}{1 - (2 e^{-R \omega_c T} \cosh( \omega_c T \sqrt{R^2 - 1}))z^{-1} + e^{-2 R\omega_c T }z^{-2}} \\ \end{aligned}$$

最後に$g$を求める。DC成分のゲインを基準にするため、$H(s)$に$s=0$を代入して$H(0)=1$を得て、$H(z)$に$z=1$を代入して

$$\begin{aligned} H(1)&=g\frac{-1+(2 \cosh(ω_cT))-1}{1 - (2 e^{-R \omega_c T} \cosh( \omega_c T \sqrt{R^2 - 1})) + e^{-2 R\omega_c T }}\\ &=g\frac{-2+2 \cosh(ω_cT)}{1 - (2 e^{-R \omega_c T} \cosh( \omega_c T \sqrt{R^2 - 1})) + e^{-2 R\omega_c T }}\\ &=1 \end{aligned}$$

より、

$$g = \frac{1 - (2 e^{-R \omega_c T} \cosh( \omega_c T \sqrt{R^2 - 1})) + e^{-2 R\omega_c T }}{-2+2 \cosh(ω_cT)}$$

サンプリング周波数を48000Hz、中心周波数を10000Hz、$R=0.2$にした時の振幅特性と位相特性を下図に示す。 青色がアナログフィルタの特性を示し、橙色が双一次変換とプリワーピングを行って作ったデジタルフィルタ、緑色がMZTによるデジタルフィルタである。MZTによるデジタルフィルタは、ナイキスト周波数付近で大きなずれはあるが、それ以外は双一次変換のものに比べて悪くなさそうである。

続いて、アナログフィルタの特性とデジタルフィルタの特性の差分をとって比較する。上段のグラフでは、双一次変換とプリワーピングで得たデジタルフィルタとアナログフィルタの振幅特性の差分を緑色のグラフで示している。下段のグラフでは、MZTで得たデジタルフィルタフィルタとアナログフィルタの差分を緑色のグラフで示している。デジタルフィルタに変換した後、何等かのフィルタを使って形状を補正する際、MZTの結果を補正したほうが簡単そうだと分かる。

補正フィルタの導出方法を示す。補正フィルタとしてFIRフィルタを用いる。元のアナログフィルタが2次である時、補正の対象となるMZTによるデジタルフィルタの周波数特性は

$$H_D(\omega) = g\frac{1+a_1e^{-j \omega T}+a_2e^{-2j \omega T}}{1+b_1e^{-j \omega T}+b_2e^{-2j \omega T}}$$

である。一般に、FIRフィルタの伝達関数は

$$H_{\mathit{FIR}}(z)=\sum_{k=0}^{M} \alpha_k z^{-k}$$

と表せるので、タップ数を決めて、いくつかの周波数でアナログフィルタの特性とMZT後のデジタルフィルタの特性が同じになるよう連立方程式をたてて、FIRフィルタの係数を求めれば良い。論文では、式を簡単にするため、ナイキスト周波数を等分するポイントで連立方程式を立式する方法を紹介している。

例えば、FIRフィルタの伝達関数が$H_{\mathit{FIR}}(z)=\alpha_0+\alpha_1z^{-1}+\alpha_2z^{-2}$の時、ナイキスト周波数を3分割した地点のデジタルフィルタの振幅特性は

$$\begin{aligned} |H_D(f_N / 3)| &= g \frac{\sqrt{1 + a_1 - a_2 + a_1^2 + a_1 a_2 + a_2^2}}{\sqrt{1 + b_1 - b_2 + b_1^2 + b_1 b_2 + b_2^2}} \\ |H_D(2f_N / 3)| &= g \frac{\sqrt{1 - a_1 - a_2 + a_1^2 - a_1 a_2 + a_2^2}}{\sqrt{1 - b_1 - b_2 + b_1^2 - b_1 b_2 + b_2^2}} \end{aligned}$$

であり、同地点の補正FIRフィルタの周波数特性を$H_1$と$H_2$とし、DC成分の周波数特性を$H_0$とすれば、

$$\begin{aligned} \alpha_0 &= H_0 - \alpha_1 - \alpha_2\\ \alpha_1 &= \frac{H_0 - \sqrt{H_0^2 - 2H_1^2 + 2H_2^2}}{2}\\ \alpha_2 &= \frac{3(H_0 - \alpha_1) - \sqrt{-3H_0^2 + 12H_1^2 - 6H_0 \alpha_1 - 3\alpha_1^2}}{6} \end{aligned}$$

この補正FIRフィルタを、MZT後のデジタルフィルタの後に掛けることで、DC成分・ナイキスト周波数の3分の1の地点・ナイキスト周波数の3分の2の地点の振幅特性がアナログフィルタの特性と一致するフィルタを得ることができる。

サンプリング周波数を48000Hz、中心周波数を10000Hz、$R=0.2$にした時の振幅特性と位相特性を下図に示す。 青色がアナログフィルタの特性を示し、橙色がMZTを行って作ったデジタルフィルタ、緑色がMZTで作ったフィルタにFIRで補正を行って作ったMZTiデジタルフィルタである。DC、ナイキスト周波数の3分の1の地点、3分の2の地点のゲインがアナログフィルタと一致するような補正フィルタを構成している為、ナイキスト周波数におけるゲインはアナログフィルタの値と異なるが、それ以外は一致していることが分かる。

その他の著名な方法(Balance Mastering Magpha EQ / DDMF GrandEQ / Kush Audio q.632)

J. Flynn and J. D. Reiss, “Improving the frequency response magnitude and phase of analogue-matched digital filters,” Journal of the Audio Engineering Society, May 2018.

MZTiはかなり優秀な手法である。MZTiで得たデジタルフィルタの特性はアナログフィルタにかなり近い。が、完璧に一致してはいない。 MZTiを改善しようと試みる時、補正FIRフィルタの算出方法を改善すれば良いのでは？と思うのはとても自然である。MZTiの補正フィルタは、ナイキスト周波数未満を等分した数点の振幅特性がアナログフィルタの振幅特性に一致するよう、連立方程式をたててFIRフィルタの係数を求めている。そうではなく、MZT後のデジタルフィルタの特性とアナログフィルタの特性の差を逆フーリエ変換してFIRフィルタの係数を求めた方が、よりアナログフィルタの特性を近似できるだろうというのは直感的に理解できる。

この改善案を採用したのが、これから紹介する手法である。

MZTiの項で述べたように、元のアナログフィルタの伝達関数が

$$\begin{aligned} H_A(s) &= \frac{a_Ms^M+\cdots+a_1s+a_0}{b_Ns^N+\cdots+b_1s+b_0} \\ &= \frac{a_M}{b_N} \frac{\prod_{i=1}^M(s-\xi_i)}{\prod_{i=1}^N(s-p_i)} \end{aligned}$$

と書かれている時、全ての零点$\xi_i$と極$p_i$を、サンプリング周期を$T$として$z=e^{sT}$という対応関係を使ってMZTしたデジタルフィルタの伝達関数は

$$H_{\mathit{MZT}}(z)=g\frac{a_M}{b_N}\frac{\prod_{i=1}^M(z-e^{\xi_iT})}{\prod_{i=1}^N(z-e^{p_iT})}$$

である。

アナログフィルタの周波数特性は$s=j\omega$を代入して、

$$\begin{aligned} H_A(j\omega) = \frac{a_M}{b_N} \frac{\prod_{i=1}^M(j\omega-\xi_i)}{\prod_{i=1}^N(j\omega-p_i)} \end{aligned}$$

デジタルフィルタの周波数特性は$z=e^{j\omega}$を代入して、

$$H_{\mathit{MZT}}(e^{j\omega})=g\frac{a_M}{b_N}\frac{\prod_{i=1}^M(e^{j\omega}-e^{\xi_iT})}{\prod_{i=1}^N(e^{j\omega}-e^{p_iT})}$$

これら周波数特性の差は、応答の比で表すと

$$H_{\mathrm{diff}}(\omega) = \frac{H_A(j\omega)}{H_{\mathit{MZT}}(e^{j\omega})}$$

この応答の比からサンプル点を$N$個取り出して離散的な応答の比とする。デジタル信号のサンプリング周期を$T$とすればナイキスト周波数は$\frac{\pi}{T}$であり、サンプル点の個数$N$を奇数として、取り出すサンプル点の周波数列$\omega_n$は

$$\omega_n = \frac{\frac{\pi}{T}}{N/2}\times n=\frac{2 \pi n}{NT}, \quad n=-\frac{N-1}{2}, \cdots,-2,-1,0,1,2,\cdots,\frac{N-1}{2}$$

これを代入した$H_{\mathrm{diff}}(\omega_n)に$逆離散フーリエ変換を行って補正FIRフィルタ$H_{\mathrm{diff}}(z)$を得る。

最終的に、

$$H_D(z)=H_{\mathit{MZT}}(z) \cdot H_{\mathrm{diff}}(z)$$

で、所望のデジタルフィルタを得る。

実際に作って見よう。MZTiでも扱ったピークフィルタを題材とする。

$$\begin{aligned} H_A(s) &= \frac{\omega_c^2-s^2}{s^2+2R\omega_cs+\omega_c^2} \\ &= -\frac{(s-\omega_c)(s+\omega_c)}{(s+R\omega_c-\omega_c\sqrt{R^2-1})(s+R\omega_c+\omega_c\sqrt{R^2-1})} \end{aligned}$$

MZTして、

$$\begin{aligned} H_{\mathit{MZT}}(z) &= -g\frac{(z-e^{\omega_cT})(z-e^{-\omega_cT})}{(z-e^{(-R\omega_c-\omega_c\sqrt{R^2-1})T})(z-e^{(-R\omega_c+\omega_c\sqrt{R^2-1})T})}\\ &= g\frac{-1+(e^{\omega_cT}+e^{-\omega_cT})z^{-1}-z^{-2}}{1 - (e^{(-R\omega_c+\omega_c\sqrt{R^2-1})T}+e^{(-R\omega_c-\omega_c\sqrt{R^2-1})T})z^{-1} + e^{-2 R\omega_c T }z^{-2}} \end{aligned}$$

そして$g$は

$$g = \frac{1 - (2 e^{-R \omega_c T} \cosh( \omega_c T \sqrt{R^2 - 1})) + e^{-2 R\omega_c T }}{-2+2 \cosh(ω_cT)}$$

$H_A(s)$と$H_{\mathit{MZT}}(z)$の周波数特性の比を考える。

$$H_{\mathrm{diff}}(\omega) = \frac{H_A(j\omega)}{H_{\mathit{MZT}}(e^{j\omega})}$$

サンプリング周波数を48000Hz、中心周波数を10000Hz、$R=0.2$にした時のの$H_A(s)$、$H_{\mathit{MZT}(z)}$および$H_{\mathrm{diff}}(\omega)$の振幅特性と位相特性は下図のようになる。橙色の特性を青色の特性に近づけるような特性、つまり緑色のグラフのようなものが$H_{\mathrm{diff}}(\omega)$になる。なお、MZTiでは振幅特性のみ補正したが、この方法では位相特性も補正できる。

（※MZTiで紹介したグラフと緑色のグラフの形が違うが、MZTiの項ではDigital - Analogを描画し、ここではAnalog - Digitalを描画しているからである。）

ここから先はグラフを用いながら説明する。

周波数サンプル点の個数を$N=63$として、$H_{\mathrm{diff}}(\omega)$をサンプリングする。先のグラフにサンプル点を描画すると次のようになる。横軸が対数なので等間隔には見えない。また、周波数サンプル点は負の周波数も含めて合計$N$個であり、下記のグラフは10Hz～24000Hzを描画しているので、サンプル点は31個しか書かれていない。

このサンプル点を逆離散フーリエ変換することで、補正FIRフィルタ$H_{\mathrm{diff}}(z)$を得る。補正FIRフィルタの特性は次のようなものである。

最後に$H_D(z)=H_{\mathit{MZT}}(z) \cdot H_{\mathrm{diff}}(z)$により、目的のデジタルフィルタを得る。青色がアナログフィルタの特性を示し、橙色が今回得たデジタルフィルタの特性を示す。アナログフィルタの特性を非常に良く近似出来ていることが分かる。

ちょっと乱暴な方法(双一次変換+オーバーサンプリング)

これまで色々試行錯誤を重ね、アナログフィルタの特性をデジタルフィルタに落とし込んできた。どの手法も、ナイキスト周波数付近で生じる特性の歪みを抑えようと工夫を凝らしてきた。しかし、オーバーサンプリングを行い、ナイキスト周波数を可聴域の外に移せば、特性の歪みも可聴域外に押し出される。そうすれば、アナログフィルタの特性により近いデジタルフィルタを、もっと簡単に実現できるのではないだろうか。

これは、正しい。双一次変換による特性の歪は、アナログフィルタの周波数を$\Omega$、デジタルフィルタの周波数を$\omega$として

$$\begin{aligned} \Omega &= \frac{2}{T} \tan \frac{\omega T}{2} \end{aligned}$$

と書けるのは前に扱った。このグラフを、異なるサンプリング周波数で描画してみよう。横軸と縦軸は10Hzから24000Hzの範囲を描画している。

このグラフから、双一次変換しただけの教科書的な実装のEQでも、サンプリング周波数が192kHzや384kHzであればほぼ理想的な関係$\Omega=\omega$が可聴域で得られることが分かる。2倍オーバーサンプリングして96kHzであったとしても、48kHzよりはマシである。

直前に扱ったMagpha EQの手法でも、サンプリング周波数を上げれば補正フィルタが掛かり始める帯域が上方向にずれていく。その結果、ほんのごく僅かながらもフィルタの形状が変わったり、数値演算上の誤差が変わって、音が理論上変わることはあり得る(知覚できるかどうかは置いておいて…)。

普通、「サンプリング周波数を変えたりオーバーサンプリングしたら音が変わった」系の話は、コンプレッサやサチュレータなど非線形歪みを発生させる処理や、エイリアシングフィルタのスペックが違うからという説明で語られることが多い。しかし、ここで話したように、非線形歪みを発生させないEQでも、音が変わることは全然ありえるという事実は頭の片隅に置いておくべきである。

Chapter 3. IIRデジタルフィルタの実装

Chapter 2.では、主にアナログフィルタの伝達関数$H(s)$からデジタルフィルタの伝達関数$H(z)$を得る方法について述べた。伝達関数$H(z)$は

$$H(z)=\frac{Y(z)}{X(z)} = \frac{\sum_{k=0}^{M} a_k z^{-k}}{1+\sum_{k=1}^{N} b_k z^{-k}}$$

のように、有理関数で表すことができる。このように表せたら、差分方程式

$$\begin{aligned} y(nT) &= a_0x(nT) + a_1x(nT-T) + \cdots + a_Mx(nT-MT) \\ & \qquad -(b_0y(nT) + b_1y(nT-T) + \cdots + b_My(nT-NT)) \\ &= \sum_{k=0}^{M} a_k x(nT-kT) - \sum_{k=1}^{N}b_ky(nT-kT) \\ \end{aligned}$$

を計算すれば、入力信号から出力信号を得ることができる。

つまり、これをプログラムで実装すればIIRデジタルフィルタの実装が完了するのだが、実はそう上手くは行かない。何故なら、コンピュータの世界では数値は有限桁数であり、計算には少なからず誤差が生じるからである。極端な場合その誤差が積み重なって、理論上安定なフィルタが不安定になることもあり得る。この章では、浮動小数点演算の誤差について説明した後、計算をブロックダイアグラムで図示するブロック法をまとめ、様々なデジタルフィルタの計算方法をブロックダイアグラムで示し、誤差を抑えたデジタルフィルタの計算方法を実験により調べる。

コンピュータの世界の数値表現（浮動小数点数）と演算誤差

コンピュータは2進数（0と1）で数値を表現する。例えば、10進数における$10$は、2進数では1010と表現される。コンピュータでは桁数を固定してメモリに値を格納することが多い。この桁数のことをbit幅と呼ぶ。例えば、44.1kHz/16bitで記録されたCD音質のWAVファイルは、1サンプルの振幅を16bitのbit幅、つまり16桁の0と1で表現している。10進数で$10$という振幅値を記録する場合、2進数では0000 0000 0000 1010となる（注1）。

次に小数の表現方法について紹介する。小数の表現には固定小数点と浮動小数点の2種類がある。

固定小数点は、小数点の位置を予め決めておく方法である。例えば、2進数で下位5bitは小数点以下の数字を表す…といった感じである。現代のパソコンでは浮動小数点で小数を扱うことが多い。 浮動小数点は、数を符号部・固定長の指数部・固定長の仮数部で表す方法である。10進数で例えるなら、$-123.456$を$-1\times1.23456\times10^{2}$と書くような感じである。

標準規格IEEE 754に則れば、32bit浮動小数点は符号部1bit・指数部8bit・仮数部23bitである。仮数部は、整数部分が1であるような2進小数の小数部分を示す。つまり、2進数で1010と表せる値を32bit浮動小数点で表現するなら、仮数部は先頭の1を除いて010となる。実質的に仮数部は24bitなので、32bit浮動小数点の10進数における有効桁数は$24\times\log_{10}2 \approx 7$桁である。

64bit浮動小数点は符号部1bit・指数部11bit・仮数部52bitである。10進数における有効桁数は$53\times\log_{10}2 \approx 16$桁である。

これから、浮動小数点演算で生じる代表的な演算誤差について、いくつか紹介する。

丸め誤差

浮動小数点数は有限桁の仮数部を持つため、無限に続く循環小数を正確に表現することはできない。例えば、10進数における$0.1$は2進数では$0.00011001100110011001100\cdots$という循環小数になるため、コンピュータでは有限のところで打ち切った値を持つことになる。

実際、0.1を32bit浮動小数点(float型)、64bit浮動小数点(double型)で定義して10進数で表示させると、

std::cout << std::setprecision(200) << 0.1f << std::endl;
// -> 0.100000001490116119384765625
std::cout << std::setprecision(200) << 0.1 << std::endl;
// -> 0.1000000000000000055511151231257827021181583404541015625

となり、どちらも0.1ではないことが分かる。これが及ぼす悪影響を見てみよう。

今、float型の0.1を1000回加算するプログラムを書いた。結果は100になるはずだが、実際はそうなっていない。

int main() {
  float x = 0.0;
  for (int i = 0; i < 1000; i++) {
    x += 0.1f;
  }
  std::cout << std::setprecision(200) << x << std::endl;
  // -> 99.99904632568359375

  return 0;
}

double型でも100にならない。

int main() {
  double x = 0.0;
  for (int i = 0; i < 1000; i++) {
    x += 0.1;
  }
  std::cout << std::setprecision(200) << x << std::endl;
  // -> 99.9999999999985931253831950016319751739501953125

  return 0;
}

情報落ち

$a+b$の演算を行う時、$a$と$b$の桁数の差が仮数部の桁数より大きい時に、小さな数の情報が失われてしまう。例えば、10進数で$12345600+1=12345601$という計算を有効桁数6桁で行うと、$1.23456\times10^{7}+1.00000=1.2345601\times10^{7}\approx1.23456\times10^{7}$となり、加算が消えてしまう。

下記のプログラムでは100000000に1を1000回加算しているが、結果は100000000のままである。

int main() {
  float x = 100000000.0;
  for (int i = 0; i < 1000; i++) {
    x += 1.0;
  }
  std::cout << std::setprecision(200) << x << std::endl;
  // -> 100000000

  return 0;
}

桁落ち

$a-b$の演算を行う時、$a$と$b$の値が極めて近い場合、有効桁数が減少してしまう。

下記のプログラムでは、$1.000001-1.0=0.000001=1\times10^{-6}$が本来あってほしい答えだが、そうなってない。

int main() {
  float a = 1.000001;
  float b = 1.0;
  float x = a - b;
  std::cout << std::setprecision(200) << x << std::endl;
  // -> 9.5367431640625e-07

  return 0;
}

ブロック法

この次にデジタルフィルタの計算方法をいくつか紹介したいが、式で表すと直感的でないことがある。そこで、計算の流れをブロックダイアグラムで示す「ブロック法」をまず紹介する。

デジタルフィルタの実現に必要な演算は加算・乗算・1サンプル遅延の3つであり、それぞれブロックで表すと次のようになる。

そして、これらブロックからなる1塊の回路を単位ブロックとみなして、それらの組み合わせを考える。単位ブロックの伝達関数を$F$と$G$とする。

縦続接続

伝達関数を考えれば、

$$\frac{X_1}{X_{\mathrm{in}}} = F, \quad \frac{X_{\mathrm{out}}}{X_1} = G$$

なので、

$$\frac{X_{\mathrm{out}}}{X_{\mathrm{in}}}=FG$$

並列接続

$X_{\mathrm{out}}=X_1+X_2=FX_{\mathrm{in}}+GX_{\mathrm{in}}$なので

$$\frac{X_{\mathrm{out}}}{X_{\mathrm{in}}}=F+G$$

帰還接続

関係式を立てると、

$$X_1=X_{\mathrm{in}}-X_2, \quad X_{\mathrm{out}}=X_1F, \quad X_2=GX_{\mathrm{out}}$$

ここから$X_1$と$X_2$を消去して

$$X_{\mathrm{out}}=F(X_{\mathrm{in}}-GX_{\mathrm{out}})$$

これを整理すれば

$$\frac{X_{\mathrm{out}}}{X_{\mathrm{in}}}=\frac{F}{1+FG}$$

直接形I (Direct Form I)

前提知識の章で、離散時間システムの入出力関係の式は差分方程式を用いて

$$\begin{aligned} y(nT) &= a_0x(nT) + a_1x(nT-T) + \cdots + a_Mx(nT-MT) \\ & \qquad -(b_0y(nT) + b_1y(nT-T) + \cdots + b_My(nT-NT)) \\ &= \sum_{k=0}^{M} a_k x(nT-kT) - \sum_{k=1}^{N}b_ky(nT-kT) \\ \end{aligned}$$

で表現でき、これの両辺をZ変換して変形を行い、

$$\begin{aligned} Y(z) &= \sum_{k=0}^{M} a_k z^{-k} X(Z) - \sum_{k=1}^{N} b_k z^{-k}Y(z) \\ \frac{Y(z)}{X(z)} &= \frac{\sum_{k=0}^{M} a_k z^{-k}}{1+\sum_{k=1}^{N} b_k z^{-k}} \end{aligned}$$

とした物が離散時間システムの伝達関数$H(z)$であった。つまり、アナログフィルタの伝達関数からデジタルフィルタの伝達関数$H(z)$を得たら、それの分母と分子の係数を使って差分方程式を実装すれば、デジタルフィルタの実装が完了する。この差分方程式をそのまま計算する方法を直接形Iと呼ぶ。 ブロックダイアグラムは下記のようになる。

直接形II (Direct Form II)

直接系Iの形を書き換える。加算の部分を境に2つに分割すれば、ブロック$F$と$G$の縦列接続だと分かるので入れ替えることが可能。最後にサンプル遅延ブロックを統合する。

縦続形構成(Biquad Direct Form I)

伝達関数は有理関数で表されるので、これを1次や2次の縦続接続と見なせるように変形を行う。

$$\begin{aligned} \frac{Y(z)}{X(z)} &= \frac{\sum_{k=0}^{M} a_k z^{-k}}{1+\sum_{k=1}^{N} b_k z^{-k}} \\ &= g \frac{1+\gamma_0z^{-1}}{1+\alpha_0z^{-1}}\prod_{i=1}^n\frac{1+\gamma_iz^{-1}+\delta_iz^{-2}}{1+a_iz^{-1}+\beta_iz^{-2}} \end{aligned}$$

この変形により、2次以上の高次フィルタを1次や2次フィルタのカスケード接続に変換している。 これを直接形Iで書く。

縦続形構成(Biquad Direct Form II)

直前に説明した構成を直接形IIで書く。

縦続形構成(Biquad Transposed Direct Form I)

回路中の信号の向きを逆にし、加算器を信号の分岐で置き換えた回路を転置回路という。実は転置回路は元の回路と伝達関数が一致する。

プログラムに実装するとBiquad Direct Form IIと同じ実装になる。

縦続形構成(Biquad Transposed Direct Form II)

同様に、Biquad Direct Form IIを転置する。

実験

実験を行い、どのフィルタ構成が最も良いか確認する。題材となるフィルタは4次のピークフィルタとする。中心周波数を$\omega_c$とする。

$$\begin{aligned} H(s) &= \left( \frac{\omega_c^2-s^2}{s^2+2R\omega_cs+\omega_c^2} \right)^2 \\ &= \frac{s^4-2\omega_c^2s^2+\omega_c^4}{s^4 + 4 R \omega_c s^3 + (2 \omega_c^2 + 4 R^2 \omega_c^2) s^2+ 4 R \omega_c^3 s + \omega_c^4} \end{aligned}$$

実験条件を単純にするために、双一次変換によりデジタルフィルタを得る。実験の際、カットオフ周波数はナイキスト周波数より十分低いものとする。

$$\begin{aligned} H(z) &= \left( \frac{\frac{\omega_c^2T^2-4}{4+4R\omega_cT+\omega_c^2T^2} + \frac{2\omega_c^2T^2+8}{4+4R\omega_cT+\omega_c^2T^2}z^{-1}+\frac{\omega_c^2T^2-4}{4+4R\omega_cT+\omega_c^2T^2}z^{-2}}{1+\frac{-8+2\omega_c^2T^2}{4+4R\omega_cT+\omega_c^2T^2}z^{-1}+\frac{4-4R\omega_cT+\omega_c^2T^2}{4+4R\omega_cT+\omega_c^2T^2}z^{-2}} \right)^2 \\ &= \frac{a_0+a_1z^{-1}+a_2z^{-2}}{1+b_1z^{-1}+b_2z^{-2}}\cdot \frac{a_0+a_1z^{-1}+a_2z^{-2}}{1+b_1z^{-1}+b_2z^{-2}} \end{aligned}$$

直接形IとIIは次式のパラメータ$A_0, A_1, A_2, A_3, A_4, B_1, B_2, B_3, B_4$を使う。

$$\begin{aligned} H(z) &= \frac{a_0+a_1z^{-1}+a_2z^{-2}}{1+b_1z^{-1}+b_2z^{-2}}\cdot \frac{a_0+a_1z^{-1}+a_2z^{-2}}{1+b_1z^{-1}+b_2z^{-2}} \\ &= \frac{{a_0^2 + 2 a_0 a_1z^{-1} + (2 a_0 a_2 + a_1^2)z^{-2} + 2 a_1 a_2z^{-3} + a_2^2z^{-4}}}{{1 + 2 b_1z^{-1} + (2 b_2 + b_1^2)z^{-2} + 2 b_1 b_2z^{-3} + b_2^2z^{-4} }} \\ &= \frac{A_0+A_1z^{-1}+A_2z^{-2}+A_3z^{-3}+A_4z^{-4}}{1+B_1z^{-1}+B_2z^{-2}+B_3z^{-3}+B_4z^{-4}} \end{aligned}$$

Biquad DFI・IIとBiquad Transposed DFI・IIは次式のパラメータ$g, \alpha_0, \beta_0, \gamma_0, \delta_0, \alpha_1, \beta_1, \gamma_1, \delta_1$を使う。

$$\begin{aligned} H(z) &= \frac{a_0+a_1z^{-1}+a_2z^{-2}}{1+b_1z^{-1}+b_2z^{-2}}\cdot \frac{a_0+a_1z^{-1}+a_2z^{-2}}{1+b_1z^{-1}+b_2z^{-2}} \\ &= a_0^2 \frac{1+\frac{a_1}{a_0}z^{-1}+\frac{a_2}{a_0}z^{-2}}{1+b_1z^{-1}+b_2z^{-2}}\cdot \frac{1+\frac{a_1}{a_0}z^{-1}+\frac{a_2}{a_0}z^{-2}}{1+b_1z^{-1}+b_2z^{-2}} \\ &=g\frac{1+\gamma_0z^{-1}+\delta_0z^{-2}}{1+\alpha_0z^{-1}+\beta_0z^{-2}} \cdot \frac{1+\gamma_1z^{-1}+\delta_1z^{-2}}{1+\alpha_1z^{-1}+\beta_1z^{-2}} \end{aligned}$$

サンプリング周波数を48000Hzとし、中心周波数を20Hzと200Hz、$R$を0.01と0.1に変化させ、各フィルタ構成のアナログフィルタに対する精度を確認する。演算は全て64bit浮動小数点で行っている。

結果を下記に示す。上から3段目まで、左側のグラフが振幅特性、右側のグラフが位相特性を表す。各1段目のグラフが特性をそのまま示し、2段目のグラフは特性の差(Analog - Digital)、3段目のグラフは特性の差を対数表示したものである。4段目左右のグラフは1000Hz振幅0.0dBFSのサイン波を各実装に処理させ、結果の振幅特性を示した物である。2段目のグラフでは、プロットが0に近い程アナログフィルタと特性が近く、3段目のグラフでは、プロットが$-\infty$に近い程アナログフィルタと特性が近い。4段目のグラフはプロットが$-\infty$に近いほど低歪で高精度である。 なお、今まで示してきたグラフとは違い、横軸が1Hz～24000Hzの範囲を描画している点に注意せよ。

全体の傾向を見ると、

• 実装方法に関わらず、
  • 中心周波数が高い方がアナログフィルタの特性を高精度に再現出来る
  • $R$が小さい(=ゲインが低く、Qが緩い)方がアナログフィルタの特性を高精度に再現出来る
• 実装方法の優劣を判断すると
  • 直接形(DFI・DFII)に比べて縦続形(Biquad○○)の方がアナログフィルタの特性再現度が高い
  • 直接形(DFI・DFII)に比べて縦続形(Biquad○○)の方がTHD性能が良い
  • 縦続形(Biquad○○)の中でも、Biquad DFIとBiquad Transposed DF2のTHD性能が良い

と分かる。 グラフ中にBiquad Direct Form IIが見えないが、Biquad Transposed Direct Form Iと完全に重なっている(同じ実装なため)。

---

---

---

Chapter 4. アナログフィルタの伝達関数

Chapter 2.ではアナログフィルタの伝達関数$H(s)$からデジタルフィルタの伝達関数$H(z)$を得る方法について述べ、Chapter 3.ではデジタルフィルタの伝達関数$H(z)$をプログラムに実装する方法を述べた。最後に、アナログフィルタの伝達関数$H(s)$を入手する方法を紹介する。

とは言っても、アナログフィルタの伝達関数はインターネット上に多く落ちているため、必要があれば検索した方が早い。ここでは特に代表的な物を紹介するに留める。

アナログ1次フィルタ

Chapter 3.の実験結果から、高次のデジタルフィルタを作る際は1次か2次フィルタのカスケード接続をしたほうが性能が良いことを示した。似たようなことはアナログ回路でも言うことができ、高次のフィルタを作る際は低次のアクティブフィルタを多段縦続接続することが多い。

低域通過フィルタ

カットオフ周波数を$\omega_c$とすると

$$H_{\mathrm{LPF}}(s) = \frac{\omega_c}{s+\omega_c}$$

振幅特性は

$$|H_{\mathrm{LPF}}(j\omega)| = \frac{\omega_c}{\sqrt{\omega^2+\omega_c^2}}$$

であり、カットオフ周波数の地点で振幅が$\frac{1}{\sqrt{2}}$、dBで表すなら$20\log_{10}\frac{1}{\sqrt{2}}\approx-3.010\ [\mathrm{dB}]$になる。

スロープは6dB/octである。

高域通過フィルタ

周波数対数スケール上で考えると、高域通過フィルタは低域通過フィルタの特性を$\omega_c$で左右反転したものである。そこで、$\frac{s}{\omega_c} \rightarrow \frac{\omega_c}{s}$という変数変換を行うことで、

$$\begin{aligned} H_{\mathrm{HPF}}(s) = \frac{s}{s+\omega_c} \end{aligned}$$

振幅特性は

$$|H_{\mathrm{HPF}}(j\omega)| = \frac{\omega}{\sqrt{\omega^2+\omega_c^2}}$$

であり、カットオフ周波数の地点で振幅が$\frac{1}{\sqrt{2}}$、dBで表すなら$20\log_{10}\frac{1}{\sqrt{2}}\approx-3.010\ [\mathrm{dB}]$になる。

スロープは6dB/octである。

低域シェルビングフィルタ

DC成分のゲインを$A_0$とすれば伝達関数は

$$H_{\mathrm{LSF}}(s) = \frac{s+A_0\omega_c}{s+\omega_c}$$

周波数特性は

$$|H_{\mathrm{LSF}}(j\omega)| = \frac{\sqrt{\omega^2+A_0^2\omega_c^2}}{\sqrt{\omega^2+\omega_c^2}}$$

形が非対称であり、傾斜も選択できないためデジタルEQではあまり見かけない。

高域シェルビングフィルタ

低域通過フィルタから高域通過フィルタを求めたように、$\frac{s}{\omega_c} \rightarrow \frac{\omega_c}{s}$という変数変換を行えば、

$$H_{\mathrm{HSF}}(s) = \frac{A_{\infty}s+\omega_c}{s+\omega_c}$$

を得る。$A_{\infty}$は高域でのゲインを表す。

周波数特性は

$$|H_{\mathrm{HSF}}(j\omega)| = \frac{\sqrt{A_{\infty}^2\omega^2+\omega_c^2}}{\sqrt{\omega^2+\omega_c^2}}$$

形が非対称であり、傾斜も選択できないためデジタルEQではあまり見かけない。

アナログ2次フィルタ

アナログ2次フィルタは、1次フィルタに比べると形とパラメータにバリエーションが出てくる。

低域通過フィルタ

カットオフ周波数を$\omega_c$、Qを$Q$とすると、伝達関数は

$$H_{\mathrm{LPF}}(s) = \frac{\omega_c^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{LPF}}(j\omega)| = \frac{\omega_c^2}{\sqrt{(\omega^2-\omega_c^2)^2+\frac{\omega_c^2\omega^2}{Q^2}}}$$

$Q=\frac{1}{\sqrt{2}}\approx0.707$の時、フィルタはButterworth特性となる。この時、カットオフ周波数地点の振幅特性は$\frac{1}{\sqrt{2}}$、dBで表すなら$20\log_{10}\frac{1}{\sqrt{2}}\approx-3.010\ [\mathrm{dB}]$になる。

また、$Q=\frac{1}{\sqrt{3}}\approx0.577$の時、フィルタはBessel特性となる。

高域通過フィルタ

物凄く寄り道してしまった。

周波数対数スケール上で考えると、高域通過フィルタは低域通過フィルタの特性を$\omega_c$で左右反転したものである。そこで、$\frac{s}{\omega_c} \rightarrow \frac{\omega_c}{s}$という変数変換を行うことで、

$$H_{\mathrm{HPF}}(s) = \frac{s^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{HPF}}(j\omega)| = \frac{\omega^2}{\sqrt{(\omega^2-\omega_c^2)^2+\frac{\omega_c^2\omega^2}{Q^2}}}$$

$Q=\frac{1}{\sqrt{2}}\approx0.707$の時、フィルタはButterworth特性となる。この時、カットオフ周波数地点の振幅特性は$\frac{1}{\sqrt{2}}$、dBで表すなら$20\log_{10}\frac{1}{\sqrt{2}}\approx-3.010\ [\mathrm{dB}]$になる。

また、$Q=\frac{1}{\sqrt{3}}\approx0.577$の時、フィルタはBessel特性となる。

帯域通過フィルタ

伝達関数は

$$H_{\mathrm{BPF}}(s) = \frac{\frac{\omega_c}{Q} s}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{BPF}}(j\omega)| = \frac{\frac{\omega_c}{Q} \omega}{\sqrt{(\omega^2-\omega_c^2)^2+\frac{\omega_c^2\omega^2}{Q^2}}}$$

伝達関数を$H_{\mathrm{BPF}}(s) = \frac{\omega_c s}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$にすることで、BPFのスカートのゲイン位置を固定したまま、頂点のゲインを変えられるフィルタにもできる。

このようなフィルタになっても、$Q=\frac{1}{\sqrt{3}}\approx0.577$の時、フィルタはBessel特性となる。

Notchフィルタ

帯域除去フィルタとも言う。

伝達関数は

$$H_{\mathrm{BRF}}(s) = \frac{s^2 + \omega_c^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{BRF}}(j\omega)| = \frac{|-\omega^2+\omega_c^2|}{\sqrt{(\omega^2-\omega_c^2)^2+\frac{\omega_c^2\omega^2}{Q^2}}}$$

位相特性のグラフを見ても分かる通り連続では無くなるため、$Q=\frac{1}{\sqrt{3}}\approx0.577$でBessel特性は持たない

全域通過フィルタ

伝達関数は

$$H_{\mathrm{APF}}(s) = \frac{s^2-\frac{\omega_c}{Q}s+\omega_c^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{APF}}(j\omega)| = 1$$

$Q=\frac{1}{\sqrt{3}}\approx0.577$でBessel特性を持つ。

ピーキングフィルタ

ピークのゲインを$A$として、伝達関数は

伝達関数は

$$H_{\mathrm{PKF}}(s) = \frac{s^2+\frac{\sqrt{A}}{Q}\omega_cs+\omega_c^2}{s^2+\frac{\omega_c}{Q\sqrt{A}}s+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{PKF}}(j\omega)| = \frac{\sqrt{(\omega^2-\omega_c^2)^2+\frac{A\omega_c^2\omega^2}{Q^2}}}{\sqrt{(\omega^2-\omega_c^2)^2+\frac{\omega_c^2\omega^2}{AQ^2}}}$$

$Q=\frac{1}{\sqrt{3}}\approx0.577$でBessel特性を持つとは言えなさそう。

低域シェルビングフィルタ

伝達関数は

$$H_{\mathrm{LSF}}(s) = \sqrt{A}\frac{s^2+\frac{\sqrt[4]{A}}{Q}\omega_cs+\omega_c^2\sqrt{A}}{\sqrt{A}s^2+\frac{\sqrt[4]{A}}{Q}\omega_cs+\omega_c^2}$$

振幅特性は

$$|H_{\mathrm{LSF}}(j\omega)| = \sqrt{A}\frac{\sqrt{(\omega^2-\omega_c^2\sqrt{A})^2+\frac{\sqrt{A}}{Q^2}\omega_c^2\omega^2}}{\sqrt{(\sqrt{A}\omega^2-\omega_c^2)^2+\frac{\sqrt{A}}{Q^2}\omega_c^2\omega^2}}$$

$Q=\frac{1}{\sqrt{2}}\approx0.707$でButterworth特性を持つ。 $Q=\frac{1}{\sqrt{3}}\approx0.577$でBessel特性を持つ。

高域シェルビングフィルタ

伝達関数は

$$H_{\mathrm{HSF}}(s) = \sqrt{A}\frac{\sqrt{A}s^2+\frac{\sqrt[4]{A}}{Q}\omega_cs+\omega_c^2}{s^2+\frac{\sqrt[4]{A}}{Q}\omega_cs+\omega_c^2\sqrt{A}}$$

振幅特性は

$$|H_{\mathrm{HSF}}(j\omega)| = \sqrt{A}\frac{\sqrt{(\sqrt{A}\omega^2-\omega_c^2)^2+\frac{\sqrt{A}}{Q^2}\omega_c^2\omega^2}}{\sqrt{(\omega^2-\omega_c^2\sqrt{A})^2+\frac{\sqrt{A}}{Q^2}\omega_c^2\omega^2}}$$

$Q=\frac{1}{\sqrt{2}}\approx0.707$でButterworth特性を持つ。 $Q=\frac{1}{\sqrt{3}}\approx0.577$でBessel特性を持つ。

アナログ高次フィルタ

HPFとLPFは高次フィルタが望まれる場面がある。そこで高次のデジタルフィルタを作るために1次か2次フィルタを単純に縦続接続すると、Qが同じなのにカットオフ周波数地点のゲインが変わってしまう。例えば3次のLPFを作る為に、下記の2つを縦続接続すると次のような結果が得られる。

$$H_{\mathrm{1,\ LPF}}(s) = \frac{\omega_c}{s+\omega_c}$$ $$H_{\mathrm{2,\ LPF}}(s) = \frac{\omega_c^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$ $$H_{\mathrm{3,\ NaiveLPF}}(s) = \frac{\omega_c}{s+\omega_c} \cdot \frac{\omega_c^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

これはあまり望ましくない。そこで、Butterworth特性を基準にQを考え直してみる。「Butterworth特性とは」の項で説明したように、Butterworth多項式を使えば、LPFがButterworth特性を持つ時の分母の係数を知ることができる。

$$B_n(s) = \begin{cases} \prod_{k=1}^{\frac{n}{2}} (s^2 - 2s \cos \left( \frac{(2k+n-1)\pi}{2n} \right) + 1) & \quad (nが偶数) \\ (s+1)\prod_{k=1}^{\frac{n-1}{2}} (s^2 - 2s \cos \left( \frac{(2k+n-1)\pi}{2n} \right) + 1) & \quad (nが奇数) \\ \end{cases}$$

3次Butterworth多項式は$B_3(s) = (s+1)(s^2+s+1)$であり、これはカットオフ周波数で正規化されていることを踏まえれば、3次Butterworth低域通過フィルタの伝達関数は

$$H_{\mathrm{3,\ ButterworthLPF}}(s) = \frac{\omega_c}{s+\omega_c} \cdot \frac{\omega_c^2}{s^2+\omega_cs+\omega_c^2}$$

である。$H_{\mathrm{2,\ LPF}}(s)$の$Q$パラメータの位置に倣って、これに$Q$パラメータを追加すると

$$H_{\mathrm{3,\ LPF}}(s) = \frac{\omega_c}{s+\omega_c} \cdot \frac{\omega_c^2}{s^2+\frac{\omega_c}{Q}s+\omega_c^2}$$

を得る。$H_{\mathrm{3,\ LPF}}(s)$は$Q=1$でButterworth特性を持つ。$H_{\mathrm{2,\ LPF}}(s)$も$Q=1$でButterworth特性を持つようにしたいので、2次Butterworth多項式は$B_2(s) = s^2+\sqrt{2}s+1$であることから、

$$H_{\mathrm{2,\ LPF}}(s) = \frac{\omega_c^2}{s^2+\sqrt{2}\frac{\omega_c}{Q}s+\omega_c^2}$$

と定義しなおせば、$H_{\mathrm{2,\ LPF}}(s)$も$Q=1$でButterworth特性を持つようになる。図示すると次のようになり、$Q$が同じならカットオフ周波数地点のゲインも同じになる。

4次以降の低域通過フィルタを作る時は$Q$パラメータを追加する位置を工夫する。例えば、4次Butterworth多項式は下記であるから、

$$B_4 = \left(s^2-2s \cos \left(\frac{5}{8}\pi\right)+1 \right)\left(s^2-2s \cos \left(\frac{7}{8}\pi\right)+1\right)$$

4次Butterworth低域通過フィルタの伝達関数は、

$$H_{\mathrm{4,\ ButterworthLPF}}(s) = \frac{\omega_c^2}{s^2-2\omega_cs \cos \left(\frac{5}{8}\pi\right)+\omega_c^2} \cdot \frac{\omega_c^2}{s^2-2\omega_cs \cos \left(\frac{7}{8}\pi\right)+\omega_c^2}$$

となる。$Q$パラメータは、Butterworth多項式で$k=1$の項に付与する。

$$H_{\mathrm{4,\ LPF}}(s) = \frac{\omega_c^2}{s^2- \frac{2\omega_cs \cos \left(\frac{5}{8}\pi\right)}{Q}+\omega_c^2} \cdot \frac{\omega_c^2}{s^2-2\omega_cs \cos \left(\frac{7}{8}\pi\right)+\omega_c^2}$$

同様に$H_{\mathrm{5,\ LPF}}(s)$と$H_{\mathrm{6,\ LPF}}(s)$は

$$\begin{aligned} H_{\mathrm{5,\ LPF}}(s) &= \frac{\omega_c}{s+\omega_c} \cdot \frac{\omega_c^2}{s^2- \frac{2\omega_cs \cos \left(\frac{3}{5}\pi\right)}{Q}+\omega_c^2} \cdot \frac{\omega_c^2}{s^2-2\omega_cs \cos \left(\frac{4}{5}\pi\right)+\omega_c^2} \\ H_{\mathrm{6,\ LPF}}(s) &= \frac{\omega_c^2}{s^2- \frac{2\omega_cs \cos \left(\frac{7}{12}\pi\right)}{Q}+\omega_c^2} \cdot \frac{\omega_c^2}{s^2-2\omega_cs \cos \left(\frac{3}{4}\pi\right)+\omega_c^2} \cdot \frac{\omega_c^2}{s^2-2\omega_cs \cos \left(\frac{11}{12}\pi\right)+\omega_c^2} \end{aligned}$$

図示すると次のようになる。$Q$に依らず、カットオフ周波数の地点のゲインが同じである。1次LPFは$Q$パラメータが存在しないので、これは例外である。

Chapter 5. まとめ

IIRなデジタルEQは

1. Chapter 4. を元にアナログフィルタの伝達関数$H(s)$を選択
2. Chapter 2. を元に$H(s)$からデジタルフィルタの伝達関数$H(z)$を獲得
3. Chapter 3. を元にプログラム実装

の手順で作ることができる。Chapter 2.～4. の背景知識としてChapter 1.の内容を使っている。

おすすめ文献

本記事を書くにあたって多くの本・サイトを参考にした。閲覧した中で特に役立ったものを下記に紹介する。

• ディジタル信号処理システムの基礎 (森北出版 2008)
  • 著者: 渡部英二
  • 大学で使っていた教科書。
• 応用解析学の基礎 (森北出版 2014)
  • 著者: 坂和正敏
  • 大学で使っていた教科書。
• やる夫で学ぶディジタル信号処理
  • http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/main.html
  • デジタル信号処理を扱った日本語Webサイトの中で最も有名なもの。
• Vadim Zavalishin. The Art of VA Filter Design 2.1.2. (2020)
  • https://www.kvraudio.com/forum/viewtopic.php?t=350246
  • 本記事では扱わなかったTPT法による$H(z)$の入手法も紹介している。
• Robert Bristow-Johnson, Cookbook formulae for audio equalizer biquad filter coefficients
  • https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
  • 通称「RBJ Audio-EQ-Cookbook」。伝達関数$H(s)$のレシピ集
• Introduction to Digital Filters: with Audio Applications (2007)
  • https://ccrma.stanford.edu/~jos/filters/
  • デジタル信号処理の教科書。無料で読めるのに文量が多く、充実している。
  • 同じWebサイトにて、音声信号処理に関する他の教科書も閲覧できる。