16. ディジタルフィルタの設計

16.1 フィルタの仕様と設計方針

やらない夫

ディジタルフィルタについていろいろと考えてきたわけだが，要するに，フィルタの伝達関数が与えられれば周波数特性とか安定性とかいった特性を知ることができるってことだった．

で，いよいよその逆の話だ．望ましい周波数特性が与えられたときに，そのような特性を実現するフィルタはどうやって作ればいいかという問題を考える．フィルタの設計問題ってやつだな．

やる夫

んー，そんな改まって考えなきゃいけない話かお? 周波数特性 $H(\omega)$ が与えられてるんだから，離散時間逆フーリエ変換すればインパルス応答 $h[n]$ がわかるお．で，それをたたみこむ $y[n] = h[n] * x[n]$ ってフィルタを作ればそれで終わりだお．簡単な話だお?

やらない夫

そう思うなら実際に具体例を計算してみようか．最も基本的なものってことで，理想的低域通過フィルタ

$$G_{-\omega_\textnormal{c},\omega_\textnormal{c}}(\omega) = \left\{ \begin{array}{ll} 1, & \vert\omega\vert \leq \omega_\textnormal{c} 0, & \text{otherwise} \end{array}\right.$$ (16.1)

を考えよう．周波数 $\omega_\textnormal{c}$ 以下の低域は完全に通すが，それより上の周波数は完全に遮断する．周波数特性が常に実数値を取っているので零位相特性になっている．

$\omega_\textnormal{c}$ のことを遮断周波数とか，カットオフ周波数とか呼んだりする．ディジタルフィルタを考えているので，角周波数としては $-\pi$ 〜 $\pi$ の範囲しか考えないことに注意しておこう．遮断周波数は $\omega_\textnormal{c}< \pi$ の値のみ取れるってことだな．

やる夫

離散時間逆フーリエ変換してみるお．

$${\text{DTFT}^{-1}}[G_{-\omega_\textnormal{c},\omega_\textnormal{c}}(\omega)] = \frac{1}{2\pi} \int_{-\pi}^{\pi} G_{-\omega_\textnormal{c},\omega_\textnormal{c}}(\omega) e^{j\omega n} d\omega$$ (16.2)

$$= \frac{1}{2\pi} \int_{-\omega_\textnormal{c}}^{\omega_\textnormal{c}} e^{j\omega n} d\omega$$ (16.3)

ってなるから，あ，なんだ，計算上は式 (3.31) で計算した逆フーリエ変換と一緒だお．あのときの $a$ が $\omega_\textnormal{c}$ になって， $t$ は $n$ になるだけだから結局

$${\text{DTFT}^{-1}}[G_{-\omega_\textnormal{c},\omega_\textnormal{c}}(\omega)] = \frac{1}{\pi n} \sin{\omega_\textnormal{c}n}$$ (16.4)

になるお．

やらない夫

sinc 関数を離散化したものだな．グラフはこんな形だった．本当は離散的に描かないといけないが，面倒なので連続的に描いている．以下同様だ．

さて考えてほしいんだが，これをどうやって実際にたたみこむ?

やる夫

どうやってって，そんなもん公式通りたたみこめばいいんだお．

$$y[n] = \left\{\frac{1}{\pi n} \sin{\omega_\textnormal{c}n}\right\} * x[n]$$ (16.5)

$$= \sum_{m = -\infty}^{\infty} \left\{\frac{1}{\pi m} \sin{\omega_\textnormal{c}m}\right\} x[n - m]$$ (16.6)

…あ，そうか，無限に総和とらなきゃならないお．

やらない夫

その通り．だから有限の時間では計算できない．ついでに言っておくと，sinc 関数は時刻が負のときも値を持っている．つまりこのフィルタは 因果的ではない．

やる夫

あー，てことは，仮に無限の計算能力を持っていたとしても，リアルタイムに動作するフィルタとしては動かせないんだお．ていうか，よく考えるとそもそも零位相特性の因果的フィルタは作れないんだったから当たり前だお．

やらない夫

といった辺りが最初の疑問に対する答えだ．周波数特性を単純に離散時間逆フーリエ変換しても，物理的に実行可能なフィルタが得られるとは限らないし，まあ普通はそううまくいかない．

やる夫

じゃあどうすればいいんだお?

やらない夫

それが今回の話題なわけだが，具体的な方法は，それこそ数え切れないほど提案されている．そのうち本当に基本的なものだけを簡単に説明していこうと思う．

やる夫

やらない夫の言う「簡単に」は信用できないお．

やらない夫

…ああ，確かにそうかも知れん．

やる夫

否定して欲しかったお．

やらない夫

基本的なものといっても，いくつかの方針に分かれる．最初の，そしておそらく最大の分かれ目は，FIR にするか IIR にするかだ．

やる夫

えーと，インパルス応答が有限で終わるか，無限に続くかだったお．

やらない夫

さっきの理想的低域通過フィルタの例でいくと，本来は無限に続くインパルス応答を何とか誤魔化して有限で打ち切ると FIR フィルタになる．結果として得られるフィルタ特性は，もちろん本来欲しかったものとは多かれ少なかれ異なったものになる．

やる夫

同じ打ち切るんでも，できるなら長く取る方が本来の特性には近づけられそうだお．

やらない夫

基本的にはその通りだ．FIR フィルタってのは時間領域の差分方程式で

$$y[n] = \sum_{m=0}^{M} h[m]x[n - m],$$ (12.4)

あるいは伝達関数で表すなら

$$H(z) = \sum_{n=0}^{M} h[n]z^{-n}$$ (16.7)

と書けたわけだが，この定数 $M$ をフィルタの次数という．次数が大きいほど打ち切りの影響は小さくできる．しかしその分，過去の入力値 $x[n - m]$ をたくさん記憶しておかないといけなくなるし，もちろん計算量も増える．それに，初期状態から始めてフィルタが定常的な状態になるまでに長い時間がかかるようになる．

やる夫

ていうことは，できるだけ次数は小さく取りたいわけだお．でもそうすると特性が要求仕様からずれてってしまうお．

やらない夫

なので，特性をどのくらいまでなら妥協できるかを明確にしておくことが重要だ．妥協できる範囲でできるだけ小さな次数を選ぶ．

やる夫

インパルス応答を有限で打ち切るから妥協が必要なんだお? じゃあ，IIR フィルタならインパルス応答を打ち切らなくていいんじゃないかお?

やらない夫

打ち切る必要はないが，だからといって所望のインパルス応答を必ず表現できるわけじゃない．IIR フィルタは差分方程式で

$$\sum_{k = 0}^{N} a_k y[n - k] = \sum_{k = 0}^{M} b_k x[n - k]$$ (12.15)

あるいは伝達関数で表して

$$H(z) = \frac{\sum_{k = 0}^{M} b_k z^{-k}}{\sum_{k = 0}^{N} a_k z^{-k}}$$ (16.8)

と書けるわけだが，$N$ や $M$ が有限じゃないとやっぱり実行できない．

やる夫

あー，IIR フィルタのインパルス応答が無限に続くからといって，無限に続くどんなインパルス応答でも実現できるというわけじゃないんだお．

やらない夫

次数が小さいとやっぱり表現できるインパルス応答の自由度が小さいので，仕様として与えられた特性そのものを実現するのは難しくなる．その辺の事情は FIR フィルタと同じだ．

やる夫

やっぱり妥協が必要だってことだお．でも，具体的にどんな点を妥協をすればいいんだお?

やらない夫

そうだな，ざっくり言うと2点ある．1つめは遮断特性だ．理想的なフィルタ特性では，遮断周波数の前後で「完全に通過」「完全に阻止」としたいわけだが，後で見るように実際には「通過域」から「阻止域」まで徐々に振幅特性が落ちていく区間が発生してしまう．フィルタの次数が小さいと，この区間を小さくするのが難しい．遮断特性を急峻にしにくいという言い方をする．

やる夫

ふーん，つまり遮断特性が鈍るのをどのくらいまで許容できるかってのが第 1 の妥協点になるわけだお．

やらない夫

2つめはリップル特性と呼ばれるものだ．理想的には，通過域では振幅を完全に 1 に，阻止域では完全に 0 にしたい．しかし実際には，これも後で実例を見るが，振幅が波打ってしまったりする．この振幅の波のことをリップルと呼んでいる．あるいは，リップルの振幅の大きさのことをリップルと呼ぶ場合もある．

やる夫

リップルも小さい方が嬉しいけど，どのくらいまで発生してもいいかってのが妥協点になるわけだお．

やらない夫

フィルタの仕様に関する用語をまとめておこう．低域通過フィルタの場合だけ説明するが，ほかのフィルタでも同じように使う．

通過域と阻止域それぞれで許容できるリップルの大きさを，それぞれ許容通過域リップル，許容阻止域リップルと呼ぶ．で，今まで通過域とか阻止域とかをちゃんと定義せずに使ってきたが，振幅の変動が許容通過域リップル内に収まっている帯域が通過域だ．同じように，許容阻止域リップル内に収まっている帯域を阻止域と呼ぶことにする．残った部分が遷移域ということになる．

やる夫

遷移域で徐々に振幅が落ちていくと，遮断周波数って曖昧にならないかお?

やらない夫

曖昧だと困るので，普通は「振幅の2乗が半分になる周波数」と定めることが多い．「振幅が $1/\sqrt{2}$ 倍になる周波数」と言い換えてもいい．

やる夫

つまり，こういう遮断特性とかリップル特性とかをどのくらい妥協できるかを考えて，その範囲でできるだけ次数の小さいフィルタを作ろうってことだお．

やらない夫

そういうことになるな．もちろん，FIR なのか IIR なのかとか，どういう設計法を使うのかとかによってそれらの特性への影響も変わってくるので，それぞれのケースでよく考える必要がある．

やる夫

んー，そういえば，FIR にするか IIR にするかが大きな分かれ目だっていう話の途中だったお．いまだにその決め手がわからないままだお．

やらない夫

後で実例を見たときに確認して欲しいんだが，大雑把に言うと，同じ次数なら IIR フィルタの方が FIR フィルタよりも急峻な遮断特性を得やすい．なので，その点では IIR フィルタの方が有利だ．

やる夫

じゃあ FIR フィルタなんて必要ないのかお?

やらない夫

そんなことないぞ．前回の話を思い出して欲しいんだが，IIR フィルタには線形位相特性を持たせることが絶対にできない．だから，線形位相特性が必要なときは FIR フィルタを選ばざるを得ない．

やる夫

あー，そうか，波形が歪まないことが重要な場合は FIR フィルタを使うことになるわけだお．

やらない夫

そのほかにも，FIR フィルタが常に安定だという点が設計上の利便性につながる場合があったりとか，数値計算上の誤算の影響も FIR の方が受けにくかったりとか，いろいろな事情が FIR フィルタを選ぶ理由になる．まあその辺の詳しい話は教科書で勉強してもらうとして，以降では FIR フィルタ，IIR フィルタのごく基本的な設計法を紹介しようと思う．

16.2 FIR フィルタの窓関数設計

やらない夫

じゃあまず FIR フィルタからだ．要求仕様として与えられた周波数特性 $H_\textnormal{d}(\omega)$ を離散時間逆フーリエ変換したときに，仮に有限長のインパルス応答が得られて，しかもその長さが十分に短いなら，そのままたたみこめば FIR フィルタの出来上がりだ．しかし一般にそうはならないというのはさっきの例で見た通りだ．

やる夫

そうだったお．理想的低域通過フィルタはインパルス応答が sinc 関数で無限に続いたお．

やらない夫

例としてはちょうど良いので，その理想的低域通過フィルタの場合で考えていこう．ともかくインパルス応答を何とか有限長で打ち切ってやらなくてはならないわけだ．さてどうするかだが，まずは安直に時刻 $\vert n\vert \leq L$ の範囲だけ残して，それ以外はばっさり捨ててみよう．

やる夫

本当に安直だお．ていうかそれじゃそもそも因果的なフィルタにすらできてないお．

やらない夫

因果性のことはしばらく置いておこう．ともかくこれで，有限長のインパルス応答が作れて，非因果的ではあるが FIR フィルタとして実行できる．それから，このインパルス応答は偶関数なので，零位相特性になっていることにも注意しておこう．

やる夫

ああ， $\vert n\vert \leq L$ ってなふうに原点を中心とした区間を切り出したのはそれが理由かお．

やらない夫

さて，こんな風に安直に作った FIR フィルタは，実はある意味で最適なものだってことを示すことができる．

やる夫

最適? ある意味で? 何のことだお?

やらない夫

理想的な周波数応答 $H_\textnormal{d}(\omega)$ と，実現できる周波数応答 $H(\omega)$ の平均2乗誤差

$$J = \int_{-\pi}^{\pi} \vert H_\textnormal{d}(\omega) - H(\omega)\vert^2 d\omega$$ (16.9)

を考えよう．$H(\omega)$ をいろいろと変化させたときに $J$ が最小になるもの，これが最適な $H(\omega)$ だ．

やる夫

えー， $H(\omega) = H_\textnormal{d}(\omega)$ で最小値 $J = 0$ になるんじゃないかお?

やらない夫

いやいや，そういうことじゃなくて，あくまでインパルス応答が $\vert n\vert > L$ で 0 になるようなもののみを考えて，そのうちで $J$ が最小になるってことだ．

$$h[n] \stackrel{\text{DTFT}}{\longrightarrow}H(\omega) \text{where} h[n] = 0 (\vert n\vert > L)$$ (16.10)

を満たす $H(\omega)$，あるいは $h[n]$ の範囲で考えるってことだな．

やる夫

あー，そういうことかお．

やらない夫

理想的な周波数応答を離散時間逆フーリエ変換したものを $h_\textnormal{d}[n]$ と書くことにしよう．

$$h_\textnormal{d}[n] \stackrel{\text{DTFT}}{\longrightarrow}H_\textnormal{d}(\omega)$$ (16.11)

すると，

$$h_\textnormal{d}[n] - h[n] \stackrel{\text{DTFT}}{\longrightarrow}H_\textnormal{d}(\omega) - H(\omega)$$ (16.12)

の関係があるのはわかるか?

やる夫

えーと，離散時間フーリエ変換の線形性から明らかだお．

やらない夫

ここでパーセバルの等式 (9.3) から

$$\sum_{n=-\infty}^{\infty} \vert h_\textnormal{d}[n] - h[n] \vert^... ...i} \int_{-\pi}^{\pi} \vert H_\textnormal{d}(\omega) - H(\omega) \vert^2 d\omega$$ (16.13)

