15. ディジタルフィルタの解析

15.1 ディジタルフィルタの周波数特性

やらない夫

さて，z 変換を使って実際にディジタルフィルタの特性を解析していこうと思う．

やる夫

んー，特性っていっても抽象的すぎてピンと来ないお．

やらない夫

一口に特性といってもいろいろあるが，前々々回話した通り，ここで我々が考えたいのは周波数応答だ．つまり，ディジタルフィルタという離散時間線形時不変システムを通すことで，各周波数の成分がそれぞれどれだけ変化するかが知りたいんだった．フィルタの分野では，周波数応答という用語の代わりに周波数特性と呼ぶことも多い．

やる夫

まあ要するに今までやってきた周波数応答の計算ができればいいんだお? z 変換してから $z = e^{j\omega}$ って置き換えればいいんだったお．

やらない夫

とりあえずはそういうことだな．ただし，いろいろ注意しなくちゃいけないことがあったり，あるいは，周波数特性の計算を省略してもいい場合があったりする．まあその辺は追々話すことにして，具体的な計算を見ていこうか．

まず，インパルス応答が

$$h[n] = \left\{\begin{array}{ll} \frac{1}{5}, & n = 0, 1, 2, 3, 4 0, & \text{otherwise} \end{array}\right.$$ (15.1)

であるようなシステムを考えてみよう．伝達関数はどうなる?

やる夫

えーと，z 変換の公式に入れるお．

$$H(z) = \sum_{k = 0}^{4} \frac{1}{5} z^{-k}$$ (15.2)

$$= \frac{1}{5} \frac{1 - z^{-5}}{1 - z^{-1}}$$ (15.3)

でいいかお．

やらない夫

いいだろう．そのまま先に進んでもいいんだが，練習のために，差分方程式を経由する方法でも伝達関数を求めておこう．システム $h[n]$ の入出力関係を差分方程式で表すとどうなる?

やる夫

んー，入力 $x[n]$ に $h[n]$ をたたみこんで出力 $y[n]$ を得るんだお．だから

$$y[n] = \sum_{k = 0}^{\infty} h[k] x[n - k]$$ (15.4)

$$= \sum_{k = 0}^{4} \frac{1}{5} x[n - k]$$ (15.5)

になるお．

やらない夫

そうだな．これは FIR フィルタだから，線形差分方程式の一般式 (12.15) で $a_1, a_2, \cdots, a_N = 0$ の場合になっていることに注意しておこう．さて，この両辺を z 変換して，伝達関数を求めてみてくれ．

やる夫

1時刻遅延が $z^{-1}$ 倍になるんだから

$$Y(z) = \sum_{k = 0}^{4} \frac{1}{5} z^{-k} X(z)$$ (15.6)

になるお．伝達関数は $H(z) = Y(z) / X(z)$ だから

$$H(z) = \sum_{k = 0}^{4} \frac{1}{5} z^{-k}$$ (15.7)

になるお．結局さっきと同じだお．

やらない夫

もちろん同じになってくれないと困るんだけどな．どっちでやっても計算の構造上全く同じことになっているのをよく確認しておいて欲しい．

さて，伝達関数が求まったから次は周波数応答を求めよう．

やる夫

$z = e^{j\omega}$ を代入するんだお．

$$H(e^{j\omega}) = \frac{1}{5} \frac{1 - e^{-j5\omega}}{1 - e^{-j\omega}}$$ (15.8)

で，えーと，これ前に計算したことある気がするお．

やらない夫

そうだな，つまりは矩形波の離散時間フーリエ変換を計算しているわけなので，窓関数の話のときに計算した通りだ．

やる夫

ああ，そうだったお．式 (11.8) で $N = 5$ にしたものを 1/5 倍すればいいので，

$$H(e^{j\omega})= \frac{1}{5} e^{-j2\omega} \frac{\sin \frac{5\omega}{2}}{\sin \frac{\omega}{2}}$$ (15.9)

になるお．

やらない夫

振幅特性の形状はこうなるな．大雑把に見ると高周波になるほど振幅が小さくなるように作用する．ただし， $\omega = \pm 2\pi/5, \pm 4\pi/5$ で完全に 0 になるのが特徴的だ．

やる夫

つまり， $x[n] = \cos{\frac{2\pi}{5}n}$ とかを入力すると出力は 0 になるってことになるお．

やらない夫

実際に入出力関係をいくつか見てみよう．いずれも，青色○印でプロットされているのが入力 $x[n]$ で，赤色×印が出力 $y[n]$ だ．

$x[n] = \cos{\frac{\pi}{10}n}$ の場合:

$x[n] = \cos{\frac{\pi}{5}n}$ の場合:

$x[n] = \cos{\frac{2\pi}{5}n}$ の場合:

$x[n] = \cos{\frac{3\pi}{5}n}$ の場合:

やる夫

なるほど，確かに振幅特性の通りに出力の振幅が変わってるお．

やらない夫

じゃあ別の例だ．

$$y[n] = 0.9 y[n - 1] - 0.81 y[n - 2] + x[n] + x[n - 2]$$ (15.10)

同じようにやってみようか．

やる夫

まずは z 変換してみるお．

$$Y(z) = 0.9 z^{-1}Y(z) - 0.81 z^{-2} Y(z) + X(z) + z^{-2}X(z)$$ (15.11)

$$Y(z) = \frac{1 + z^{-2}}{1 - 0.9 z^{-1} + 0.81 z^{-2}} X(z)$$ (15.12)

で， $z = e^{j\omega}$ を代入すれば周波数特性がわかるお．

$$H(e^{j\omega}) = \frac{Y(e^{j\omega})}{X(e^{j\omega})}$$ (15.13)

$$= \frac{1 + e^{-j2\omega}}{1 - 0.9 e^{-j\omega} + 0.81 e^{-j2\omega}}$$ (15.14)

…えー，こんなのどんな形のグラフになるのかわからんお．

やらない夫

ちょっと難しいだろうな．計算機を使ってプロットすると，振幅特性はこうなる．

こっちもいくつか入出力関係を見ておこうか．

$x[n] = \cos{\frac{\pi}{10}n}$ の場合:

$x[n] = \cos{\frac{\pi}{3}n}$ の場合:

$x[n] = \cos{\frac{\pi}{2}n}$ の場合:

$x[n] = \cos{\frac{3\pi}{5}n}$ の場合:

やる夫

うーん，こういうややこしい式の場合の特性は，手計算じゃ把握しようがないのかお?

やらない夫

いや，そうとも限らない．周波数特性 $H(e^{j\omega})$ を真面目に計算しなくても，特性のおおまかな形状は目星がつく場合がある．どうするかというと，極と零点の配置を考える．

15.2 極と零点

やらない夫

ラプラス変換の場合と同じで，極ってのは伝達関数が無限大に発散するところで，零点ってのは伝達関数が 0 になるところだ．複素平面でのこれらの配置を考えることで，フィルタの特性についていろいろと知ることができる．

具体例で進めよう．最初の例に戻って極と零点を計算してみるか．

$$H(z) = \frac{1}{5} \frac{1 - z^{-5}}{1 - z^{-1}}$$ (15.15)

やる夫

えーと，ラプラス変換のときと同じように考えればいいんだお? だったら，分子 = 0 という方程式の解が零点で，分母 = 0 の解が極になるお．だから…

やらない夫

ちょっと待った．そのまま進むと見落としが発生するぞ．

やる夫

えっ? どういうことだお?

やらない夫

$z = 0$ のときに分子も分母も無限大になるだろ．ということは，分子 = 0 の解だけが零点とは限らないし，分母 = 0 の解だけが極とは限らなくなってしまう．今の場合は分子の方が $z^{-1}$ の次数が高いから，$z = 0$ で伝達関数は無限大に発散するはずだ．でも $z = 0$ は 分母 = 0 の解ではない．

やる夫

あー，じゃあ，どうすればいいのかお…

やらない夫

分子・分母に $z$ を適当な回数かけて，$z$ の有理式にしてやってから考えるといい．

やる夫

んーと，今の場合は $z^5$ をかけると

$$H(z) = \frac{1}{5} \frac{z^5 - 1}{z^4(z - 1)}$$ (15.16)

になるお．なるほど，これならさっきみたいなことはややこしいことは起きないお．

やらない夫

この辺はちょっと注意しないといけないところだな．部分分数展開のときとかは $z^{-1}$ の有理式として見て式変形したんだった．でも極や零点を考えるときは $z$ の有理式として見ないといけない．

やる夫

じゃあ，まあこれで計算してみるお．まず，分母 = 0 の解は $z = 0$ (4重解) と $z = 1$ だお．これらが極ってことになるお．

それから 分子 = 0 の解は，えーと，$z$ を極座標で表して $z = r e^{j\theta} (r \geq 0, -\pi \leq \theta < \pi)$ とすると

$$r^5 e^{j5\theta} = 1 \cdot e^{j2\pi n} (n: )$$ (15.17)

なので $r = 1$， $\theta = 2\pi n / 5$ ってことになるお．だから，

$$z = 0, \pm e^{j2\pi/5}, \pm e^{j4\pi/5}$$ (15.18)

が零点だお．

やらない夫

複素平面上に極と零点を表示するときは，慣例として極を×印，零点を○印で表すことが多い．原点の極は4重解だということを明示するために 4 と付記しておこう．

やる夫

あれ? $z = 1$ のところはどうなったんだお? 極にも零点にもなっているはずだお．

やらない夫

1 重解の極と零点が同じ場所にあるので，これらは分子・分母で約分されてしまって，極としても零点としてもはたらかない．なので図からも省略してしまっている．実際にシステムの入出力関係に影響を与えるのは，そのほかの極と零点の配置だ．

やる夫

で，その極と零点の配置から何がわかるんだお?

やらない夫

伝達関数と周波数特性の関係をもう一度思い出してみると，単位円 $z = e^{j\omega}$ 上の伝達関数の値を拾っていったものが周波数特性だった．今の極零配置の図で単位円上を見ていくと， $\omega = 0, \pm 2\pi/5, \pm 4\pi/5$ で振幅が 0 になることがわかる．

やる夫

あー，そうか，真面目に計算した振幅特性と一致しているお．…でも，真面目に計算するのとそんなに手間は変わってない気がするお．

やらない夫

じゃあ，2つ目の例だとどうだ?

やる夫

ええと，2つ目は

$$H(z) = \frac{1 + z^{-2}}{1 - 0.9 z^{-1} + 0.81 z^{-2}}$$ (15.19)

だったお．分子・分母を $z^2$ 倍して

$$H(z) = \frac{z^2 + 1}{z^2 - 0.9 z^{1} + 0.81}$$ (15.20)

になるお．だから，極は 分母 = 0 の解で

$$z = 0.9(\frac{1}{2} \pm j\frac{\sqrt{3}}{2}),$$ (15.21)

零点は 分子 = 0 の解で

$$z = \pm j$$ (15.22)

になるお．

やらない夫

この極零配置を見ながら周波数特性を考えてみよう．まず $\omega = \pm \pi/2$ で振幅が 0 になるのはいいだろう．

やる夫

さっきの振幅特性のグラフでも確かにそうなってるお．

やらない夫

次に極の方だが，単位円上にはないがかなり近いところにある．だから，無限大に発散はしないが， $\omega = \pm \pi/3$ で強いピークが出そうだと想像できる．実際に振幅特性のグラフはそうなっているわけだ．

やる夫

んー，振幅特性のグラフを先に見てるからそれで合ってるってわかるけど，「強いピークが出そう」なんて，そんな曖昧なことでいいのかお?

やらない夫

真面目に考えるならこうなる．極と零点を使って伝達関数を表すと

$$H(z) = \frac{(z - j)(z + j)}{(z - 0.9(\frac{1}{2} + j\frac{\sqrt{3}}{2}))(z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2}))}$$ (15.23)

