ディジタルフィルタ設計 — FIR/IIRフィルタの理論と実装

音楽を聴くとき、イコライザーで低音を強調したり高音を抑えたりした経験はないでしょうか。あるいは心電図のノイズを除去して心拍のパターンだけを取り出す場面を想像してみてください。これらはすべて、信号の中から特定の周波数成分だけを取り出したり除去したりする「フィルタリング」という操作です。アナログの世界ではコンデンサやコイルでフィルタを構成しますが、コンピュータの中では数値演算だけでフィルタを実現できます。これが ディジタルフィルタ です。

ディジタルフィルタは信号処理のあらゆる分野で中核的な役割を果たしており、その理解は多くの応用に直結します。

通信: 受信信号から目的の帯域だけを抽出するバンドパスフィルタ、隣接チャネルの干渉を除去するノッチフィルタ
音声処理: ノイズ除去、エコーキャンセル、音声認識の前処理
計測・制御: センサデータの平滑化、振動解析、フィードバック制御系のコントローラ設計
画像処理: エッジ検出（微分フィルタ）、ぼかし（平滑化フィルタ）、シャープニング
生体信号処理: 心電図・脳波からのアーティファクト除去、特定の周波数帯の抽出

本記事の内容

ディジタルフィルタの基礎（差分方程式、伝達関数 $H(z)$）
FIR フィルタ: 直線位相特性と窓関数法（ハミング窓、ブラックマン窓等）による設計
IIR フィルタ: バターワース・チェビシェフのアナログプロトタイプと双一次変換
FIR vs IIR の比較（安定性、計算量、位相特性）
フィルタの周波数応答、インパルス応答、零点-極点配置
Python 実装: FIR/IIR フィルタ設計 + ノイズ除去デモ

前提知識

この記事を読む前に、以下の知識があると理解が深まります。

フーリエ変換の基礎
z変換の基礎知識（$z = e^{j\omega T}$ の関係、z平面における極と零点）

ディジタルフィルタの基礎

差分方程式: フィルタの「レシピ」

ディジタルフィルタとは、入力信号の数列 $x[n]$ から出力信号の数列 $y[n]$ を生成する演算システムです。アナログフィルタが微分方程式で記述されるのに対し、ディジタルフィルタは 差分方程式 で記述されます。

最も一般的な線形時不変（LTI）ディジタルフィルタの差分方程式は

$$ y[n] = \sum_{k=0}^{M} b_k \, x[n-k] – \sum_{k=1}^{N} a_k \, y[n-k] $$

と書けます。ここで $b_k$（$k = 0, 1, \dots, M$）は フィードフォワード係数、$a_k$（$k = 1, 2, \dots, N$）は フィードバック係数 です。$M$ はフィードフォワード部の次数、$N$ はフィードバック部の次数です。

この式を料理のレシピにたとえると、$b_k$ が「材料（入力信号）の配合比」、$a_k$ が「前回の出来具合（過去の出力）を参考にする度合い」です。$a_k = 0$（$\forall k \geq 1$）なら過去の出力を参照しない「前向き」なレシピ、$a_k \neq 0$ なら過去の出力も取り入れる「フィードバック」付きのレシピです。

伝達関数 $H(z)$

差分方程式にz変換を適用すると、時間領域の畳み込みが周波数領域の積として表現できます。z変換の遅延特性 $\mathcal{Z}\{x[n-k]\} = z^{-k} X(z)$ を使うと

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

が得られます。$Y(z)$ について整理すると

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

となり、伝達関数 $H(z) = Y(z) / X(z)$ は

$$ H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{N} a_k z^{-k}} = \frac{b_0 + b_1 z^{-1} + \cdots + b_M z^{-M}}{1 + a_1 z^{-1} + \cdots + a_N z^{-N}} $$

と表されます。分子を零にする $z$ の値が 零点（zero）、分母を零にする $z$ の値が 極（pole） です。分子・分母をそれぞれ因数分解すると

$$ H(z) = b_0 \frac{\prod_{k=1}^{M} (1 – z_k z^{-1})}{\prod_{k=1}^{N} (1 – p_k z^{-1})} $$

と書けます。ここで $z_k$ は零点、$p_k$ は極です。

周波数応答

ディジタルフィルタの周波数応答は、伝達関数 $H(z)$ に $z = e^{j\omega}$ を代入して得られます（ここで $\omega$ は正規化角周波数、$0 \leq \omega \leq \pi$）。

$$ H(e^{j\omega}) = |H(e^{j\omega})| \, e^{j\angle H(e^{j\omega})} $$

$|H(e^{j\omega})|$ が 振幅応答（ゲイン特性）、$\angle H(e^{j\omega})$ が 位相応答 です。振幅応答はどの周波数をどれだけ通すか/減衰させるかを、位相応答はどの周波数をどれだけ遅延させるかを表します。

理想的なローパスフィルタは、カットオフ周波数 $\omega_c$ 以下の周波数成分を無変化で通過させ、$\omega_c$ より上の成分を完全に遮断します。

$$ H_{\text{ideal}}(e^{j\omega}) = \begin{cases} 1 & (|\omega| \leq \omega_c) \\ 0 & (\omega_c < |\omega| \leq \pi) \end{cases} $$

しかし、この理想フィルタのインパルス応答 $h[n]$ は

$$ h[n] = \frac{\sin(\omega_c n)}{\pi n} $$

で、これは $n \to \pm\infty$ で初めてゼロになる無限長の系列（sinc関数）です。実際のシステムでは無限長のフィルタは実現できないため、これをどのように有限長で近似するかがフィルタ設計の核心です。

ディジタルフィルタの基礎的な枠組みが整いました。ここから、フィルタの2大分類である FIR フィルタと IIR フィルタについて、それぞれの理論と設計法を詳しく見ていきましょう。まず FIR フィルタから始めます。

FIR フィルタ

FIR フィルタとは

FIR（Finite Impulse Response：有限インパルス応答）フィルタ は、フィードバック係数がすべてゼロ（$a_k = 0$, $\forall k \geq 1$）のフィルタです。伝達関数は分子多項式のみで構成されます。

$$ H(z) = \sum_{k=0}^{M} b_k z^{-k} = b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots + b_M z^{-M} $$

フィードバックがないため、入力信号の有限個のサンプル $x[n], x[n-1], \dots, x[n-M]$ だけから出力が決まります。入力にインパルス $\delta[n]$ を与えると、出力は $h[n] = b_n$（$n = 0, 1, \dots, M$）となり、$n > M$ でゼロになります。つまり、インパルス応答が文字通り「有限」の長さ $M+1$ で終了します。

