ディジタル信号処理第09回ディジタル IIR フィルタ

💡

IIRのアナログフィルタを設計できた。ディジタル計算機でこのフィルタを使いたいが、どうしたらいいだろうか？ 2 つの方法を紹介する。

0. アナログフィルタ vs. ディジタルフィルタ

復習１：フィルタの種類

低域通過フィルタ (low pass filter, LPF)：低い周波数を通過 (pass) させ，高い周波数を遮断 (cut) するフィルタ．高い周波数を遮断することを強調してハイカットフィルタとも言う

高域通過フィルタ (high pass filter, HPF)：高い周波数を通過させ，低い周波数を遮断するフィルタ。低い周波数を遮断することを強調してローカットフィルタとも言う

帯域通過フィルタ (band pass filter, BPF)：特定の中間帯域を通過させるフィルタ

帯域遮断フィルタ (band elimination filter, BEF)：特定の中間帯域を遮断させるフィルタ

全域通過フィルタ (all pass filter, APF)：全周波数帯域において同一ゲインで通貨させるフィルタ

アナログとディジタルの違い：フィルタにおいて

https://www.youtube.com/watch?v=5tc14WHUoMw

アナログフィルタ (動画では6秒あたりのペダル)。比較的安価な部品で作れ、連続信号のまま扱える（遅延が非常に小さく）。

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

ディジタルフィルタ。変更をかけやすい。数値誤差はのるが、設計の自由度が高いため誤差を小さくしやすい

💡

アナログフィルタで良いものができたとして、それと同じ効果を持つディジタルフィルタを作れないか？今回は IIR フィルタに限定して議論。

2つの方法を紹介。ざっくりいうと

インパルス不変法：めちゃくちゃ簡単だけど使用に注意が必要

双一次変換：安全に使えるけどちょっと難しい

復習２：インパルス応答

　あるシステムを考える。その特性はインパルス応答 (impulse response) で記述できる。インパルス応答とは、デルタ関数で記述されるインパルス信号に対して、そのシステムが出力する時間波形である。

　デルタ関数を離散時間で表現すると