になるな．その絶対値は

$$\vert H(z)\vert = \frac{\vert z - j\vert \vert z + j\vert}{\vert z - 0.9(\frac{... ...rac{\sqrt{3}}{2})\vert \vert z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2})\vert}$$ (15.24)

と表せる．

やる夫

えーと，で，振幅特性は，この $\vert H(z)\vert$ の単位円上の値を拾っていったものになるんだお．

やらない夫

そこまでわかればほとんど答えにたどりついたようなもんだ．$\vert H(z)\vert$ は， $\vert z - j\vert$ と $\vert z + j\vert$ をかけたものを $\vert z - 0.9(\frac{1}{2} + j\frac{\sqrt{3}}{2})\vert$ と $\vert z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2})\vert$ で割ることで計算できる．これらのうち，例えば $\vert z - j\vert$ ってのは複素平面での $z$ と $j$ の間の距離だ．$z$ を単位円上で動かしながらこの長さを考えていくことになる．

やる夫

あー，なるほど，だから例えば $\omega = \pi/3$ のときには分母の $\vert z - 0.9(\frac{1}{2} + j\frac{\sqrt{3}}{2})\vert$ が小さな値になって，それで伝達関数にピークが現れるんだお．

やらない夫

こんな風に，極や零点のそれぞれの影響を分解して考えてやることができる．もちろん正確なグラフを描こうと思うとちょっと無理があるけどな．