FIR フィルタを日常のたとえで言うなら、「直近の $M+1$ 日間の気温を加重平均して、今日の”平均気温”を出す」操作です。過去の計算結果は使わず、元のデータ（入力）だけから毎回計算します。

FIR フィルタの3つの優れた性質

FIR フィルタには、IIR フィルタにはない重要な利点があります。

安定性の保証: FIR フィルタの伝達関数 $H(z)$ の極はすべて $z = 0$（$N$ 重極）にあり、z平面の単位円内に位置します。したがって FIR フィルタは 常に安定 です。フィードバックがないため発散する心配がありません。

直線位相特性: フィルタ係数が対称（$b_k = b_{M-k}$）または反対称（$b_k = -b_{M-k}$）の場合、FIR フィルタは 厳密な直線位相応答 を持ちます。すなわち

$$ \angle H(e^{j\omega}) = -\frac{M}{2} \omega + \phi_0 $$

の形になります。位相が周波数に対して直線的に変化するということは、すべての周波数成分が同じ時間遅延 $M/2$ サンプルを受けることを意味します。これは 群遅延が一定 であることと等価で、波形の歪みが生じません。

音声信号や画像など、波形の形状が重要な応用では、直線位相は極めて重要です。IIR フィルタでは厳密な直線位相は実現できないため、これはFIRフィルタの決定的な強みです。

実装の単純さ: FIR フィルタの演算は入力の加重和（畳み込み）のみであり、フィードバックループが不要です。ハードウェア実装（FPGA, DSP）でもパイプライン化が容易で、並列処理に適しています。

窓関数法による FIR フィルタ設計

窓関数法は、FIR フィルタ設計の最も直感的で広く使われる手法です。基本的なアイデアは、理想フィルタの無限長インパルス応答 $h_d[n]$ を有限の 窓関数 $w[n]$ で切り取ることです。

ステップ1 — 理想インパルス応答の計算: 目標とする理想周波数応答 $H_d(e^{j\omega})$ の逆離散時間フーリエ変換（IDTFT）を計算します。理想ローパスフィルタ（カットオフ角周波数 $\omega_c$）の場合

$$ h_d[n] = \frac{\omega_c}{\pi} \text{sinc}\left(\frac{\omega_c n}{\pi}\right) = \frac{\sin(\omega_c n)}{\pi n} $$

です。$n = 0$ のとき $h_d[0] = \omega_c / \pi$ と定義します。

ステップ2 — 窓関数による切り出し: 無限長の $h_d[n]$ を窓関数 $w[n]$（長さ $M+1$）で打ち切ります。

$$ h[n] = h_d[n – M/2] \cdot w[n], \quad n = 0, 1, \dots, M $$

ここで $M/2$ だけシフトしているのは、因果的なフィルタ（$h[n] = 0$ for $n < 0$）にするためです。

ステップ3 — 窓関数の選択: 窓関数の形状が、設計されたフィルタの性能（遷移帯域幅とサイドローブ抑圧のトレードオフ）を決めます。

代表的な窓関数

窓関数法の性能は窓関数の選択で大きく変わります。以下に代表的な窓関数をまとめます。

矩形窓（Rectangular）:

$$ w[n] = 1, \quad 0 \leq n \leq M $$

最も単純な窓で、遷移帯域幅は最も狭いですが、サイドローブが $-13$ dB と大きく、リップルが目立ちます。

ハミング窓（Hamming）:

$$ w[n] = 0.54 – 0.46 \cos\left(\frac{2\pi n}{M}\right), \quad 0 \leq n \leq M $$

サイドローブが $-43$ dB に抑圧され、実用上最も広く使われる窓関数の1つです。遷移帯域幅は矩形窓の約2倍になりますが、通過帯域のリップルが大幅に改善されます。

ハニング窓（Hanning / von Hann）:

$$ w[n] = 0.5 – 0.5 \cos\left(\frac{2\pi n}{M}\right), \quad 0 \leq n \leq M $$

ハミング窓とよく似ていますが、両端で厳密にゼロになるため、スペクトル解析でよく使われます。サイドローブは $-31$ dB です。

ブラックマン窓（Blackman）:

$$ w[n] = 0.42 – 0.5 \cos\left(\frac{2\pi n}{M}\right) + 0.08 \cos\left(\frac{4\pi n}{M}\right), \quad 0 \leq n \leq M $$

3項のコサイン和で構成され、サイドローブが $-58$ dB と非常に低いですが、遷移帯域幅はハミング窓より広くなります。高いサイドローブ抑圧が必要な応用（隣接チャネル干渉の除去など）に適しています。

カイザー窓（Kaiser）:

$$ w[n] = \frac{I_0\left(\beta \sqrt{1 – \left(\frac{2n/M – 1}{1}\right)^2}\right)}{I_0(\beta)}, \quad 0 \leq n \leq M $$

ここで $I_0$ は第1種0次の修正ベッセル関数、$\beta$ はパラメータです。$\beta$ を調整することで、遷移帯域幅とサイドローブ抑圧のトレードオフを連続的に制御できます。$\beta = 0$ で矩形窓、$\beta$ を大きくするほどサイドローブが低くなる代わりに遷移帯域が広がります。

各窓関数のトレードオフをまとめます。

窓関数	サイドローブ [dB]	遷移帯域幅（概算）	最小阻止域減衰 [dB]
矩形	$-13$	$4\pi / (M+1)$	$-21$
ハニング	$-31$	$8\pi / (M+1)$	$-44$
ハミング	$-43$	$8\pi / (M+1)$	$-53$
ブラックマン	$-58$	$12\pi / (M+1)$	$-74$

窓関数法の設計手順のまとめ

窓関数法でフィルタを設計する手順を整理します。

仕様の決定: カットオフ周波数 $\omega_c$、阻止域減衰 $A_s$ [dB]、遷移帯域幅 $\Delta\omega$ を決める
窓関数の選択: $A_s$ の要件を満たす窓関数を上の表から選ぶ
フィルタ次数の決定: $M \approx (A_s – 8) / (2.285 \Delta\omega)$（カイザーの近似公式）で次数を見積もる
理想インパルス応答の計算: $h_d[n-M/2]$ を計算
窓関数の適用: $h[n] = h_d[n-M/2] \cdot w[n]$ を計算
周波数応答の検証: 設計されたフィルタが仕様を満たすか確認