であることを使うと，

$$J = 2\pi \sum_{n = -\infty}^{\infty} \vert h_\textnormal{d}[n] - h[n] \vert^2$$ (16.14)

$$= 2\pi\sum_{\vert n\vert \leq L} \vert h_\textnormal{d}[n] - h[n] \vert^2 + 2\pi\sum_{\vert n\vert > L} \vert h_\textnormal{d}[n] - h[n] \vert^2$$ (16.15)

$$= 2\pi\sum_{\vert n\vert \leq L} \vert h_\textnormal{d}[n] - h[n] \vert^2 + 2\pi\sum_{\vert n\vert > L} \vert h_\textnormal{d}[n] \vert^2$$ (16.16)

と書ける．

やる夫

えーと，パーセバルの等式で周波数特性の平均2乗誤差をインパルス応答の平均2乗誤差に置き換えてるんだお．で，えーと，あ，総和の範囲を $\vert n\vert \leq L$ とそれ以外に分けたんだお．

やらない夫

$\vert n\vert > L$ についての総和の部分は定数なので， $\vert n\vert \leq L$ についての総和の部分が最小になるときに $J$ が最小ってことになるな．

やる夫

あー，だから， $\vert n\vert \leq L$ の範囲で $h[n] = h_\textnormal{d}[n]$ とすると $J$ が最小になるわけだお．