やる夫

でも，大まかな傾向をつかむには役に立ちそうだお．

やらない夫

今は振幅特性のみを考えたが，位相特性も同じように把握することができる．$H(z)$ の偏角が

$$\angle{H(z)} = \angle{(z - j)} + \angle{(z + j)} - \angle{(\vert z - 0.9(\frac... ...+ j\frac{\sqrt{3}}{2}))} - \angle{(z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2}))}$$ (15.25)

になるのはわかるか?

やる夫

えーと，例えば $z - j$ は $\vert z - j\vert e^{\angle{(z - j)}}$ って表せるお．絶対値の $\vert z - j\vert$ の部分は位相には関係ないから無視すると， $\angle{H(z)}$ は

$$\angle{H(z)} = \angle \frac{e^{\angle(z - j)} e^{\angle(z + j)}}{e^{\angle(z -... ... j\frac{\sqrt{3}}{2}))} e^{\angle(z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2}))}}$$ (15.26)

$$= \angle e^{\angle(z - j) + \angle(z + j) - \angle(z - 0.9(\frac{... ...} + j\frac{\sqrt{3}}{2})) - \angle(z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2}))}$$ (15.27)

$$= \angle(z - j) + \angle(z + j) - \angle(z - 0.9(\frac{1}{2} + j\frac{\sqrt{3}}{2})) - \angle(z - 0.9(\frac{1}{2} - j\frac{\sqrt{3}}{2}))$$ (15.28)

に確かになるお．

やらない夫

$\angle(z - j)$ ってのは，複素平面で $j$ から $z$ まで引いた線が水平方向となす角度だ．そういう角度を，分子のものについては足し合わせて，分母のものについては引くことで $\angle H(z)$ は計算される． $z$ を単位円上で動かしながらそういう量を考えていくことで，位相特性が得られることになる．

15.3 安定性

やらない夫

ちょっとまた別の例を考えてみよう．

$$y[n] = 1.2 y[n-1] + x[n]$$ (15.29)

こんな差分方程式で表されるフィルタの特性はどうなるか．

やる夫

同じように考えればいいんだお?

$$Y(z) = 1.2 z^{-1}Y(z) + X(z)$$ (15.30)

$$Y(z) = \frac{1}{1 - 1.2 z^{-1}} X(z)$$ (15.31)

だから伝達関数が $H(z) = \frac{1}{1 - 1.2 z^{-1}}$ になるお．周波数特性は $z = e^{j\omega}$ を代入して

$$H(e^{j\omega}) = \frac{1}{1 - 1.2 e^{-j\omega}}$$ (15.32)

になるわけだお．

うーん，グラフ描くの面倒くさそうだから，計算機でプロットしてもらっていいかお?

やらない夫

いまの式の絶対値 $\vert H(e^{j\omega})\vert$ をプロットするとこうなるな．

やる夫

なるほど，直流成分は 5 倍くらいに増幅されて，周波数が増えると急激に増幅率が落ちていって， $\omega = \pi$ 付近では 0.5 倍くらいに落ち着くって感じだお．

やらない夫

じゃあ実際にいくつか入出力関係を見てみるか．

$x[n] = \cos{\frac{\pi}{10}n}$ の場合:

$x[n] = \cos{\frac{\pi}{3}n}$ の場合:

やる夫

えっ，こりゃ一体何なんだお…．盛大に発散してるお．

やらない夫

さっきの差分方程式をよく見てみると，$y$ を計算するのに前の時刻の $y$ の値を 1.2 倍しているわけだ．そんなことを毎時刻やってれば，そりゃ発散するだろ，常識的に考えて…．

やる夫

うーん，でも，ちゃんと z 変換で求めた $H(z)$ に $z = e^{j\omega}$ を代入して周波数特性を計算したはずだお．その計算によると振幅は 0.5 倍からせいぜい 5 倍くらいの間になるはずで，こんな風に発散するなんて聞いてないお．