FIR フィルタの理論と設計法がわかりました。FIR は安定性と直線位相が保証される優れた方式ですが、急峻な遮断特性を実現するにはフィルタ次数 $M$ を大きくする必要があり、計算量が増えます。少ない次数で急峻な特性を実現したい場合は、フィードバックを利用する IIR フィルタが有力な選択肢になります。

IIR フィルタ

IIR フィルタとは

IIR（Infinite Impulse Response：無限インパルス応答）フィルタ は、フィードバック係数 $a_k \neq 0$（少なくとも1つ）を持つフィルタです。フィードバックにより、入力にインパルスを与えると出力は理論的には無限に続きます（ただし安定なフィルタでは指数関数的に減衰します）。

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

フィードバックの存在は、音響空間におけるエコーや反響に似ています。部屋の中で手を叩くと、壁に反射した音が繰り返し聞こえます。IIR フィルタの内部でも、出力が再び入力に「反射」されて繰り返し処理されるため、少ない係数でも効果的なフィルタリングが実現できます。

IIR フィルタの最大の利点は 効率性 です。同じ遮断特性（阻止域減衰、遷移帯域幅）を実現するのに、FIR フィルタの $1/5$ から $1/10$ の次数で済むことが多いです。これはフィードバックの「増幅効果」により、少数の極で急峻な特性を作れるためです。

アナログプロトタイプフィルタ

IIR フィルタの設計は、「よく研究されたアナログフィルタをディジタルに変換する」というアプローチが主流です。まず連続時間のアナログフィルタ（プロトタイプ）を設計し、それを離散時間に変換してディジタルフィルタを得ます。

代表的なアナログプロトタイプフィルタを2つ紹介します。

バターワースフィルタ

バターワース（Butterworth）フィルタは、通過帯域で最大限に平坦な振幅応答 を持つフィルタです。$N$ 次バターワースフィルタの振幅応答の二乗は

$$ |H_a(j\Omega)|^2 = \frac{1}{1 + (\Omega / \Omega_c)^{2N}} $$

で与えられます。ここで $\Omega$ はアナログ角周波数、$\Omega_c$ はカットオフ角周波数（$-3$ dB 点）です。

この式のポイントは、$\Omega = 0$ の近傍で $|H_a|^2$ の $\Omega$ に関する導関数が $2N-1$ 次まで全てゼロになることです。これが「最大平坦（maximally flat）」と呼ばれる所以で、通過帯域にリップルが全くありません。

$\Omega = \Omega_c$ で $|H_a|^2 = 1/2$（$-3$ dB）、$\Omega \gg \Omega_c$ では $|H_a|^2 \approx (\Omega_c / \Omega)^{2N}$ なので、阻止域では $20N$ dB/decade（$6N$ dB/octave）の傾斜で減衰します。次数 $N$ を上げるほど、理想フィルタの急峻な特性に近づきます。

バターワースフィルタの極は、s平面上の左半平面に等間隔で配置されます。$N$ 次フィルタの極は

$$ s_k = \Omega_c \, e^{j\pi(2k+N-1)/(2N)}, \quad k = 0, 1, \dots, N-1 $$

のうち、$\text{Re}(s_k) < 0$ となる $N$ 個です。これらの極は半径 $\Omega_c$ の円上に等間隔で並びます。

チェビシェフI型フィルタ

チェビシェフ（Chebyshev）I型フィルタは、通過帯域に 等リップル を許容する代わりに、バターワースフィルタよりも急峻な遮断特性を実現します。$N$ 次チェビシェフI型フィルタの振幅応答の二乗は

$$ |H_a(j\Omega)|^2 = \frac{1}{1 + \epsilon^2 T_N^2(\Omega / \Omega_c)} $$

です。ここで $\epsilon$ は通過帯域リップルの大きさを制御するパラメータ、$T_N(x)$ は $N$ 次のチェビシェフ多項式です。

チェビシェフ多項式 $T_N(x)$ は再帰的に定義されます。

$$ T_0(x) = 1, \quad T_1(x) = x, \quad T_{N+1}(x) = 2x T_N(x) – T_{N-1}(x) $$

$|x| \leq 1$（通過帯域）では $|T_N(x)| \leq 1$ なので振幅応答は $1/(1+\epsilon^2)$ 以上に保たれますが、$T_N$ が $\pm 1$ の間を振動するため等リップルが生じます。$|x| > 1$（阻止域）では $T_N$ が急速に増大し、急峻な減衰を実現します。

通過帯域リップル $R_p$ [dB] と $\epsilon$ の関係は

$$ \epsilon = \sqrt{10^{R_p/10} – 1} $$

です。$R_p = 1$ dB ならば $\epsilon \approx 0.509$ です。

同じ次数 $N$ で比較すると、チェビシェフI型はバターワースよりも阻止域での減衰が大きく、遷移帯域が狭くなります。この急峻さの代償として通過帯域にリップルが生じます。

チェビシェフII型フィルタ

チェビシェフII型（逆チェビシェフ）フィルタは、通過帯域が最大平坦（バターワースと同様）で、阻止域に等リップルを持つフィルタです。通過帯域のリップルが許容できない場合に使われます。

これらのアナログプロトタイプフィルタをディジタルフィルタに変換する必要があります。その変換法として最も広く使われているのが双一次変換です。

双一次変換（Bilinear Transform）

双一次変換は、アナログフィルタの伝達関数 $H_a(s)$ をディジタルフィルタの伝達関数 $H(z)$ に変換する手法です。s平面からz平面への写像を

$$ s = \frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} $$

と定義します。ここで $T$ はサンプリング周期です。$H_a(s)$ のこの式を代入して $H(z)$ を得ます。

$$ H(z) = H_a\left(\frac{2}{T} \cdot \frac{1 – z^{-1}}{1 + z^{-1}}\right) $$

双一次変換がなぜうまく機能するかを理解するために、この写像の性質を見てみましょう。

安定性の保存: s平面の左半平面（$\text{Re}(s) < 0$）がz平面の単位円内部（$|z| < 1$）に写されます。$j\Omega$ 軸（虚軸）は単位円上に写されます。したがって、安定なアナログフィルタは安定なディジタルフィルタに変換されます。

周波数の歪み（ワーピング）: s平面の虚軸 $s = j\Omega$ とz平面の単位円 $z = e^{j\omega}$ の関係は

$$ j\Omega = \frac{2}{T} \cdot \frac{1 – e^{-j\omega}}{1 + e^{-j\omega}} = j\frac{2}{T} \tan\left(\frac{\omega}{2}\right) $$

となり、アナログ周波数 $\Omega$ とディジタル周波数 $\omega$ の間に非線形な関係