\begin{align} \delta(n) = \left\{\begin{array}{ll}1 & n = 0 \\0 & n \neq 0\end{array}\right. \end{align}

同様に連続時間で表すと

\begin{align} \delta(t) = \left\{\begin{array}{ll}\infty & t = 0 \\0 & t \neq 0\end{array}\right.,\;\; \int_{t=-\infty}^{\infty} \delta(t) dt = 1 \end{align}

である。ちなみに、上式の定義は厳密には不正確である。厳密な定義を知りたい方は「デルタ関数　超関数」などで調べるとよい。

💡

インパルス信号は、すべての周波数成分が 1 である（すなわち、全周波数の成分を同程度に持っている）信号である。

　入力信号が与えられたときの出力信号は畳み込みで計算できる。離散時間の入力信号を $x(n)$ 、出力信号を $y(n)$ とすると、その関係は

\begin{align} y(n) = h(n) * x(n) \end{align}

で計算できる。 $x(n) = \delta(n)$ 、すなわちインパルス信号が入力されたとき、

y(n) = h(n) * \delta(n) = h(n)

となる。このときの出力 $h(n)$ がインパルス応答と呼ばれる。連続時間においても同様。

復習３：FIR フィルタと IIR フィルタ

伝達関数の形によってフィルタ（システム）は 2 種類に分類される。式 (3) を z 変換すると

\begin{align} Y(z) = H(z) X(z),\;\; H(z) = \frac{Y(z)}{X(z)} \end{align}

となる。連続時間信号の場合はラプラス変換 (s変換) であり、上式の z が s になった式が導かれる。この $H(z)$ がシステムの伝達関数である。この伝達関数の式の形によって分類すると以下のようになる。

IIR (infinite impulse response) フィルタ

過去の時刻の出力が入力側に戻って、現在の時刻の出力を決める形

H(z) = \frac{b_0 + \sum\limits_{k = 1}^{J}b_k z^{-k}}{1+\sum\limits_{k = 1}^{I}a_kz^{-k}}

FIR (finite impulse response) フィルタ

出力は入力側に戻らず、入力のみで現在の時刻の出力が決まる形

H(z) = b_0 + \sum\limits_{k = 1}^{J}b_k z^{-k}

1. インパルス不変法 (impulse invariant method)

アナログフィルタからディジタルフィルタを求める

　インパルス不変法とは、アナログフィルタとディジタルフィルタで、インパルス信号を入力したときの出力が変わらないよう（不変になるよう）に設計する方法である。不変であることを図を用いて説明する。

上図において、 $x(n)$ はサンプリング周期 $T$ で $x(t)$ を離散化したものとする。すなわち、 $x(n) = x(nT)$ である。二つのシステムがありそれぞれの出力を $y(t), y(n)$ としたとき、 $y(n) = y(nT)$ であることをここでは不変と呼ぶ。 $h_D(n) = h_A(nT)$ が不変を満たすことのは自明である。

　上記を前提に式を展開する。連続時間システムのインパルス応答 $h_A(t)$ の伝達関数を $H_A(s)$ とし、この伝達関数が

\begin{align} H_A(s) = \frac{a_1}{s-s_1} + \frac{a_2}{s-s_2} + \ldots + \frac{a_N}{s-s_N} \end{align}

で書かれるとする。逆ラプラス変換を用いるとインパルス応答は

h_A(t) = a_1 e^{s_1t} + a_2 e^{s_2t} + \ldots + a_N e^{s_Nt}

で記述される。離散時間システムのインパルス応答 $h_D(n)$ は、このインパルス応答をサンプリング周期 $T$ で離散化すればよいため、

h_D(n) = h_A(nT)= a_1 e^{s_1nT} + a_2 e^{s_2nT} + \ldots + a_N e^{s_NnT}

である。これを z 変換すると、離散時間システムの伝達関数 $H_D(z)$ が求められる。伝達関数を求めるために、インパルス応答の第 $i$ 項に着目すると、

\begin{align} \mathcal{Z}[e^{s_i nT}] = \sum\limits_{n = 0}^{\infty}e^{s_i nT} z^{-n} = \frac{1}{1 - e^{s_iT}z^{-1}} \end{align}

である。これをすべての項に適用すると、

H_D(z) = \frac{a_1}{1 - e^{s_1T}z^{-1}} + \frac{a_2}{1 - e^{s_2T}z^{-1}} + \ldots + \frac{a_N}{1 - e^{s_NT}z^{-1}}

が得られる。

💡

H_A(s) = \frac{a_1}{s-s_1}

のアナログフィルタがある。インパルス不変法を用いてディジタルフィルタに変換し、伝達関数を求めよ。また、その周波数特性の実部と虚部を求めよ。

💡

H_A(s) = \frac{b}{(s+a)^2+b^2}

のアナログフィルタがある。インパルス不変法を用いてディジタルフィルタに変換し、伝達関数を求めよ。また、その周波数特性の実部と虚部を求めよ。ヒント：式(5)の形に変形できればよさそう。

求めたディジタルフィルタは安定なのか？

　フィルタの安定性は極 (pole) で判定できる。いま、アナログフィルタの $i$ 番目の極を

s_i = \alpha_i + j\beta_i

とする。ラプラス変換で全ての $\alpha_i$ について $\alpha_i < 0$ であれば、そのフィルタは安定である。式(6) より、対応するディジタルフィルタの伝達関数の極 $z_i$ は、

z_i = e^{s_iT} = e^{(\alpha_i + j\beta_i)T}

である。z 変換で全ての $z_i$ について $\left| z_i \right| < 1$ であれば、そのフィルタは安定である。 $\left| z_i \right|$ は

\left|z_i\right| = \left|e^{(\alpha_i + j\beta_i)T}\right| = \left|e^{\alpha_iT}\right|\cdot \left|e^{j\beta_iT}\right| = \left|e^{\alpha_iT}\right|

である。 $\alpha_i < 0$ のとき $\left| z_i \right| < 1$ である。すなわち、元のアナログフィルタが安定であれば、インパルス不変法により導出されるディジタルフィルタも同様に安定である。

💡

同様の展開により、s 平面における虚数軸 (

\alpha = 0

) は、z 平面における単位円 (

|z| = 1

) に写像されることになる。

どんなアナログフィルタでも変換できるのか？

　インパルス不変法は、アナログフィルタのインパルス応答の離散化に基づいている。そのため、サンプリング定理を満たさなければならない。元のアナログフィルタの最大周波数が $f_{\mathrm max}$ [Hz] であるとき、その 2 倍以上のサンプリング周波数を用いなければならない。逆に、サンプリング周波数があらかじめ決まっている場合は、使用できるアナログフィルタの最大周波数は限定されることになる。

　この理由より、インパルス不変法で変換できるフィルタはLPF、BPFに限られる。また、その最大周波数（ここでは、阻止域に至る周波数ととらえてよい）がサンプリング定理を満たしていなければならない。

2. 双一次変換 (Bilinear transform)

復習：ラプラス変換 (s 変換) と z 変換の関係

　双一次変換の前にラプラス変換と z 変換の関係を思い出そう。ラプラス変換と z 変換はそれぞれ

\begin{align} X(s) &= \int_{0}^{\infty}x(t)s^{-st}dt \\ X(z) &= \sum\limits_{n = 0}^{\infty}x(nT) z^{-n} \end{align}

この 2 つの式において、離散時間信号は連続時間信号をサンプリング周期 $T$ でサンプリングしたものである。この 2 式が等価となるためには、 $z = e^{sT}$ でなければならない。（この論理はやや突飛しているため、正確な論理展開は z 変換の回を参照されたい）

　この関係式を、さらにテイラー展開による1次近似で表すと、

\begin{align} z &= e^{sT} \\ &= \frac{e^{sT/2}}{e^{-sT/2}} \\ &\simeq \frac{1 + sT/2}{1 - sT/2} \end{align}

となる。この式を s について解くと

\begin{align} s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}} \end{align}