やらない夫

それってのはつまり，さっき安直に作ったインパルス応答と同じものだ．$\vert n\vert > L$ の範囲で 0 になるようなインパルス応答を持つあらゆるフィルタの中で，与えられた周波数特性の仕様からの平均2乗誤差を最小にするものが，この安直な手順で得られることになる．

やる夫

あ，そういえばそうだお．じゃあこれで FIR フィルタの設計は終わりかお?

やらない夫

いや，非因果的なままでいいのであればこれでいいんだが，因果的なフィルタが欲しい場合はまだ続きがある．

やる夫

そっか，まだ因果的になってなかったお．

やらない夫

でも後はもう簡単だ． $L$ 時刻だけ時間シフトしてやればいい．前回やったように，これで振幅特性は変わらず，線形位相特性を持つフィルタが作れる．インパルス応答が $n = 0, \cdots 2L$ の間続くので，$2L$ 次の FIR フィルタだ．

やる夫

なるほど，それなら因果的になるお．

やらない夫

というわけで，ようやく実際にリアルタイムに走らせることができる FIR フィルタが設計できたわけだ．元の $H_\textnormal{d}(\omega)$ として他の形状のものを与えても，それが零位相であれば，今のやり方で線形位相特性を持つフィルタを設計できる．

やる夫