$$ \Omega = \frac{2}{T} \tan\left(\frac{\omega}{2}\right) $$

が成り立ちます。アナログ周波数 $\Omega \in [0, \infty)$ がディジタル周波数 $\omega \in [0, \pi)$ に圧縮されます。この非線形写像を 周波数ワーピング と呼びます。

周波数ワーピングの影響は、ディジタルフィルタのカットオフ周波数がアナログの仕様からずれることです。これを補正するために、設計前にアナログのカットオフ周波数を プリワーピング します。

$$ \Omega_c = \frac{2}{T} \tan\left(\frac{\omega_c}{2}\right) $$

ここで $\omega_c$ は目標とするディジタルカットオフ周波数です。プリワーピングした $\Omega_c$ でアナログフィルタを設計し、双一次変換すると、正確に $\omega_c$ にカットオフが配置されたディジタルフィルタが得られます。

IIR フィルタの設計手順

IIR フィルタ設計の全体的な流れをまとめます。

ディジタル仕様の決定: カットオフ周波数 $\omega_c$、通過帯域リップル $R_p$、阻止域減衰 $A_s$、遷移帯域幅 $\Delta\omega$
プリワーピング: $\Omega_c = (2/T)\tan(\omega_c/2)$ でアナログカットオフ周波数を計算
プロトタイプ選択と次数決定: バターワースなら $N \geq \log(10^{A_s/10}-1)/(2\log(\Omega_s/\Omega_c))$、チェビシェフなら $N \geq \cosh^{-1}(\sqrt{(10^{A_s/10}-1)/\epsilon^2})/\cosh^{-1}(\Omega_s/\Omega_c)$ で最小次数を計算
アナログフィルタの設計: $H_a(s)$ の極（と零点）を配置
双一次変換: $s = (2/T)(1-z^{-1})/(1+z^{-1})$ を代入して $H(z)$ を得る
検証: 周波数応答が仕様を満たすか確認

IIR フィルタの注意点: 安定性

IIR フィルタにはフィードバックがあるため、不安定 になる可能性があります。安定性の条件は、$H(z)$ の全ての極が z 平面の単位円の内部（$|p_k| < 1$, $\forall k$）にあることです。

双一次変換で設計した IIR フィルタは自動的に安定ですが、量子化（固定小数点演算）により係数が丸められると、極が単位円の外に移動して不安定になることがあります。特に高次のフィルタで顕著です。この問題を回避するため、高次フィルタは 2次セクション（SOS: Second-Order Sections） の縦続接続として実装するのが一般的です。

FIR フィルタと IIR フィルタの理論がそろいました。ここで、両者を様々な観点から比較し、それぞれの長所と短所を整理しましょう。

FIR vs IIR の比較

2つのフィルタ方式を、設計者が判断に迫られる主要な観点で比較します。

安定性

FIR フィルタは 常に安定 です。極がすべて原点にあるため、係数の量子化による不安定化も起こりません。一方、IIR フィルタはフィードバックの存在により安定性の検証が不可欠です。高次フィルタを直接形（Direct Form）で実装すると、係数量子化で極が単位円外に出るリスクがあります。

位相特性

FIR フィルタは 厳密な直線位相（群遅延一定）を実現できます。音声、画像、通信の波形整形など、位相歪みが許されない応用では FIR が必須です。

IIR フィルタの位相応答は非線形であり、群遅延が周波数によって変化します。特にカットオフ周波数付近で群遅延が大きくなります。位相の線形化が必要な場合は、全域通過フィルタ を縦続接続して群遅延を等化する方法がありますが、設計が複雑になります。

計算量

同じ遮断特性（阻止域減衰、遷移帯域幅）を実現するのに必要な乗算回数は、IIR の方が大幅に少ないのが一般的です。

たとえば、$-60$ dB の阻止域減衰を持つローパスフィルタを設計する場合、バターワース IIR では 4 – 6 次程度で十分ですが、FIR（ハミング窓）では 50 – 100 次程度が必要です。1サンプルあたりの乗算回数は、IIR が $2N$（$N$ は次数）、FIR が $M+1$（$M$ はフィルタ長）なので、IIR は $10 – 20$ 回、FIR は $50 – 100$ 回となります。

ただし FIR フィルタは畳み込みのみで構成されるため、FFT を使った高速畳み込み（overlap-save 法など）で大幅に高速化できます。

設計の容易さ

FIR は窓関数法やパークス・マクレラン法（等リップル最適設計）など、直感的で体系的な設計法が確立されています。仕様から一意にフィルタが決まり、試行錯誤が少ないのが利点です。

IIR はアナログプロトタイプからの変換が必要で、設計の自由度が高い反面、安定性の確認や係数量子化への配慮が必要です。

用途別の選択指針

用途	推奨	理由
音声信号処理	FIR	直線位相で波形歪みなし
リアルタイム制御	IIR	低次数で低遅延
通信の波形整形	FIR	直線位相 + ISI最小化
スペクトル解析の前処理	IIR	計算効率 + 位相は重要でない
生体信号のオフライン解析	FIR	前方向+後方向フィルタで零位相を実現可能
組込みシステム	IIR	メモリ・計算量の制約

このように、FIR と IIR はそれぞれ明確な得意分野を持っています。実際の設計では、要件に応じてどちらを使うか（あるいは組み合わせるか）を判断します。

比較の理論を踏まえた上で、零点-極点配置がフィルタの周波数応答をどのように決定するかを詳しく見ていきましょう。

零点-極点配置と周波数応答

零点と極の直感的な意味

z平面上の零点と極は、フィルタの周波数応答を幾何学的に理解するための強力なツールです。

$z = e^{j\omega}$ を単位円上の点とすると、周波数 $\omega$ における周波数応答の振幅は

$$ |H(e^{j\omega})| = |b_0| \frac{\prod_{k=1}^{M} |e^{j\omega} – z_k|}{\prod_{k=1}^{N} |e^{j\omega} – p_k|} $$

で与えられます。$|e^{j\omega} – z_k|$ は単位円上の点から零点 $z_k$ までの距離、$|e^{j\omega} – p_k|$ は極 $p_k$ までの距離です。

この式から、直感的な理解が得られます。

零点が単位円に近いとき: 単位円上の点が零点に近づくと、分子のある項 $|e^{j\omega} – z_k|$ が小さくなり、$|H(e^{j\omega})|$ が減少します。零点が単位円上にあれば、その角度の周波数成分は完全に遮断されます（ゲイン = 0）。

