FIR フィルタの設計 — 最小2乗近似と窓関数

信号処理論入門　第11章 §11.1　★ 理想LPFの打切り（最小2乗）と窓（Kaiser）

仕上げの章である。フィルタとは、入力信号に含まれる必要な成分を通過させ、不要な成分を阻止するシステムのことで、本書では離散時間 LTI システムで記述できる周波数選択フィルタを扱う。フィルタの実体は $h[n]$ であり $H(z)$ であり$H(\mathrm{e}^{\mathrm{j}\Omega})$ であり、差分方程式の係数 $\{a_p\},\{b_q\}$（フィルタ係数）である。フィルタを設計するとは、所望の周波数特性を持つようにフィルタ係数を定めることにほかならない。FIR フィルタでは係数は $\{b_q\}$ のみ、IIR フィルタでは $\{a_p\},\{b_q\}$ の両方を定める。

理想の周波数選択特性は4種類ある：LPF（ローパス：$\left\vert \Omega\right\vert \le\Omega_C$ を通す）、HPF（ハイパス）、BPF（バンドパス：$\Omega_1\le\left\vert \Omega\right\vert \le\Omega_2$）、BSF（バンドストップ）。いずれも通過域で振幅比 1（複素数として）、阻止域で 0 の矩形が理想である。境界の周波数 $\Omega_C$ を遮断周波数（カットオフ周波数）と呼ぶ。

手順1：理想応答の逆変換。理想の周波数応答 $H_{\mathrm{ideal}}(\mathrm{e}^{\mathrm{j}\Omega})$を実現するインパルス応答は、逆 DTFT で逆算できる。LPF なら

$$h_{\mathrm{ideal}}[n] =\frac{1}{2\pi}\int_{-\Omega_C}^{\Omega_C}\mathrm{e}^{\mathrm{j}\Omega n}\mathop{}\!\mathrm{d}\Omega =\frac{1}{\pi}\int_0^{\Omega_C}\cos(\Omega n)\mathop{}\!\mathrm{d}\Omega =\frac{\sin(\Omega_C n)}{n\pi}$$

（$n=0$ では $\Omega_C/\pi$）。これは第3章の「矩形 $\leftrightarrow$ $\operatorname{sinc}$」の離散時間版で、無限に長く、しかも両側に伸びる（非因果的）。そのままでは実現できない。

手順2：打ち切り＝最小2乗近似。そこで長さを $n\in[-M/2,M/2]$（$M$：次数、タップ数）に制限して近似する。実際の応答と理想応答の平均2乗誤差

$$J=\frac{1}{2\pi}\int_{-\pi}^{\pi} \left\vert H(\mathrm{e}^{\mathrm{j}\Omega})-H_{\mathrm{ideal}}(\mathrm{e}^{\mathrm{j}\Omega})\right\vert ^2\mathop{}\!\mathrm{d}\Omega =\sum_{\left\vert n\right\vert \le M/2}\left\vert h[n]-h_{\mathrm{ideal}}[n]\right\vert ^2 +\sum_{\left\vert n\right\vert >M/2}\left\vert h_{\mathrm{ideal}}[n]\right\vert ^2$$

（Parseval の等式で周波数領域から時間領域に書き換えた）を最小化する $h[n]$ は、第2項が $h$ に依らない定数であることから、

$$h[n]=h_{\mathrm{ideal}}[n]\qquad(n\in[-M/2,M/2])$$

—つまり単純に切り出したものが最小2乗最適である。これは第2章で「Fourier 係数での打ち切りが最良近似」を示した議論の再演であり、実際 DTFT は「周期 $2\pi$ の関数 $H(\mathrm{e}^{\mathrm{j}\Omega})$ の Fourier 級数」だから同じ定理がそのまま適用できる。

手順3：因果化。切り出した $h[n]$ は偶関数なので $H(\mathrm{e}^{\mathrm{j}\Omega})$ は実数、すなわち位相 0 の零位相フィルタである（当然、非因果的）。$M/2$ サンプル右へずらせば

$$h_{\mathrm{causal}}[n]=h[n-M/2],\qquad H_{\mathrm{causal}}(\mathrm{e}^{\mathrm{j}\Omega})=\mathrm{e}^{-\mathrm{j} M\Omega/2}H(\mathrm{e}^{\mathrm{j}\Omega})$$

となり、因果的かつ線形位相（係数偶対称＝相反多項式、第9章）のFIR フィルタが完成する。

手順4：窓掛け。ただし打ち切りは急峻な切断なので、遮断周波数のまわりにリップル（Gibbs 現象、第2章・第6章）が現れる。最小2乗の意味で最適でも、実用上はリップルが小さいことの方が重要な場合が多い。そこで切り出しに滑らかな窓関数を掛けてリップルを抑える。実用的な選択がKaiser 窓

$$w[n]=\frac{I_0\Bigl(\beta\sqrt{1-\bigl(\frac{n-N/2}{N/2}\bigr)^2}\Bigr)}{I_0(\beta)} \quad(n\in[0,N]),\qquad I_0(x)=\sum_{m=0}^{\infty}\Bigl(\frac{x^m}{m!\,2^m}\Bigr)^2$$

である（$I_0$：0次第1種変形 Bessel 関数）。パラメータ $\beta$ で形が連続的に変わり、$\beta=0$ で方形窓、$\beta\approx5$ でほぼ Hamming 窓、$\beta\approx6$ でほぼ Hann 窓、$\beta\approx8.6$ でほぼ Blackman 窓を再現する—第6章の窓たちの「連続化」と思えばよい。$\beta$ を大きくすると阻止域のリップルは小さくなるが遷移域は広がる、というおなじみのトレードオフを1本のつまみで調整できる。

設計仕様は次の用語で与えられる：通過域の振幅を $1\pm\delta$、阻止域の振幅を$\delta$ 以下に収め、通過域と阻止域の間の遷移域の幅を $\Delta\Omega$ とする。減衰量を $A=-20\log_{10}\delta$（dB）と書くと、Kaiser 窓のパラメータは経験式

$$\beta=\begin{cases} 0.1102(A-8.7), & A\ge50\\ 0.5842(A-21)^{0.4}+0.07886(A-21), & 21<A<50\\ 0, & A\le21 \end{cases} \qquad N\approx\frac{A-7.95}{2.285\,\Delta\Omega}$$

で決める。なお、最小2乗（＋窓）ではなく誤差の最大値を最小化（minimax 近似）するのが Parks–McClellan 法で、リップルが所望特性の上下に等幅で並ぶ等リップル特性のフィルタが得られる（設計には数値最適化を用いる）。