実際，インパルス応答を途中で打ち切った影響ってどんなもんなのかお?

やらない夫

具体例を見てみようか．遮断周波数 $\omega_\textnormal{c}= 0.4\pi$ の低域通過フィルタを，次数 $M = 2L$ を変えながら設計してみた例だ．

やる夫

なるほど，次数が小さいと遮断特性が鈍ってるお．あと，リップルが結構大きいんだお．

やらない夫

そうだな．リップルが大きくて，しかも次数を大きくしても消える気配がない．

やる夫

ちょっと気持ち悪いお．

やらない夫

少しその原因を考えておこう．$\vert n\vert > L$ の部分を単純に打ち切るっていうことは，インパルス応答に矩形窓をかけているということだ．

やる夫

ああ，周波数分析のときに考えた話だお．

やらない夫

あのときと同じように，時間領域で矩形窓をかけることの意味を周波数領域で考えてみよう．

やる夫

えーと，式 (11.8) で計算したように，長さ $N$ の矩形窓の振幅スペクトルは $\sin \frac{\omega N}{2} / \sin \frac{\omega}{2}$ になるんだったお．sinc 関数っぽい形状だったお．

やらない夫

時間領域で矩形窓を乗算するんだから，周波数領域ではそのスペクトルをたたみこむことになるわけだ．メインローブの影響で，近い周波数どうしが平滑化される効果が出てくる．そのせいで遮断周波数前後の振幅が平滑化されるので，遷移域がなだらかになってしまう．次数が大きいほうがメインローブが細くなるので，遷移域はより急峻になる．これが遮断特性に対する説明だ．

やる夫

なるほど，同じように考えて，サイドローブが影響した結果がリップルなんだお．

やらない夫

その通りだ．で，矩形窓のスペクトルはサイドローブが大きくて，しかも窓長を大きくしても消えてくれるわけじゃない．リップルがどんどん細かくはなるが，特に遮断周波数の前後ではリップルの高さがいつまでたっても消えてくれない．このことをギブスの現象と呼んだりする．数学の授業でフーリエ級数を習ったときに聞いたかもしれない．

やる夫