極が単位円に近いとき: 単位円上の点が極に近づくと、分母のある項 $|e^{j\omega} – p_k|$ が小さくなり、$|H(e^{j\omega})|$ が増大します。極が単位円に近いほど鋭い共振ピークが現れます。極が単位円上に到達すると振幅が無限大（不安定）になります。

零点-極点配置の設計例

具体例として、単純なフィルタの零点-極点配置と周波数応答の関係を見てみましょう。

ノッチフィルタ（特定周波数の除去）: 角周波数 $\omega_0$ を除去するには、$z = e^{j\omega_0}$ と $z = e^{-j\omega_0}$ に零点を配置します。これにより $\omega = \omega_0$ で振幅がゼロになります。ただし零点だけでは帯域全体の特性が変わるため、零点の近くに極を配置（$p = r e^{\pm j\omega_0}$, $r < 1$）して、$\omega_0$ の近傍以外のゲインを1に近づけます。$r$ を1に近づけるほどノッチの幅が狭くなります。

共振器（特定周波数の強調）: $\omega_0$ を強調するには、$z = r e^{\pm j\omega_0}$（$r$ は1に近い正の実数）に極を配置します。$r$ を1に近づけるほど鋭い共振が得られます。

このような幾何学的な理解は、フィルタの動作を直感的に把握する上で非常に有用です。特にフィルタの特性を微調整する際に、零点と極の位置を視覚的に操作できます。

零点-極点配置の理論を把握したところで、いよいよ Python でフィルタを設計し、実際の信号にフィルタリングを適用してみましょう。

Python での FIR フィルタ設計

窓関数法によるローパスフィルタ

まず、窓関数法で FIR ローパスフィルタを設計し、各窓関数の効果を比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# フィルタ仕様
fs = 1000          # サンプリング周波数 [Hz]
fc = 100           # カットオフ周波数 [Hz]
M = 51             # フィルタ長（タップ数）
wc = 2 * np.pi * fc / fs  # 正規化カットオフ角周波数

# 理想ローパスフィルタのインパルス応答
n = np.arange(M)
n_center = (M - 1) / 2
h_ideal = np.where(n == n_center,
                   wc / np.pi,
                   np.sin(wc * (n - n_center)) / (np.pi * (n - n_center)))

# 各窓関数
windows = {
    'Rectangular': np.ones(M),
    'Hamming': np.hamming(M),
    'Hanning': np.hanning(M),
    'Blackman': np.blackman(M),
}

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 窓関数の波形
ax = axes[0, 0]
for name, w in windows.items():
    ax.plot(n, w, linewidth=1.5, label=name)
ax.set_xlabel('Sample index n')
ax.set_ylabel('Window value')
ax.set_title('Window Functions')
ax.legend()
ax.grid(True, alpha=0.3)

# FIRフィルタのインパルス応答
ax = axes[0, 1]
for name, w in windows.items():
    h = h_ideal * w
    ax.plot(n, h, linewidth=1.2, label=name, alpha=0.8)
ax.set_xlabel('Sample index n')
ax.set_ylabel('Amplitude')
ax.set_title('FIR Filter Impulse Response')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# 周波数応答（振幅）
ax = axes[1, 0]
for name, w in windows.items():
    h = h_ideal * w
    freq, H = signal.freqz(h, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H) + 1e-15), linewidth=1.5, label=name)
ax.axvline(x=fc, color='gray', linestyle='--', alpha=0.5, label=f'fc={fc}Hz')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('Frequency Response (Magnitude)')
ax.set_xlim(0, fs/2)
ax.set_ylim(-100, 5)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# 周波数応答（位相）
ax = axes[1, 1]
for name, w in windows.items():
    h = h_ideal * w
    freq, H = signal.freqz(h, worN=4096, fs=fs)
    phase = np.unwrap(np.angle(H))
    ax.plot(freq, phase * 180 / np.pi, linewidth=1.2, label=name, alpha=0.8)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Phase [degrees]')
ax.set_title('Frequency Response (Phase)')
ax.set_xlim(0, fs/2)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fir_window_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

この4つのグラフから、窓関数法による FIR フィルタ設計の核心が読み取れます。

左上（窓関数の波形）: 矩形窓は全区間で値が1で、端部での不連続（打ち切り）が大きいです。ハミング窓とハニング窓は両端でスムーズにゼロに近づき、ブラックマン窓はさらに滑らかです。窓の端部でのスムーズさが、周波数領域でのサイドローブ抑圧に直結しています。
右上（インパルス応答）: 理想 sinc 関数を各窓関数で切り出した結果です。矩形窓は sinc をそのまま切り取るため端部の値が大きく残り、これがサイドローブの原因になります。ブラックマン窓では端部が強く抑圧され、sinc の「裾」が切り詰められています。
左下（振幅応答）: 最も重要なグラフです。矩形窓は遷移帯域幅が狭い（急峻な遮断）反面、阻止域で $-21$ dB 程度の大きなサイドローブが観測されます。ハミング窓は阻止域で $-53$ dB 程度の抑圧を実現していますが、遷移帯域はやや広がっています。ブラックマン窓は阻止域で $-74$ dB 以上の優れた抑圧を示しますが、遷移帯域がさらに広くなっています。遷移帯域幅とサイドローブ抑圧がトレードオフの関係にあることが明確に確認できます。
右下（位相応答）: すべての窓関数で位相が直線的に変化しています。これは FIR フィルタの対称係数による厳密な直線位相特性で、フィルタ長 $M = 51$ に対して群遅延は $(M-1)/2 = 25$ サンプルで一定です。

フィルタ次数の効果

次に、窓関数を固定してフィルタ次数 $M$ を変化させた場合の効果を確認します。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000
fc = 100
wc = 2 * np.pi * fc / fs

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 異なるフィルタ次数で比較（ハミング窓）
orders = [11, 21, 51, 101, 201]

ax = axes[0]
for M in orders:
    n = np.arange(M)
    n_center = (M - 1) / 2
    h_ideal = np.where(n == n_center,
                       wc / np.pi,
                       np.sin(wc * (n - n_center)) / (np.pi * (n - n_center)))
    h = h_ideal * np.hamming(M)
    freq, H = signal.freqz(h, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H) + 1e-15), linewidth=1.5, label=f'M={M}')