が得られる。この式(11)(12) を、s と z の双一次変換と呼ぶ。双一次変換を用いてアナログフィルタをディジタルフィルタに変換する場合は、アナログフィルタの伝達関数の周波数 s をこれらの式を用いて z に変換すればよい（詳しい手順は後述）。

双一次変換の意味するところ：周波数の対応

　式 (11) に示すように、連続時間信号のすべての周波数は離散時間信号の周波数に対応することになる。では、どのような対応になるだろうか。いま、連続時間信号の角周波数を $\omega_A$ 、離散時間信号の角周波数を $\omega_D$ とし、双一次変換の式に $s = j\omega_A, z = e^{j\omega_DT}$ を代入すると、

\begin{align} j\omega_A &= \frac{2}{T}\frac{1-e^{-j\omega_DT}}{1 + e^{j\omega_DT}} \\ \omega_A &= \frac{2}{T} \tan\left( \omega_D T /2 \right) \end{align}

となる。この式が意味するのは 2 つである。

連続時間信号の角周波数は、離散時間信号の角周波数と非線形に対応する。なお、インパルス不変法では線形に対応に対応する。

連続時間信号のすべての周波数は、arctan によって一定範囲に変換され、離散時間信号の周波数となる。時間の都合上省略するが、この効果によってインパルス不変法で触れたないキスと周波数の問題は発生しない。

求めたディジタルフィルタは安定なのか？

結論から言えば、インパルス不変法と同様に、元のアナログフィルタが安定であれば、双一次変換で求めたディジタルフィルタも安定である。

💡

元のアナログフィルタが安定であれば、双一次変換で求めたディジタルフィルタも安定であることを証明せよ。具体的には、式(11) に

s = \sigma + j\omega

を代入し、

\sigma < 0

であれば

|z| < 1

であることを導け。

アナログフィルタからディジタルフィルタを求める

具体的な数字を用いて設計してみよう。カットオフ角周波数（遮断角周波数）が $\omega_A$ のアナログフィルタ

H(s) = \frac{\omega_A}{s + \omega_A}

がある。これをもとに、サンプリング周波数 8 kHz、カットオフ周波数（遮断周波数）が 1 kHz のディジタルフィルタを作りたい。この時のアナログフィルタの特性（すなわちカットオフ周波数 $\omega_A$ ）と、ディジタルフィルタの伝達関数を求めよう。

step 1: アナログフィルタの特性を求める

式(14) を用いる。 $T= 1/8000, \omega_D = 2\pi\cdot1000$ (カットオフ角周波数) を代入すると、

\omega_A = 2\cdot8000\cdot\tan\left(2\pi\frac{1000}{8000\cdot2} \right) \simeq 6627

故に、アナログフィルタは

\begin{align} H(s) = \frac{6627}{s + 6627} \end{align}

となる。なお、このカットオフ角周波数からカットオフ周波数を求めると

6627 \mathrm{[rad/s]} / 2\pi \mathrm{[rad]} \simeq 1055\mathrm{[Hz]}

となり、アナログフィルタのカットオフ周波数とディジタルフィルタのカットオフ周波数が異なることを確認できる。

step 2: 双一次変換を用いてディジタルフィルタの伝達関数を求める

式(15)に式(12)を代入し整理すればよい。

H_D(z) = \left.H_A(s)\right|_{s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}} = \frac{6627(1 + z^{-1})}{22627 - 9373z^{-1}}

3. まとめ

アナログIIRフィルタをディジタルフィルタに変換したい

インパルス不変法は、アナログフィルタをサンプリングする方法である。単純な方法で変換できるが、サンプリング定理を満たす範囲のフィルタのみを扱える。

双一次変換は、アナログフィルタの周波数をディジタルフィルタの周波数に対応させる方法である。サンプリング定理を自動的に満たすが、周波数の対応が非線形であることに注意する。