…聞いたかもしれないけど，聞いてないかもしれないお．

やらない夫

まあともかく，矩形窓で打ち切っている以上，次数をどんなに増やしてもリップルの発生からは逃れられないということだ．じゃあ，リップルの発生を避けたい場合にはどうしたらいいだろう?

やる夫

話が見えてきたお．周波数分析のときと同じように，別の窓関数を使うのかお?

やらない夫

そういうことになるな．ハミング窓を使った設計例はこうなる．

やる夫

あー，リップルがほとんど目立たないお．

やらない夫

その代わり，遮断特性は矩形窓の場合よりも鈍っていることにも注意して欲しい．矩形窓と比べて，ハミング窓はサイドローブが小さい代わりにメインローブは太くなるんだった．リップルを抑えるために遮断特性を犠牲にしていると考えるといいかな．

やる夫

結局，目的に応じて使い分ける必要があるわけだお．

やらない夫

このような設計方法を総称して，FIR フィルタの窓関数設計と呼ぶ．所望の遮断特性やリップル特性が得られるようにいろいろと試行錯誤しないといけないのが厄介なところだが，FIR フィルタの設計は他の方法でもだいたい試行錯誤が必要だ．

16.3 IIR フィルタの間接設計

やる夫

じゃあ，IIR フィルタの方はどうするのかお?

やらない夫

FIR の場合のようにインパルス応答を定めてたたみこめばいいってわけにはいかない．インパルス応答が無限に続くわけだからな．あくまで，差分方程式を定めてやる必要がある．言い換えると伝達関数を $z$ の有理多項式で表してやる必要があるってことだ．

やる夫

なんか難しそうだお．

やらない夫

これまたたくさんの方法があるわけだが，一番基本的なものとして挙げられるのは，アナログフィルタから変換することでディジタルフィルタを作る方法だろう．間接設計と呼ばれる．

やる夫

いったんアナログフィルタの設計を経由するから間接ってことかお．

やらない夫

そうだな．アナログフィルタに関しては，古くから電気回路の分野で詳細に検討されている．こういう特性のフィルタが作りたければこうやって作る，みたいな話がかなり確立されているわけだ．その結果を利用してディジタルフィルタを作ってやろうということだな．

やる夫

アナログフィルタってどんなものなのかお?

やらない夫

代表的なものとして，バタワースフィルタを挙げておこう．遮断周波数 $\Omega_\textnormal{c}$ の $N$ 次バタワース低域通過フィルタは，伝達関数が

$$H_\textnormal{bw}(s) = \frac{1}{\prod_{k=0}^{N-1} (s - s_k)}$$ (16.17)

$$s_k = \Omega_\textnormal{c}e^{j\left\{\frac{\pi}{2} + \frac{(2k + 1)\pi}{2N}\right\}} (k = 0, \cdots ,N - 1)$$ (16.18)

と書ける．その振幅特性は

$$\vert H_\textnormal{bw}(j\Omega)\vert = \frac{1}{\sqrt{1 + \left(\frac{\Omega}{\Omega_\textnormal{c}}\right)^{2N}}}$$ (16.19)

になることが知られている．

やる夫

こりゃまた盛大な天下り式だお．なんだかさっぱりわからんお．

やらない夫

アナログフィルタが本題ではないからな．もし興味があるなら付録を参照して欲しい．伝達関数の式の解釈をざっとだけ説明しておくと， $s$ 領域の原点を中心とした半径 $\Omega_\textnormal{c}$ の円周上の左半分を $N$ 等分して，その各区間の真ん中に極を配置する．これが $N$ 次バタワースフィルタの極配置だ．極が左半面にしかないので，安定であることに注意しておこう．

やる夫

ん? あー，そうか，複素平面に円が描かれるとつい $z$ 領域かと思ってしまうけど，アナログフィルタだから z 変換じゃなくてラプラス変換で表されているんだお．だから左半平面に極があれば安定なんだお．

やらない夫

振幅特性はグラフで見ておこう．FIR の例と同じく $\Omega_\textnormal{c}= 0.4\pi$ としたものだ．その辺りが遮断周波数になっていて，次数が増えるごとに遮断特性が急峻になっていくのがわかってもらえればいい．あと，振動的なリップルが全くないという点がバタワースフィルタの特徴的なところだ．

やる夫

でも結局天下りなのでよくわからないのは変わらんお．

やらない夫

まあそうなんだが，重要なのはこういう振幅特性を持つ低域通過フィルタの伝達関数が，$s$ の有理式として表せているということだ．

バタワースフィルタのほかにも，チェビシェフフィルタとか楕円フィルタとかいろいろなものがあって，それぞれ異なる特徴の特性を持っていて，それぞれ $s$ の式で表すことができる．そういう $s$ の式を，$z$ の式にうまいこと変換しようという話の流れになる．一般に，変換の元となるアナログフィルタのことをプロトタイプフィルタとか呼んだりする．プロトタイプってのは「原型」ってことだな．