ax.axvline(x=fc, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('FIR LPF: Effect of Filter Order (Hamming window)')
ax.set_xlim(0, fs/2)
ax.set_ylim(-100, 5)
ax.legend()
ax.grid(True, alpha=0.3)

# 群遅延
ax = axes[1]
for M in orders:
    n = np.arange(M)
    n_center = (M - 1) / 2
    h_ideal = np.where(n == n_center,
                       wc / np.pi,
                       np.sin(wc * (n - n_center)) / (np.pi * (n - n_center)))
    h = h_ideal * np.hamming(M)
    freq, gd = signal.group_delay((h, 1), w=4096, fs=fs)
    ax.plot(freq, gd, linewidth=1.5, label=f'M={M}, delay={(M-1)/2:.0f}')

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Group delay [samples]')
ax.set_title('FIR LPF: Group Delay')
ax.set_xlim(0, fs/2)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fir_order_effect.png', dpi=150, bbox_inches='tight')
plt.show()

この結果から、フィルタ次数の影響が明確に読み取れます。

左図（振幅応答）: フィルタ次数 $M$ を増やすほど、遷移帯域幅が狭くなり、理想的な矩形特性に近づいていきます。$M = 11$ では遷移が非常に緩やかで、$M = 201$ ではかなり急峻な遮断が実現されています。阻止域減衰はハミング窓の特性（約 $-53$ dB）で決まり、$M$ を変えてもほぼ一定です。これは窓関数法の重要な性質で、遷移帯域幅はフィルタ次数で制御し、サイドローブ抑圧は窓関数の種類で制御するという明確な分離が成り立っています。
右図（群遅延）: 群遅延は全周波数でほぼ一定であり、$(M-1)/2$ サンプルに等しいことが確認できます。$M = 201$ では群遅延が 100 サンプル（$0.1$ 秒 @ 1000 Hz）になります。急峻なフィルタ特性を求めるほど遅延が大きくなるというトレードオフが存在します。リアルタイム性が重要な応用では、この遅延が制約になります。

Python での IIR フィルタ設計

バターワースフィルタとチェビシェフフィルタの比較

SciPy を使って IIR フィルタを設計し、バターワースとチェビシェフI型の特性を比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000
fc = 100
wn = fc / (fs / 2)  # 正規化カットオフ周波数（ナイキスト周波数に対する比）

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# バターワースフィルタ（次数を変化）
ax = axes[0, 0]
for N in [2, 4, 6, 8, 10]:
    b, a = signal.butter(N, wn, btype='low')
    freq, H = signal.freqz(b, a, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H) + 1e-15), linewidth=1.5, label=f'N={N}')
ax.axvline(x=fc, color='gray', linestyle='--', alpha=0.5)
ax.axhline(y=-3, color='gray', linestyle=':', alpha=0.5, label='-3 dB')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('Butterworth LPF')
ax.set_xlim(0, fs/2)
ax.set_ylim(-80, 5)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# チェビシェフI型フィルタ（次数を変化、リップル1dB）
ax = axes[0, 1]
rp = 1  # 通過帯域リップル [dB]
for N in [2, 4, 6, 8, 10]:
    b, a = signal.cheby1(N, rp, wn, btype='low')
    freq, H = signal.freqz(b, a, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H) + 1e-15), linewidth=1.5, label=f'N={N}')
ax.axvline(x=fc, color='gray', linestyle='--', alpha=0.5)
ax.axhline(y=-rp, color='gray', linestyle=':', alpha=0.5, label=f'-{rp} dB ripple')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title(f'Chebyshev Type I LPF (Rp={rp} dB)')
ax.set_xlim(0, fs/2)
ax.set_ylim(-80, 5)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# バターワース vs チェビシェフ（同次数N=4）
ax = axes[1, 0]
N = 4
b_butter, a_butter = signal.butter(N, wn, btype='low')
b_cheby1, a_cheby1 = signal.cheby1(N, 1, wn, btype='low')
b_cheby2, a_cheby2 = signal.cheby2(N, 40, wn, btype='low')

for (b, a, name) in [(b_butter, a_butter, 'Butterworth'),
                      (b_cheby1, a_cheby1, 'Chebyshev I (1dB)'),
                      (b_cheby2, a_cheby2, 'Chebyshev II (40dB)')]:
    freq, H = signal.freqz(b, a, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H) + 1e-15), linewidth=1.5, label=name)
ax.axvline(x=fc, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title(f'Filter Type Comparison (N={N})')
ax.set_xlim(0, fs/2)
ax.set_ylim(-80, 5)
ax.legend()
ax.grid(True, alpha=0.3)

# 群遅延の比較
ax = axes[1, 1]
for (b, a, name) in [(b_butter, a_butter, 'Butterworth'),
                      (b_cheby1, a_cheby1, 'Chebyshev I'),
                      (b_cheby2, a_cheby2, 'Chebyshev II')]:
    freq, gd = signal.group_delay((b, a), w=4096, fs=fs)
    ax.plot(freq, gd, linewidth=1.5, label=name)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Group delay [samples]')
ax.set_title(f'Group Delay Comparison (N={N})')
ax.set_xlim(0, fs/2)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('iir_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

この比較から、IIR フィルタの特性について重要な知見が得られます。

左上（バターワース）: 次数 $N$ を増やすほど遮断が急峻になり、$N = 10$ では非常に急峻な特性が実現されています。通過帯域は完全に平坦で、全ての次数で $-3$ dB 点がカットオフ周波数 $100$ Hz に正確に一致しています。これがバターワースフィルタの最大平坦特性です。
右上（チェビシェフI型）: 同じ次数 $N$ でもバターワースより急峻な遮断が実現されていますが、通過帯域に $1$ dB のリップルが観察されます。次数が上がるとリップルの振動回数が増えますが、各リップルの振幅は $1$ dB（設定値）で等しくなっています。この等リップル性がチェビシェフ多項式の性質を反映しています。
左下（3種類の比較）: 同じ4次で比較すると、チェビシェフI型が最も急峻な遮断を示しています。バターワースは通過帯域が完全に平坦ですが遮断は最も緩やかです。チェビシェフII型は通過帯域は平坦で阻止域にリップルがあり、バターワースとチェビシェフI型の中間的な特性です。
右下（群遅延）: IIR フィルタの群遅延は周波数によって変化しています。特にカットオフ周波数近傍で群遅延が大きくなるピークが見られます。チェビシェフI型はバターワースよりもカットオフ付近の群遅延変動が大きく、位相特性の非線形性がより顕著です。FIR フィルタの一定群遅延とは対照的であり、位相歪みが問題になる応用では注意が必要です。

零点-極点配置の可視化

FIR フィルタと IIR フィルタの零点-極点配置をz平面上で比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))