やらない夫

そこに落とし穴がある．$H(z)$ に $z = e^{j\omega}$ を代入して周波数特性を求めることができるのは，$H(z)$ が単位円 $z = e^{j\omega}$ 上で収束するときだけだ．今は $H(z) = \frac{1}{1 - 1.2 z^{-1}}$ だったわけだが，これは z 変換表によると，時間領域の $h[n] = 1.2^n u_0[n]$ に対応するんだった．一方で，$h[n]$ の z 変換である

$${\cal Z}[h[n]] = \sum_{n = 0}^{\infty} 1.2^n z^{-n}$$ (15.33)

が $\frac{1}{1 - 1.2 z^{-1}}$ に収束するのは $\vert z\vert > 1.2$ のときだけだ．

やる夫

あー，単位円は $\vert z\vert = 1$ だから，収束範囲には入ってないお．

やらない夫

そういう場合は $H(z)$ に $z = e^{j\omega}$ を代入して考えることに意味がない．というか，$h[n]$ の離散時間フーリエ変換がそもそも収束しないわけで，周波数特性を考えること自体が不可能だ．

やる夫

意味のない計算をしていたってことだお…．

やらない夫

そうならないように常に気をつけておく必要があるわけだな．

やる夫

でもどうやって気をつければいいんだお? 何かうまく判断する方法はないのかお?

やらない夫

ああ，実はある．今みたいに出力が発散してしまうような場合，そのフィルタは不安定であるという．不安定なフィルタを避けなければいけない．言い換えれば，フィルタの安定条件がわかれば，それが「うまく判断する方法」だ．

やる夫

なるほど，じゃあ，それを教えて欲しいお．

やらない夫

まあ慌てるな．順番に話していこう．まずは安定性の定義からだ．

やる夫

制御工学でも安定性って概念は出てきたお．それとは違うのかお?

やらない夫

本質的にはほぼ同じことなんだが，違う定義をしてやらなくてはならない．何故かというと，制御ってのは基本的には何らかの量を目標値に近づけようとするものだ．だから目標値からの偏差が時間とともに 0 に十分近づくことを安定性の定義とすればいい．もちろん「十分近づく」の意味がいろいろとありえるのでややこしいんだが，大筋としてはそういうことだ．

やる夫

あー，そうか，信号処理の場合は目標値ってのがないから同じようには考えられないお．

やらない夫

だから，信号処理の場合はこんな風に安定性を定義することが多い．任意の有界な入力，つまり $\vert x[n]\vert < \infty$ に対して，出力も有界であるならば，つまり $\vert y[n]\vert < \infty$ ならば，そのシステムは安定であると定義する．特に区別する場合は，これを有界入力有界出力安定性，略して BIBO (Bounded-Input Bounded-Output) 安定性と呼ぶ．

やる夫

何か急に難しくなったお．

やらない夫

いや，言葉がイカツいだけで，大して難しい話じゃないぞ．つまり，無限大に発散しない入力を入れている限りは，出力も発散しないようにしろってことだ．

やる夫

ああ，まあ，そう言われればわかるお．でも，発散しないあらゆる入力についてテストするなんて現実的には無理だお．

やらない夫

ああ，そりゃそうだが，その点は心配無用だ．実際にあらゆる入力をテストしなくても，インパルス応答が以下の条件を満たせば安定になることを証明できる．

$$\sum_{n = -\infty}^{\infty} \vert h[n]\vert < \infty$$ (15.34)

実はこれは必要十分条件になっている．つまり，安定ならば必ずこの条件が成り立つことも示せる．

やる夫

えーと，インパルス応答の絶対和が発散しないのが条件だってことだお．

やらない夫

まず必要性はほとんど明らかだろう．

やる夫

まあ，絶対和が発散するようなものたたみこんだら出力は簡単に発散しちゃうと思うお．具体的には，えーと， $y[n] = \sum_k h[k] x[n-k]$ なわけだから，例えば $x[-k] = \mathop{\rm sgn}\nolimits h[k]$ になるように入力を選ぶと

$$y[0] = \sum_k h[k] x[0 - k]$$ (15.35)

$$= \sum_k h[k] \mathop{\rm sgn}\nolimits h[k]$$ (15.36)

$$= \sum_k \vert h[k]\vert$$ (15.37)

となって出力が発散するお．だから $\sum_k \vert h[k]\vert$ が発散しないことは BIBO 安定のために必要だお．

やらない夫

十分性も同じようにたたみこみの計算から示せる．入力が有界ということは何か有限な値 $\bar{x}$ があって $\vert x[n]\vert < \bar{x}$ だってことだ．このとき出力の絶対値は

$$\vert y[n]\vert = \vert\sum_k h[k] x[n - k]\vert$$ (15.38)

$$\leq \sum_k \vert h[k]\vert \vert x[n - k]\vert$$ (15.39)

$$\leq \bar{x} \sum_k \vert h[k]\vert$$ (15.40)

$$< \infty$$ (15.41)

になるので有界だ．

やる夫

$\sum_k \vert h[k]\vert$ が発散しなければ，BIBO 安定のためには十分だってことだお．

やらない夫

このことから容易にわかる事実は，FIR フィルタは常に BIBO 安定だということだ．

やる夫

んー，ああ，インパルス応答が有限の長さしか続かないんだから，まあ確かに絶対和は発散しないお．ということは，IIR フィルタだけ注意すればいいんだお．

やらない夫

それからもう一つ．因果的なシステムを考えている場合は， BIBO 安定性を極の配置から判定することができる．

やる夫

ん? どうやって判定するんだお?

やらない夫

すべての極が単位円の内側にあることが，BIBO 安定であることの必要十分条件だ．