やる夫

要求仕様に合わせてよさげなフィルタのよさげな次数を選べばいいわけだお．

やらない夫

そうなんだが，FIR フィルタのときとはちょっと違う点がある．プロトタイプフィルタとして使うようなアナログフィルタは，要求仕様から必要な次数を決めるときに，数式一発で求められるものが多い．具体例は 付録 で見てもらうことにしよう．まあコンピュータが発達していなかった時代から使われてきたものなので，そうじゃないと使い物にならなかったんだろうな．基本的に試行錯誤が必要な FIR フィルタ設計とはその点で大きく違うってことを知っておいて欲しい．

16.3.1 インパルス不変変換

やる夫

で，その $s$ から $z$ への変換ってのはどうやるんだお?

やらない夫

一番基本的なのは，プロトタイプフィルタのインパルス応答をサンプリングして離散時間化するという考え方だ．

やる夫

サンプリングすれば，離散時間のインパルス応答が得られて，それをたたみこめばディジタルフィルタになるってことかお? でも，IIR フィルタなんだからインパルス応答は無限に続くんだお．有限の手順でサンプリングできないお．

やらない夫

そう，だから，インパルス応答を実際に数値的にサンプリングするってのは不可能だ．そうじゃなくて，インパルス応答をサンプリングした場合に得られるであろう離散時間の伝達関数 $H(z)$ を，プロトタイプフィルタの伝達関数から，数式の変換で求めてやることを考える．

やる夫

えーと，何が何だかわかりませんお．

やらない夫

バタワースフィルタをプロトタイプにする場合を考えよう．プロトタイプフィルタの伝達関数 (16.19) を部分分数展開すると，

$$H_\textnormal{bw}(s) = \sum_{k = 0}^{N-1} \frac{w_k}{s - s_k}$$ (16.20)

な形に書ける．だからインパルス応答は，ラプラス変換表から

$$h_\textnormal{bw}(t) = u_0(t) \sum_{k = 0}^{N-1} w_k e^{s_k t}$$ (16.21)

と書けることになる．

やる夫

指数関数の和になるんだお．ラプラス変換の話のときに聞いたことだお．

やらない夫

このそれぞれの指数関数を，サンプリング周期 1 でサンプリングして離散化することを考える．つまり，

$$h[n] = u_0[n] \sum_{k = 0}^{N-1} w_k e^{s_k n}$$ (16.22)

$$= u_0[n] \sum_{k = 0}^{N-1} w_k (e^{s_k})^n$$ (16.23)

というインパルス応答を考えることになる．これを z 変換するとどうなる?

やる夫

えーと，これなら z 変換表を見ればいいんだお．

$$H(z) = \sum_{k = 0}^{N-1} \frac{w_k}{1 - e^{s_k} z^{-1}}$$ (16.24)

になるお．

やらない夫

これでインパルス応答をサンプリングできたわけだが，結果的には，プロトタイプフィルタの伝達関数を部分分数展開した式 (16.22) に

$$\frac{1}{s - s_k} \longrightarrow \frac{1}{1 - e^{s_k}z^{-1}}$$ (16.25)

という置き換えを機械的に施しただけだ．

やる夫

あー，確かにそうだお．

やらない夫

つまり，実際に時間領域でサンプリングしなくても，部分分数展開された伝達関数に式 (16.26) の変換をするだけで，ディジタルフィルタ $H(z)$ が得られる．インパルス不変変換と呼ばれる方法だ．

やる夫

なるほど．今はバタワースフィルタをプロトタイプにして低域通過フィルタを設計したけど，この手順さえ覚えておけば，あとは要求仕様に合わせてその都度よさげなプロトタイプフィルタを選んで来れば，いつでもディジタルフィルタが作れるってわけだお．

やらない夫

…うーん，常にそう言い切れないのが残念なところだ．

やる夫

えっ，何かまずいのかお?

やらない夫

ちょっと注意しなくてはならない点がある．インパルス不変変換は，インパルス応答をサンプリングして離散時間化するものだ．なので，サンプリング特有の問題が発生し得る．

やる夫

えー，サンプリング特有の問題って何だお．

やらない夫

サンプリング定理の話をしたときのことを思い出してくれ．サンプリングする際に注意しなくちゃいけないことがなかったか?

やる夫

注意? えーと，サンプリングの話のときだから…，ナイキスト何とかを超えるとエイリアシングがどうのこうのってやつかお?

やらない夫

どうのこうのってのはどうかと思うが，ポイントはそこだ．サンプリングされる連続時間信号がナイキスト周波数以上の周波数成分を持っていると，エイリアシングが起きてしまうんだった．

やる夫

そうそう，それが言いたかったお．

やらない夫

つまり，インパルス応答がサンプリング周波数の半分までに帯域制限されていないと，プロトタイプフィルタの周波数特性とは違った特性のフィルタができてしまう．

やる夫