# FIRフィルタ（ハミング窓、M=21）
M_fir = 21
fs = 1000
fc = 100
wc = 2 * np.pi * fc / fs
n = np.arange(M_fir)
n_center = (M_fir - 1) / 2
h_ideal = np.where(n == n_center, wc / np.pi,
                   np.sin(wc * (n - n_center)) / (np.pi * (n - n_center)))
h_fir = h_ideal * np.hamming(M_fir)

ax = axes[0]
zeros_fir, poles_fir, _ = signal.tf2zpk(h_fir, [1])
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5, alpha=0.5)
ax.plot(np.real(zeros_fir), np.imag(zeros_fir), 'bo', markersize=6, label='Zeros')
ax.plot(np.real(poles_fir), np.imag(poles_fir), 'rx', markersize=8,
        markeredgewidth=2, label='Poles')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title(f'FIR LPF (M={M_fir}, Hamming)')
ax.set_aspect('equal')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-2.5, 2.5)

# バターワースIIRフィルタ（N=4）
wn = fc / (fs/2)
b_iir, a_iir = signal.butter(4, wn, btype='low')

ax = axes[1]
zeros_iir, poles_iir, _ = signal.tf2zpk(b_iir, a_iir)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5, alpha=0.5)
ax.plot(np.real(zeros_iir), np.imag(zeros_iir), 'bo', markersize=6, label='Zeros')
ax.plot(np.real(poles_iir), np.imag(poles_iir), 'rx', markersize=8,
        markeredgewidth=2, label='Poles')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('Butterworth IIR LPF (N=4)')
ax.set_aspect('equal')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)

# チェビシェフI型IIRフィルタ（N=4）
b_cheby, a_cheby = signal.cheby1(4, 1, wn, btype='low')

ax = axes[2]
zeros_cheby, poles_cheby, _ = signal.tf2zpk(b_cheby, a_cheby)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5, alpha=0.5)
ax.plot(np.real(zeros_cheby), np.imag(zeros_cheby), 'bo', markersize=6, label='Zeros')
ax.plot(np.real(poles_cheby), np.imag(poles_cheby), 'rx', markersize=8,
        markeredgewidth=2, label='Poles')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('Chebyshev I IIR LPF (N=4, Rp=1dB)')
ax.set_aspect('equal')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)

plt.tight_layout()
plt.savefig('zpk_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

零点-極点配置から、FIR と IIR の構造的な違いが一目瞭然です。

FIR（左図）: 極は全て原点（$z = 0$）に集中しており、$M$ 重極です。単位円の内部にあるため安定性は自動的に保証されます。零点は単位円の内外に分布しており、阻止域を実現するために単位円の $z = -1$（$\omega = \pi$、ナイキスト周波数）付近に零点が密集しています。零点の分布パターンがフィルタの周波数応答を決定しています。
バターワース IIR（中央図）: 極が単位円の内部に配置され、カットオフ周波数に対応する角度付近に位置しています。極が単位円に近いため、カットオフ付近で急峻な遮断が実現されています。零点は $z = -1$ に集中しており、ナイキスト周波数で完全な遮断を保証しています。極と零点の数が少ないことが、IIR フィルタの計算効率の良さを反映しています。
チェビシェフI型 IIR（右図）: バターワースと同様の構造ですが、極の配置が楕円上にあり（バターワースは円上）、単位円により近い極が存在します。この近接した極が、通過帯域のリップルと急峻な遮断特性の両方を生み出しています。

ノイズ除去デモ

実践的なフィルタリング

最後に、FIR と IIR フィルタを使った実践的なノイズ除去デモを行います。低周波の信号に高周波ノイズが重畳した合成信号をフィルタリングし、元の信号を復元します。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 合成信号の生成
fs = 1000       # サンプリング周波数 [Hz]
T = 0.5         # 信号長 [s]
t = np.arange(0, T, 1/fs)
n_samples = len(t)

# 元の信号: 10Hzと30Hzの正弦波の合成
signal_clean = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 30 * t)

# ノイズ: 高周波ノイズ（200Hzと350Hzの正弦波 + ランダムノイズ）
np.random.seed(0)
noise = 0.8 * np.sin(2 * np.pi * 200 * t) + 0.3 * np.sin(2 * np.pi * 350 * t)
noise += 0.5 * np.random.randn(n_samples)

# ノイズ混入信号
signal_noisy = signal_clean + noise

# FIRフィルタ（ハミング窓、カットオフ50Hz）
M_fir = 101
wn_fir = 50 / (fs / 2)
h_fir = signal.firwin(M_fir, wn_fir, window='hamming')

# IIRフィルタ（バターワース4次、カットオフ50Hz）
b_iir, a_iir = signal.butter(4, wn_fir, btype='low')

# フィルタリング
signal_fir = signal.lfilter(h_fir, 1, signal_noisy)
signal_iir = signal.lfilter(b_iir, a_iir, signal_noisy)

# 零位相フィルタリング（前方向+後方向、IIR）
signal_iir_zero = signal.filtfilt(b_iir, a_iir, signal_noisy)

fig, axes = plt.subplots(4, 1, figsize=(14, 14))

# 元の信号とノイズ混入信号
ax = axes[0]
ax.plot(t * 1000, signal_clean, 'b-', linewidth=1.5, label='Clean signal', alpha=0.8)
ax.plot(t * 1000, signal_noisy, 'r-', linewidth=0.5, alpha=0.5, label='Noisy signal')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('Original and Noisy Signals')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 200)

# FIRフィルタの結果
ax = axes[1]
ax.plot(t * 1000, signal_clean, 'b-', linewidth=1.5, label='Clean', alpha=0.5)
ax.plot(t * 1000, signal_fir, 'g-', linewidth=1.5, label=f'FIR (M={M_fir})')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('FIR Filtered Signal')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 200)

# IIRフィルタの結果（通常 vs 零位相）
ax = axes[2]
ax.plot(t * 1000, signal_clean, 'b-', linewidth=1.5, label='Clean', alpha=0.5)
ax.plot(t * 1000, signal_iir, 'm-', linewidth=1.5, label='IIR (forward only)')
ax.plot(t * 1000, signal_iir_zero, 'g-', linewidth=1.5, label='IIR (zero-phase, filtfilt)')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('IIR Filtered Signal')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 200)

# スペクトル比較
ax = axes[3]
freqs = np.fft.rfftfreq(n_samples, 1/fs)
spec_clean = np.abs(np.fft.rfft(signal_clean)) / n_samples
spec_noisy = np.abs(np.fft.rfft(signal_noisy)) / n_samples
spec_fir = np.abs(np.fft.rfft(signal_fir)) / n_samples
spec_iir = np.abs(np.fft.rfft(signal_iir_zero)) / n_samples