今までの具体例を見てみると，最初にやった2つの例は極がすべて単位円内だ．実際フィルタとして正しく動作していた．

しかしついさっきの例は極が $z = 1.2$ にあるので単位円外だ．つまりこのシステムは BIBO 安定ではない．

やる夫

あー，だから出力が発散してしまったんだお．なるほどこれなら簡単だけど，どうしてこんな単純な方法で判断できるんだお?

やらない夫

完全な証明というわけではないが，直観的に説明しておこうか．伝達関数の極ってのがいったい何だったかというと，インパルス応答を

$$h[n] = (w_1 \alpha_1^n + w_2 \alpha_2^n + \cdots + w_n\alpha_N^n)u_0[n]$$ (15.42)

てな具合に指数関数に分解したときの $\alpha_1, \cdots, \alpha_N$ の値，これらが極だ．

やる夫

そうだったお．ただし何か，重解がないこととか分子と分母の次数がどうのこうのとかいろいろ条件があった気がするお．

やらない夫

ああ，完全な証明じゃないといったのはその辺の事情だが，話の大勢に影響がないのでこのケースで続けよう．それらの条件が満たされない場合については 別途議論するが，結論は変わらない．

で，それらの条件が満たされる場合については，既にほとんど結論にたどり着いている．極がすべて単位円内にあるってことは，今の分解されたインパルス応答についてはどういうことになる?

やる夫

$\vert\alpha_i\vert < 1$ なんだから，どの指数関数も時間とともに減衰していって収束するお．

やらない夫

その場合，絶対和 $\sum_{n = 0}^{\infty} \vert\alpha_i^n\vert$は，公比の絶対値が 1 未満の等比数列の和なので，有限の値に落ち着く．インパルス応答はこれらの指数関数を高々 $N$ 個足し合わせたものなので，やっぱり絶対和は有限で押さえられる．

やる夫

逆に，極が 1 個でも単位円外にあったら，その分の指数関数は発散するので，インパルス応答の絶対和は収束しようがないお．

やらない夫

というわけで，極がすべて単位円内にあることが安定性の条件になっていることが，大雑把ではあるけど理解してもらえたんじゃないかと思う．

やる夫

…んー，何か制御工学で習った話と似ているお．制御でも，連続時間システムのインパルス応答をラプラス変換して，その極の配置から安定性を判断できたんだったお．えーと，すべての極の実部が負であることが安定条件だったお．

やらない夫

いいところに気づいた．$z = e^s$ であることを思い出すと，その辺はさらに明確になるぞ．

やる夫

そうだったお．$z = e^s$ で，$s$ 領域の左半平面は $z$ 領域の単位円内に移されるんだったお．さっきの話とぴったり符合しているお．

やらない夫

結局，制御工学で考えた安定性も，これまで説明してきた BIBO 安定性も，本質的なところは同じだってことだ．

15.4 線形位相特性

やる夫

ディジタルフィルタの周波数特性っていいながら，ここまでは振幅特性ばっかり見てきた気がするお．位相の方は気にしなくていいのかお?

やらない夫

もちろん大いに気にする必要あるぞ．そもそも位相特性ってどういう意味だった?

やる夫

ある周波数の単振動を入力したとき，その位相がどのくらいずれて出てくるかを表すものだったお．

やらない夫

そうだな．位相を時間に換算して，入力信号と出力信号の間にどれだけ時間のずれがあるかを表していると考えてもいい．特に因果的なシステムの場合は，システムを通過することによる周波数成分ごとの時間遅れだと考えることができる．因果的な場合，出力が入力より進んでいることはあり得ないからな．

やる夫

ていうか，遅れなんてないに越したことはないわけだお．遅れ時間が 0，つまり位相特性は常に 0 になるようにしよう，ってことじゃダメなのかお?

やらない夫

そういう特性のことを零位相特性という．ある意味理想の特性なわけだが，実現には大きな制約がある．まずそのことを見ていこう．

インパルス応答が $h[n]$ であるシステムを考えよう．インパルス応答を離散時間フーリエ変換したものが周波数特性 $H(\omega)$ だ．このシステムが零位相特性を持つとはどういうことだろう?

やる夫

すべての $\omega$ について $\angle{H(\omega)}$ が 0 だってことだお．だから，$H(\omega)$ は実数値の関数だってことになるお．

やらない夫

よし，じゃあ次の質問だ．$H(\omega)$ が実数値であるためには， $h[n]$ はどういう関数じゃないといけないだろうか?

やる夫

え? 何かつかみどころのない質問だお．

やらない夫

わかりにくければ，時間領域と周波数領域を逆にして類推してみるといい．

やる夫

えー，どういうことだお．逆にするってことは，時間領域で実数値の信号は周波数領域でどうなるかってことかお? …あー，複素指数関数型のフーリエ級数の話のときにそんなことを考えた気がするお．実数値関数 $f(t)$ をフーリエ級数展開すると，正の周波数と負の周波数成分どうしで虚数部分を打ち消しあわないといけないから，振幅は偶関数で，位相は奇関数じゃないといけなかったんだお．

やらない夫

その話を，時間-周波数をそっくりひっくり返して考えるといい．

やる夫

ということは，零位相特性を持つということは，インパルス応答 $h[n]$ の振幅が偶関数で，位相が奇関数だってことになるかお．…って， $h[n]$ の振幅と位相って何だお?

やらない夫

いや，別に，一般には時間信号だって複素数値を取っていいんだから，振幅と位相があっても構わないだろ．まあ多くの応用では実数値の時間信号を考えることが多いわけで，その場合に限定して考えると，零位相特性を持つための条件は $h[n]$ が偶関数であることだということになる．

やる夫