どんな風に違っちゃうのかお?

やらない夫

周波数領域でどうなるかちょっと見てみよう．プロトタイプフィルタのインパルス応答を周期 1 でサンプリングするわけだから，周波数スペクトルは，元のスペクトルを $2\pi$ ごとに繰り返して足し合わせたものになる．

やる夫

だから，周波数 $\pi$ 以上の成分があると，スペクトルの裾どうしが重なってエイリアシングが起きるんだお．

やらない夫

さて，例えばさっきのバタワースフィルタはどうだった?

やる夫

うーん，式 (16.20) の上では $\Omega$ をどんなに増やしても周波数成分は 0 にはならないお．でもまあ，グラフを見ると次数によってはほとんど 0 と見なせそうだったりそうでもなさそうだったり，いろいろだお．

やらない夫

厳密に言うと常にエイリアシングが起きているわけだな．それでも影響が十分小さければまあそれほど問題にならない．しかし，次数が低い場合や，遮断周波数がナイキスト周波数に近い場合は影響が大きく出てしまう．

やる夫

ナイキスト周波数を超えても，まだだいぶ振幅が残ってしまってるので，その部分が折り返されてしまうわけだお．

やらない夫

実際の周波数特性を見てみようか．さっきの $\omega_\textnormal{c}= 0.4\pi$ のバタワースフィルタのうち次数 $N = 2, 4, 8, 16$ のそれぞれをインパルス不変変換でディジタルフィルタ化したものの振幅特性だ．

$N = 2$ のときがわかりやすいかな． $\omega = \pi$ で折り返された結果，遮断域で振幅が下がりにくくなっている．

やる夫

なるほど…，ん? あれれ? 何かおかしくないかお? $N = 2$ のときは遮断域の振幅が大きくなっているけど，$N = 4$ のときはむしろやや小さくなってないかお? ていうか $N = 2$ のときも，通過域の振幅が小さくなってるお．折り返された分が足し合わされるのに，どうして振幅が減るのかお?

やらない夫

いや，別にそれはおかしくないぞ．振幅が足し合わされるんじゃなくて，振幅と位相からなる複素数の加算が起きるんだ．その結果振幅は増えることもあれば減ることもある．

やる夫

あー，そうか．足し合わせても常に振幅が大きくなるわけじゃないんだお．

やらない夫

まあいずれにせよ，低域通過フィルタはまだマシな方だ．それ以上に問題なのは，高域通過フィルタとか，帯域阻止フィルタとかの場合だな．何しろある周波数以降，無限大までずっと 1 に近い振幅を持つようなプロトタイプフィルタを使おうってことだ．エイリアシングが起きるというか，もはやエイリアシングの影響しか残らないくらいだ．

やる夫

ああ…，確かにそれは困るお．

やらない夫

というわけで，インパルス不変変換が使える状況は割と限られている．まあその限られた状況の範囲内であれば，プロトタイプフィルタの特性を少なくともサンプリング定理の意味で忠実に再現できる手法だということもできるかな．

16.3.2 双線形変換

やる夫

じゃあエイリアシングが無視できないようなときはどうするんだお?

やらない夫

インパルス不変変換とは別のやり方を考えよう．方針はこうだ．プロトタイプフィルタの $\Omega = -\infty$ 〜$\infty$ の範囲を，ディジタルフィルタの $\omega = -\pi$ 〜 $\pi$ の範囲に対応させるような変数変換を考える．もしそういうことができれば， $\Omega = -\infty$ 〜$\infty$ まで広がったプロトタイプフィルタの周波数特性を，ギュっと縮めてディジタルフィルタの $\omega = -\pi$ 〜 $\pi$ の範囲の周波数特性にすることができるだろう．これならエイリアシングの心配はない．

やる夫

んー，ピンと来ないお．

やらない夫

具体的な話を先にしてしまおうか．

$$s = 2 \frac{1 - z^{-1}}{1 + z^{-1}}$$ (16.26)

という式で $s$ から $z$ への変換を行うことを考えよう．

やる夫

なんだ，こりゃまた唐突だお．

やらない夫

この式はこの式でいろいろと意味づけができるんだが，その辺の話は置いといて，この変換でどういう周波数変換がされるかを見ていこう．

プロトタイプフィルタの伝達関数が $H(s)$ だとすると，その周波数特性は $s = j\Omega$ を代入して $H(j\Omega)$ と書ける．一方，式 (16.27) で変数変換して作ったディジタルフィルタの伝達関数は

$$H(2 \frac{1 - z^{-1}}{1 + z^{-1}})$$ (16.27)

になって，その周波数特性は $z = e^{j\omega}$ を代入して

$$H(2 \frac{1 - e^{-j\omega}}{1 + e^{-j\omega}})$$ (16.28)

と書ける．

やる夫

まあ淡々と代入してっただけだお．

やらない夫

これらが何を意味しているかというと