ax.plot(freqs, spec_noisy, 'r-', linewidth=0.8, alpha=0.5, label='Noisy')
ax.plot(freqs, spec_clean, 'b-', linewidth=1.5, label='Clean', alpha=0.7)
ax.plot(freqs, spec_fir, 'g--', linewidth=1.5, label='FIR filtered')
ax.plot(freqs, spec_iir, 'm:', linewidth=2, label='IIR filtered (zero-phase)')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude')
ax.set_title('Frequency Spectrum Comparison')
ax.set_xlim(0, fs/2)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('noise_removal_demo.png', dpi=150, bbox_inches='tight')
plt.show()

このノイズ除去デモから、FIR と IIR フィルタの実践的な振る舞いが明確に読み取れます。

最上段（元信号 vs ノイズ混入信号）: 10 Hz と 30 Hz の低周波信号に、200 Hz と 350 Hz の高周波ノイズとランダムノイズが重畳しています。ノイズ混入信号（赤）からは元の波形パターンをほとんど識別できません。
2段目（FIR フィルタ結果）: 101次のFIRフィルタがノイズを効果的に除去し、元の信号をよく再現しています。ただし、フィルタの群遅延（$(M-1)/2 = 50$ サンプル = 50 ms）分だけ出力が遅延していることに注目してください。この遅延は FIR フィルタの本質的な特性であり、リアルタイム処理では考慮が必要です。
3段目（IIR フィルタ結果）: 前方向のみのフィルタリング（マゼンタ）は位相歪みにより波形が若干変形しています。一方、filtfiltによる零位相フィルタリング（緑）は、フィルタを前方向と後方向の2回適用することで位相歪みを完全に除去し、遅延もゼロです。ただし filtfilt はオフライン処理（全データが手元にある場合）にしか使えません。
最下段（スペクトル比較）: ノイズ成分（200 Hz, 350 Hz, ランダム雑音）がフィルタリング後にほぼ完全に除去されています。FIR と IIR の両方とも、目的の信号成分（10 Hz, 30 Hz）はほぼ無損失で通過させ、カットオフ（50 Hz）より上の成分を効果的に除去しています。

SOS 形式による安定な IIR フィルタ実装

高次の IIR フィルタを直接形で実装すると、係数の数値精度の問題で不安定になることがあります。この問題を解決する標準的な手法が、2次セクション（SOS）形式 です。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000
fc = 100
wn = fc / (fs / 2)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 高次バターワースフィルタ（直接形 vs SOS形式）
orders = [4, 8, 12, 16]

ax = axes[0]
for N in orders:
    # 直接形
    b, a = signal.butter(N, wn, btype='low')
    freq, H_direct = signal.freqz(b, a, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H_direct) + 1e-15),
            linewidth=1, alpha=0.6, linestyle='--')

    # SOS形式
    sos = signal.butter(N, wn, btype='low', output='sos')
    freq, H_sos = signal.sosfreqz(sos, worN=4096, fs=fs)
    ax.plot(freq, 20 * np.log10(np.abs(H_sos) + 1e-15),
            linewidth=1.5, label=f'N={N} (SOS)')

ax.axvline(x=fc, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('Direct Form (dashed) vs SOS (solid)')
ax.set_xlim(0, fs/2)
ax.set_ylim(-120, 5)
ax.legend()
ax.grid(True, alpha=0.3)

# 高次フィルタの極配置比較
ax = axes[1]
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5, alpha=0.5)

colors = ['blue', 'red', 'green', 'purple']
for N, color in zip(orders, colors):
    b, a = signal.butter(N, wn, btype='low')
    z, p, k = signal.tf2zpk(b, a)
    ax.plot(np.real(p), np.imag(p), 'x', color=color, markersize=8,
            markeredgewidth=2, label=f'N={N} poles')

ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('Pole Locations (Higher Order = Closer to Unit Circle)')
ax.set_aspect('equal')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-1.2, 1.2)

plt.tight_layout()
plt.savefig('sos_implementation.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、高次 IIR フィルタの実装上の注意点が読み取れます。

左図（直接形 vs SOS）: 低次（$N = 4, 8$）では直接形（破線）と SOS 形式（実線）はほぼ同じ結果を示しています。しかし高次（$N = 12, 16$）では、直接形の周波数応答に不自然な振動や数値的な乱れが現れる場合があります。SOS 形式はどの次数でも安定で滑らかな特性を維持しています。実用上は常に SOS 形式を使うことが推奨されます。
右図（極の配置）: 次数が上がるほど極が単位円に接近しています。$N = 16$ の極は単位円のごく近傍にあり、わずかな係数の丸め誤差で極が単位円外に出る（不安定化する）リスクがあります。SOS 形式では各2次セクションの極が個別に管理されるため、この問題が緩和されます。

まとめ

本記事では、ディジタルフィルタの2大分類である FIR フィルタと IIR フィルタについて、理論から設計、実装まで体系的に解説しました。

ディジタルフィルタの基礎: 差分方程式 $y[n] = \sum b_k x[n-k] – \sum a_k y[n-k]$ と伝達関数 $H(z) = B(z)/A(z)$ が全ての出発点。零点と極がフィルタの周波数応答を幾何学的に決定する
FIR フィルタ: フィードバックなし（$a_k = 0$）で常に安定。対称係数により厳密な直線位相を実現。窓関数法では、遷移帯域幅をフィルタ次数 $M$ で、サイドローブ抑圧を窓関数の種類で制御する
IIR フィルタ: フィードバックにより少ない次数で急峻な遮断を実現。バターワース（最大平坦）、チェビシェフ（等リップル）のアナログプロトタイプを双一次変換 $s = (2/T)(1-z^{-1})/(1+z^{-1})$ でディジタル化する
FIR vs IIR: 直線位相・安定性保証なら FIR、計算効率・低遅延なら IIR。用途に応じた選択が重要
実装の注意点: 高次 IIR は SOS 形式で実装し、係数量子化による不安定化を防ぐ。零位相フィルタリング（filtfilt）によりオフライン処理で位相歪みを除去できる

ディジタルフィルタは信号処理の基礎中の基礎であり、通信、音声、制御、画像処理など、あらゆるディジタルシステムで活用されています。本記事で学んだ FIR/IIR の理論と設計法は、より高度な適応フィルタ（LMS, RLS）、多レートフィルタ（デシメーション、インターポレーション）、ウェーブレット変換などへの入口となります。

次のステップとして、以下の記事も参考にしてください。