なるほど，インパルス応答が時刻 0 を境に左右対称じゃなきゃいけないわけだお．…ん，つまり $n < 0$ のでも値を持つようなインパルス応答じゃないとダメだってことかお．

やらない夫

そういうことだ．つまり，因果的なシステムは零位相特性を持つことができない．厳密にいうと例外があって， $h[n] = \delta[n]$ の場合，あるいはその定数倍の場合だけは因果的でありながら偶関数にすることができる．まあ何もしないかせいぜい定数倍するだけのフィルタなので，零位相特性なのは当たり前だな．

やる夫

そんなフィルタ面白くも何ともないお．

やらない夫

だから，実質的には，零位相特性を持つ因果的なフィルタは作ることができないと考えてしまっていい．

やる夫

なるほど，そりゃ残念だお．

やらない夫

というわけで何か別のものを考えなければいけないわけだ．多くの場合に次善の策であると考えられているのは，すべての周波数成分が等しい時間遅れを持つような位相特性だ．

やる夫

遅れを 0 にするのはあきらめるわけだお．で，各周波数成分の時間遅れが等しいと何が嬉しいのかお?

やらない夫

波形が歪むのを避けることができる．例えば低域通過フィルタを作ったとしよう．そのフィルタを十分に通過できるような，低周波成分のみからなる信号を入力したとしよう．当然，入力信号にできるだけ近い信号が出力されて欲しいわけだ．

やる夫

そりゃそうだお．低周波の信号はそのまま通してくれなきゃ低域通過フィルタの名がすたるお．

やらない夫

ところが，そのフィルタを通すことで各周波数がバラバラに遅れてしまうと，出力の波形はバラバラに遅れたものの重ね合わせだから，元の波形とは違った歪んだものになってしまう．

やる夫

あー，そっか，時間遅れが一定であれば，どの周波数も同じだけずれてから重ね合わされるので，波形は，遅れはするけど少なくとも歪みはしないんだお．

やらない夫

零位相特性を持つようなインパルス応答 $h_\textnormal{z}[n]$ を $n_\textnormal{d}$ だけシフトして $h_\textnormal{z}[n - n_\textnormal{d}]$ にすれば，$n < 0$ で値を持たなくなって因果的になるとしよう．あらゆる周波数成分に一定の遅れ $n_\textnormal{d}$ を持たせるってことだ．

周波数領域では，時間シフトの公式 (7.3) から，零位相の周波数特性 $H_\textnormal{z}(\omega)$ に $e^{-j\omega n_\textnormal{d}}$ をかけて $H(\omega) = e^{-j\omega n_\textnormal{d}} H_\textnormal{z}(\omega)$ という周波数特性を持つことになる．このような特性のことを，線形位相特性と呼ぶ．

やる夫

んー，何が線形なんだお?

やらない夫

周波数特性の偏角 $\angle{H(\omega)} = \angle\{ e^{-j\omega n_\textnormal{d}} H_\textnormal{z}(\omega) \} = - n_\textnormal{d}\omega$ のことだな．周波数 $\omega$ の成分が，位相としてどれだけずれて出てくるかを表している．これが $\omega$ に対して線形に変化するということだ．

やる夫

なるほど，位相特性をグラフに描くと傾き $-n_\textnormal{d}$ の直線だってことだお．この傾きが，時間領域での遅れ時間を表すんだお．

やらない夫

ちょうど時間領域に言及したところで，線形位相特性の場合のインパルス応答の形状についてもう少し考えておこう．どういう形状になるかというと，さっき見たとおり零位相特性を持つようなインパルス応答を時間シフトするんだった．

やる夫

零位相特性の場合インパルス応答は偶関数で，それを時間シフトして因果的なインパルス応答にしてやるんだったお．

やらない夫

このことからわかるのは，線形位相特性を持つということは，インパルス応答が必ず左右対称な形状をしているということだ．

やる夫

まあそりゃそうだお．対称の軸が原点じゃなくなるだけで，対称な形状はそのままだお．

やらない夫

そしてそこからわかる重要な帰結は，因果的な IIR フィルタは線形位相特性を持ち得ないということだ．

やる夫

えーと，あー，そうか，因果的な IIR フィルタのインパルス応答は時刻 0 より前では値が 0 で，でも時刻が無限大に進んでもえんえんと値が残るんだお．そんなもん左右対称になりようがないお．

15.5 群遅延

やる夫

線形位相特性を持たせることができないってことは， IIR フィルタって使い物にならないのかお?

やらない夫

そんなことはないぞ．全周波数で完全に線形位相でなくても，必要な帯域でだけ近似的に線形位相と見なせれば十分な場合もあったりする．例えば周波数選択フィルタの場合，信号を通過させる帯域だけ考えれば十分なわけだ．

やる夫

ああ，確かに阻止される帯域の出力信号は，波形が歪んでようが大して問題じゃないお．

やらない夫

さらに，応用によっては，個々の周波数成分の遅延はそれほど重要じゃなく，むしろ群遅延と呼ばれるものの方が重要な場合がある．

やる夫

ん? どういうことだお?

やらない夫

具体例で考えよう．前に，振幅変調という概念について説明したと思う．

やる夫

覚えてるお．時間シフトの公式を，時間と周波数を逆にして考えたものだったお．

やらない夫

時間信号 $x[n]$ に $e^{j\omega_\textnormal{c}n}$ をかけると，周波数領域では $X(\omega)$ から $X(\omega - \omega_\textnormal{c})$ に周波数シフトするんだった．

やる夫

例えば音声信号 $x[n]$ を周波数 $\omega_\textnormal{c}$ の搬送波に乗せて通信する場合なんかを表しているんだったお．

やらない夫