$$j\Omega = 2 \frac{1 - e^{-j\omega}}{1 + e^{-j\omega}}$$ (16.29)

が成り立つような周波数 $\Omega$ と $\omega$ が対応しているってことだな． $\Omega$ について整理すると

$$\Omega = \frac{2}{j} \frac{1 - e^{-j\omega}}{1 + e^{-j\omega}}$$ (16.30)

$$= \frac{2}{j} \frac{e^{-j\omega/2}}{e^{-j\omega/2}} \frac{e^{j\omega/2} - e^{-j\omega/2}}{e^{j\omega/2} + e^{-j\omega/2}}$$ (16.31)

$$= 2 \frac{(e^{j\omega/2} - e^{-j\omega/2})/2j}{(e^{j\omega/2} + e^{-j\omega/2})/2}$$ (16.32)

$$= 2 \frac{\sin\omega/2}{\cos\omega/2}$$ (16.33)

$$= 2 \tan\frac{\omega}{2}$$ (16.34)

になる．

やる夫

$\Omega = 2\tan\frac{\omega}{2}$ のところでプロトタイプフィルタとディジタルフィルタの周波数特性値が同じになるってことだお．

やらない夫

この関係をグラフにするとこんな感じだ． $\Omega = -\infty$ 〜$\infty$ が $\omega = -\pi$ 〜 $\pi$ に対応しているのがわかるだろう．

やる夫

あー，こういうことかお．しかし盛大にねじ曲がった関係だお．

やらない夫

無限の範囲を有限の範囲に対応付けようとするわけだから，対応関係に歪みが生じるのは無理もない．その点を認めてしまえるなら，エイリアシング無しでディジタルフィルタを作れるとても便利なものだと考えることができる．この変数変換を，双線形変換とか，双1次変換とか呼んでいる．

やる夫

この場合，遮断周波数とかはどうやって決めればいいのかお?

やらない夫

ちゃんと逆算してやらないとダメだな．遮断周波数が $\omega_\textnormal{c}$ のディジタルフィルタを作りたかったら， $\Omega_\textnormal{c}= 2\tan\frac{\omega_\textnormal{c}}{2}$ が遮断周波数になるようなプロトタイプフィルタを用意してやる必要がある．

やる夫

なるほど，で，あとはその伝達関数を改めて双線形変換してやればいいんだお．

やらない夫

そういうことだな．ただし，周波数が歪んで対応付けられる分，例えば遮断特性とかは正確に保存されないので注意が必要だ．その辺り，インパルス不変変換のようなわけにはいかない．

やる夫

まあそりゃそうだお．

やらない夫

実際の計算例を見てみようか．プロットされてる内容はインパルス不変変換のときと全く同じだ．エイリアシングの影響が出てないことがわかると思う．

やる夫

その代わり遮断特性はプロトタイプフィルタとはちょっと違った感じになってるお．

やらない夫

それが $\Omega$ と $\omega$ の間のねじ曲がった対応関係の影響だ．

まあこんなところで，ディジタルフィルタ設計の基本の基本の話はおしまいだ．まとめよう．

• 所望の特性を持つフィルタのインパルス応答を適当な窓関数で打ち切ることで FIR フィルタを設計することができる．FIR フィルタの窓関数設計法と呼ぶ．
  • 周波数特性を最小2乗近似する意味で最適なフィルタは矩形窓により設計できるが，リップルが大きい．
  • 両端をなめらかに打ち切る窓関数を使えば，遮断特性は犠牲になるもののリップルを抑えることができる．
• 適当なアナログフィルタをプロトタイプとして IIR ディジタルフィルタに変換することができる．IIR フィルタの間接設計法と呼ぶ．
  • プロトタイプフィルタのインパルス応答をサンプリングしたものをディジタルフィルタのインパルス応答とするような変換方法を，インパルス不変変換と呼ぶ．プロトタイプフィルタの周波数特性がナイキスト周波数以上の成分を多く含んでいる場合，エイリアシングによって大きな誤差が生じる．
  • エイリアシングを避けることのできる変換方法として，双線形変換がよく用いられる．変換の際に周波数特性は $\Omega = 2\tan\frac{\omega}{2}$ に従って歪んでしまうので，その点をあらかじめ考慮して仕様を与える必要がある．
• フィルタの次数を大きくするほど理想的な特性に近づけやすくなるが，コストや遅延時間が発生する．
• 同じ次数であれば FIR より IIR フィルタの方が遮断特性を急峻にしやすい．しかし IIR フィルタは線形位相特性を持てない
• 要求仕様を満たすフィルタを設計する際，FIR フィルタの場合は試行錯誤が必要となる．IIR フィルタの場合は試行錯誤の要らない設計法もある．

やる夫

基本の基本ってことは，まだ他にもいろんな設計法があるってことかお．

やらない夫

もちろんだ．俺たちの戦いはまだ始まったばかりだ．

(やらない夫先生の次回作にご期待ください)