変調された後の信号 $x[n] e^{j\omega_\textnormal{c}n}$ には，周波数 $\omega_\textnormal{c}$ やその近くの周波数成分が含まれている．この信号を何らかのフィルタに通すとそれらの周波数成分がそれぞれ遅れて出てくることになるが，その遅延にはあまり興味はない．

やる夫

ん? そうなのかお?

やらない夫

いま送ろうとしている音声信号は，変調信号の振幅の変化として表現されているわけだ．その振幅変化の波形，つまり包絡線がどのくらい遅れるかの方がずっと重要だ．

やる夫

あー，まあ確かにそうかも知れないお．こういう場合に搬送波の位相がどのくらい遅れているかは割とどうでもよさそうだお．

やらない夫

こういう状況は，通信のために変調している場合だけに限らない．例えば 440 Hz の音叉をポーンと鳴らしたときの音声信号を考えよう．このときの信号は，440 Hz のサイン波の振幅，つまり音量が，音叉を叩いた瞬間に急激に立ち上がって，その後ゆっくりと減衰していくような波形になるだろう．

やる夫

あー，そうか，音叉の音量変化の波形が $x[n]$ で，440 Hz のサイン波が搬送波だとすれば同じ話になるお．この場合も，フィルタを通したときに重要なのは包絡線の遅れの方で，440 Hz の振動の位相のずれの方じゃないお．

やらない夫

群遅延ってのは，こういう場合に包絡線の波形がどのくらい遅れるかを表す量だ．結論から言ってしまうと，群遅延は位相特性曲線の傾きで表される．正確にいうと傾きにマイナスをつけたもの，つまり群遅延 $\tau_\textnormal{gd}$ は

$$\tau_\textnormal{gd}(\omega) = - \frac{d}{d\omega} \angle H(\omega)$$ (15.43)

と表すことができる．

ついでに，これまで考えてきた，各周波数成分がそれぞれどのくらい遅れるかという意味での遅れ時間

$$\tau_\textnormal{pd}(\omega) = -\frac{\angle{H(\omega)}}{\omega}$$ (15.44)

のことを位相遅延と呼んだりする．線形位相の場合は，群遅延も位相遅延も等しい定数だ．一般にはいずれも $\omega$ の関数になる．

やる夫

あれ? $\omega$ によって変わるんなら，包絡線波形の遅れ時間ってのはどの $\omega$ のときの値を考えればいいんだお?

やらない夫

ああ，それは，搬送波周波数 $\omega_\textnormal{c}$ のときの値を使うことになる．

やる夫

ふーん，そういうもんなのかお．ていうか何故そうなるのかがわからないから，理解のしようがないお．直観的にでいいから説明してくれお．

やらない夫

…うーん，これに関しては直観的に説明する方がむしろ難しいな．ちゃんと数式を展開する方が説明しやすい．

やる夫

(あー，またヤブヘビをつついた気がするお)

やらない夫

周波数 $\omega_\textnormal{c}$ の搬送波が時間信号 $x[n]$ で振幅変調されているとしよう．さっきから話している通り，そのスペクトルは $X(\omega - \omega_\textnormal{c})$ だ．

やる夫

$X(\omega)$ が $\omega_\textnormal{c}$ だけ周波数シフトされたんだお．

やらない夫

この信号をあるフィルタ $H(\omega)$ に通す．位相特性にだけ注目したいので，振幅特性は常に 1 だとしよう．位相特性は，搬送波周波数 $\omega_\textnormal{c}$ のときに $\angle{H(\omega_\textnormal{c})} = -\phi_\textnormal{c}$ で，その周囲の $X(\omega - \omega_\textnormal{c})$ の成分が存在する範囲では傾き $-\tau_\textnormal{c}$ で一定と見なせるとしよう．

やる夫

えーと， $\omega_\textnormal{c}$ の近辺では線形位相と見なせるってことかお?

やらない夫

いや，そうとは限らない．位相特性の $\omega_\textnormal{c}$ の近辺の直線部分を延長したときに原点を通るならば，確かにその通りだ．ただし今はこの範囲で「傾きが」一定だといっているだけなので，延長線が原点を通る保証は全くない．

やる夫

そっか，なので，この範囲の位相遅延は必ずしも一定というわけではないんだお．

やらない夫

そういうことだ．それにも関わらず，包絡線の波形は歪まずに $\tau_\textnormal{c}$ 遅れるだけになることをこれから示そう．

まず，いま着目している $\omega_\textnormal{c}$ の近くの範囲で，このフィルタの周波数特性は

$$H(\omega) = e^{-j(\phi_\textnormal{c}+ (\omega - \omega_\textnormal{c})\tau_\textnormal{c})}$$ (15.45)

と表せる．

やる夫

いきなり何かややこしいお．えーと，振幅特性が常に 1 ってのは満たしてるお． $\omega = \omega_\textnormal{c}$ のときの位相はというと $H(\omega_\textnormal{c}) = e^{-j\phi_\textnormal{c}}$ になっているのでこれも大丈夫だお．で， $-(\phi_\textnormal{c}+ (\omega - \omega_\textnormal{c})\tau_\textnormal{c})$ を $\omega$ で微分したら $-\tau_\textnormal{c}$ になるので，この近辺で傾きが $-\tau_\textnormal{c}$ だってのも OK だお．合ってそうだお．

やらない夫

変調信号 $x[n] e^{j\omega_\textnormal{c}n}$ をこのフィルタに通す．その出力 $y[n]$ を周波数領域で見ると，

$$Y(\omega) = H(\omega) X(\omega - \omega_\textnormal{c})$$ (15.46)

$$= e^{-j(\phi_\textnormal{c}+ (\omega - \omega_\textnormal{c})\tau_\textnormal{c})} X(\omega - \omega_\textnormal{c})$$ (15.47)

$$= e^{-j\phi_\textnormal{c}} e^{-j(\omega - \omega_\textnormal{c})\tau_\textnormal{c}} X(\omega - \omega_\textnormal{c})$$ (15.48)

になる．

やる夫

うー，面倒くさそうな式だお．

やらない夫

一つずつ見ていけばそうでもないぞ．まず最初の $e^{-j\phi_\textnormal{c}}$ は $\omega$ を含んでない，ただの定数倍だ．残りの部分は複雑そうに見えるが，

$$e^{-j\omega\tau_\textnormal{c}} X(\omega)$$ (15.49)

を $\omega_\textnormal{c}$ だけ周波数シフトしたものだと考えればいい．

やる夫

ふーん，そう思ってみるとそんなに複雑でもなかったお．

やらない夫

さて，これを時間領域で見てみよう．定数倍の $e^{-j\phi_\textnormal{c}}$ はそのままだ． $\omega_\textnormal{c}$ だけ周波数シフトした分は，時間領域では $e^{j\omega_\textnormal{c}n}$ をかけることに対応する．そして， $e^{-j\omega\tau_\textnormal{c}} X(\omega)$ は時間シフトの公式から $x[n - \tau_\textnormal{c}]$ になる．これらをまとめると

$$y[n] = e^{-j\phi_\textnormal{c}} e^{j\omega_\textnormal{c}n} x[n - \tau_\textnormal{c}]$$ (15.50)

という信号がフィルタから出てくることになる．

やる夫

あー，$x[n]$ が $\tau_\textnormal{c}$ だけ遅れてるお．

やらない夫

そこが核心だ．包絡線波形 $x[n]$ を $\tau_\textnormal{c}$ 遅らせたものを搬送波 $e^{-j\phi_\textnormal{c}}e^{j\omega_\textnormal{c}n} = e^{j(\omega_\textnormal{c}n) - \phi_\textnormal{c}}$ に乗せたものが出力されるということだ．搬送波の位相が $\phi_\textnormal{c}$ だけずれているが，そっちには興味はない．包絡線波形が群遅延 $\tau_\textnormal{c}= \tau_\textnormal{gd}(\omega_\textnormal{c})$ の分だけ遅れることが確認できたわけだ．

やる夫

…あれ? ちょっと待って欲しいお．$x[n]$ って離散時間信号なんだお? $\tau_\textnormal{c}$ が整数じゃないときは， $x[n - \tau_\textnormal{c}]$ ってどうなるんだお?

やらない夫

鋭いな．その通り， $x[n - \tau_\textnormal{c}]$ という表記は $\tau_\textnormal{c}$ が整数のときしか意味を持たない．が，一般に群遅延は実数なので困ったことになる．

詳しい説明は省くが，こんな風に理解して欲しい．離散時間信号 $x[n]$ を，サンプリング定理のときに出てきた理想的低域通過フィルタに通して，いったん連続時間信号 $\tilde{x}(t)$ を作る．次に，これを $\tau_\textnormal{c}$ だけ遅らせた $\tilde{x}(t - \tau_\textnormal{c})$ を作る．最後にこれをサンプリングし直す．このようにして得られる離散時間信号のことを $x[n - \tau_\textnormal{c}]$ という記法で表すと約束することにしよう．

やる夫

んー，つまり $x[n]$ を，何か背後にある連続時間信号をサンプリングしたものだと考えて，遅延させる必要があるときはその連続時間信号の方を遅延させるって話になるお．天下りの説明だからちゃんとは理解できないけど，辻褄は合ってそうな気がするお．

やらない夫

ちゃんと理解したかったら，今のような手順で遅延させた信号のスペクトルと，元の信号のスペクトルに $e^{-j\omega\tau_\textnormal{c}}$ をかけたものとが一致することを計算で確かめてみるといい．

さて，また長くなりそうなので，この辺で切り上げてまとめよう．

• ディジタルフィルタの周波数選択特性は，周波数特性 (= 周波数応答) のうち特に振幅特性 (= 振幅応答) を調べることで把握できる．
  • 周波数 $\omega$ の単振動は，振幅が $\vert H(\omega)\vert$ 倍されて出力される．
• 位相特性 (= 位相応答) を調べることで，入力信号がどのくらい遅れて出てくるかを把握できる．
  • $-\angle H(\omega) / \omega$ を位相遅延といい，周波数 $\omega$ の単振動がどのくらい遅れて出力されるかを表す．
  • $- d\angle H(\omega) / d\omega$ を群遅延といい，周波数 $\omega$ の搬送波を振幅変調した際の包絡線波形がどのくらい遅れて出力されるかを表す．
  • 位相遅延 = 群遅延 = 一定 であるような特性を線形位相特性と呼ぶ．線形位相特性を持つにはインパルス応答が左右対称な形状をしていなくてはならず，したがって因果的な IIR フィルタは線形位相特性を持ち得ない．
• 伝達関数の極と零点の配置を考えることで，フィルタの大まかな特性を把握できる．
  • $z$ 領域の単位円の近くに零点があれば，振幅特性はその偏角の周波数の付近で小さくなる．特に単位円上に零点があるなら，その周波数で振幅特性が 0 になる．
  • $z$ 領域の単位円の近くに極があれば，振幅特性はその偏角の周波数の付近で大きくなる．
• 入力が有界である限り出力も有界でことが保証されるフィルタは BIBO 安定であるという．
  • FIR フィルタは常に安定である．
  • 因果的なフィルタが BIBO 安定であるための必要十分条件は，すべての極が単位円の内側にあることである．
  • 不安定なフィルタの伝達関数 $H(z)$ に $z = e^{j\omega}$ を代入して形式的に $H(e^{j\omega})$ を求めても，周波数応答にはならない．