ディジタルフィルタは、離散時間信号の周波数成分を選択的に通過・遮断する処理であり、音声処理、通信、画像処理、制御システムなど、あらゆるディジタル信号処理の根幹をなす技術です。
ディジタルフィルタは大きく FIR（Finite Impulse Response）フィルタと IIR（Infinite Impulse Response）フィルタに分類されます。本記事では、差分方程式とz変換から出発して、両者の特徴、周波数応答の導出、安定性条件、線形位相条件を数学的に導出し、Pythonによる設計と可視化を行います。
本記事の内容
- ディジタルフィルタの定義（差分方程式）
- z変換とシステム関数 $H(z)$ の導出
- FIRフィルタの特徴（有限インパルス応答、常に安定、線形位相可能）
- IIRフィルタの特徴（無限インパルス応答、低次で急峻な遮断）
- 周波数応答 $H(e^{j\omega})$ の導出
- 零点と極の配置と周波数特性の関係
- 安定性条件（極が単位円内）の証明
- 線形位相条件の導出
- Pythonでのフィルタ設計と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ディジタルフィルタの定義
差分方程式
ディジタルフィルタは、入力信号 $x[n]$ と出力信号 $y[n]$ の関係を定差分方程式で記述します。最も一般的な形は
$$ \sum_{k=0}^{N} a_k \, y[n-k] = \sum_{m=0}^{M} b_m \, x[n-m] $$
です。$a_0 = 1$ と正規化すると
$$ y[n] = \sum_{m=0}^{M} b_m \, x[n-m] – \sum_{k=1}^{N} a_k \, y[n-k] $$
右辺の第1項は入力信号の現在値と過去値の重み付き和、第2項は出力信号の過去値のフィードバックです。
LTIシステムとしてのフィルタ
ディジタルフィルタは線形時不変（LTI: Linear Time-Invariant）システムです。LTIシステムは以下の性質を持ちます。
線形性: 入力 $\alpha x_1[n] + \beta x_2[n]$ に対する出力は $\alpha y_1[n] + \beta y_2[n]$
時不変性: 入力を $x[n-n_0]$ とシフトすると、出力も $y[n-n_0]$ とシフトする
LTIシステムはインパルス応答 $h[n]$ で完全に特徴づけられ、出力は畳み込みで得られます。
$$ y[n] = \sum_{k=-\infty}^{\infty} h[k] \, x[n-k] = h[n] * x[n] $$
z変換とシステム関数
z変換の定義
離散時間信号 $x[n]$ のz変換は
$$ X(z) = \sum_{n=-\infty}^{\infty} x[n] \, z^{-n} $$
と定義されます。ここで $z$ は複素変数です。
基本性質
時間シフト性質: $x[n-k]$ のz変換は
$$ \mathcal{Z}\{x[n-k]\} = z^{-k} X(z) $$
証明:
$$ \begin{align} \mathcal{Z}\{x[n-k]\} &= \sum_{n=-\infty}^{\infty} x[n-k] \, z^{-n} \\ &= \sum_{m=-\infty}^{\infty} x[m] \, z^{-(m+k)} \quad (m = n-k) \\ &= z^{-k} \sum_{m=-\infty}^{\infty} x[m] \, z^{-m} \\ &= z^{-k} X(z) \end{align} $$
$\square$
畳み込み定理: $y[n] = h[n] * x[n]$ ならば
$$ Y(z) = H(z) \cdot X(z) $$
伝達関数（システム関数）の導出
差分方程式
$$ y[n] + \sum_{k=1}^{N} a_k \, y[n-k] = \sum_{m=0}^{M} b_m \, x[n-m] $$
の両辺をz変換します。時間シフト性質を用いると
$$ Y(z) + \sum_{k=1}^{N} a_k z^{-k} Y(z) = \sum_{m=0}^{M} b_m z^{-m} X(z) $$
$$ Y(z) \left(1 + \sum_{k=1}^{N} a_k z^{-k}\right) = X(z) \sum_{m=0}^{M} b_m z^{-m} $$
伝達関数（システム関数）$H(z)$ は
$$ \begin{align} H(z) &= \frac{Y(z)}{X(z)} \\ &= \frac{\sum_{m=0}^{M} b_m z^{-m}}{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}} \end{align} $$
分子・分母に $z^{\max(M,N)}$ を掛けて整理すると
$$ H(z) = \frac{b_0 z^N + b_1 z^{N-1} + \cdots + b_M z^{N-M}}{z^N + a_1 z^{N-1} + \cdots + a_N} $$
零点と極
$H(z)$ の分子 = 0 の根を零点（zeros）$z_1, z_2, \dots, z_M$、分母 = 0 の根を極（poles）$p_1, p_2, \dots, p_N$ と呼びます。$H(z)$ は零点と極を用いて
$$ H(z) = b_0 \cdot \frac{(1 – z_1 z^{-1})(1 – z_2 z^{-1}) \cdots (1 – z_M z^{-1})}{(1 – p_1 z^{-1})(1 – p_2 z^{-1}) \cdots (1 – p_N z^{-1})} $$
と因数分解できます。零点と極の配置がフィルタの特性を決定します。
FIRフィルタ
定義と特徴
FIR（Finite Impulse Response）フィルタは、フィードバック項を持たないフィルタです。差分方程式は
$$ y[n] = \sum_{m=0}^{M} b_m \, x[n-m] $$
インパルス応答は有限長で
$$ h[n] = \begin{cases} b_n & 0 \leq n \leq M \\ 0 & \text{otherwise} \end{cases} $$
フィルタ次数は $M$、タップ数（係数の数）は $M + 1$ です。
伝達関数
FIRフィルタの伝達関数は
$$ H(z) = \sum_{m=0}^{M} b_m z^{-m} = b_0 + b_1 z^{-1} + \cdots + b_M z^{-M} $$
分母は1（定数）なので、$M$ 個の零点と $z = 0$ に $M$ 重の極を持ちます。
FIRフィルタの利点
- 常に安定: 極が全て $z = 0$（単位円内）にあるため、常にBIBO安定（後述）
- 線形位相が実現可能: 係数の対称性・反対称性により、正確な線形位相が得られる
- 設計が容易: 窓関数法、周波数サンプリング法、最適設計法（Parks-McClellan）など
FIRフィルタの欠点
- 高次数が必要: 急峻な遮断特性を得るにはフィルタ次数が大きくなる（計算量が増加）
- 遅延が大きい: 線形位相FIRの群遅延は $M/2$ サンプル
IIRフィルタ
定義と特徴
IIR（Infinite Impulse Response）フィルタは、フィードバック項を持つフィルタです。差分方程式は
$$ y[n] = \sum_{m=0}^{M} b_m \, x[n-m] – \sum_{k=1}^{N} a_k \, y[n-k] $$
フィードバック（再帰）があるため、有限長の入力に対してもインパルス応答は無限に続きます。
具体例: 1次IIRフィルタ
最も単純な1次IIRフィルタを考えます。
$$ y[n] = x[n] + a \, y[n-1], \quad |a| < 1 $$
インパルス応答を求めます。$x[n] = \delta[n]$ を入力すると
$$ \begin{align} y[0] &= \delta[0] + a \cdot 0 = 1 \\ y[1] &= \delta[1] + a \cdot y[0] = a \\ y[2] &= \delta[2] + a \cdot y[1] = a^2 \\ &\vdots \\ y[n] &= a^n, \quad n \geq 0 \end{align} $$
したがって $h[n] = a^n u[n]$（$u[n]$ は単位ステップ関数）です。$|a| < 1$ なら $h[n] \to 0$ ($n \to \infty$) で安定ですが、インパルス応答は厳密にはゼロにはならず無限に続きます。
伝達関数
$$ H(z) = \frac{1}{1 – a z^{-1}} = \frac{z}{z – a} $$
極は $z = a$ にあります。
IIRフィルタの利点
- 低次で急峻な遮断: フィードバックにより、少ない係数で急峻な周波数特性を実現
- アナログフィルタからの変換: バターワース、チェビシェフ、楕円フィルタなどの設計法が確立
IIRフィルタの欠点
- 不安定になりうる: 極の配置によっては不安定（発振）
- 正確な線形位相は不可能: IIRフィルタでは厳密な線形位相を実現できない
- 有限精度の影響: 固定小数点実装でのリミットサイクルなど
周波数応答の導出
$H(e^{j\omega})$ の導出
ディジタルフィルタの周波数応答は、伝達関数 $H(z)$ の $z = e^{j\omega}$（単位円上）での値として得られます。
$$ H(e^{j\omega}) = H(z)\big|_{z = e^{j\omega}} = \frac{\sum_{m=0}^{M} b_m e^{-j\omega m}}{1 + \sum_{k=1}^{N} a_k e^{-j\omega k}} $$
ここで $\omega$ は正規化角周波数（$\omega = 2\pi f / f_s$、範囲 $[0, 2\pi]$ または $[-\pi, \pi]$）です。
なぜ $z = e^{j\omega}$ なのか: z変換はフーリエ変換の一般化であり、$z = re^{j\omega}$ と極座標表示したとき $r = 1$（単位円上）でDTFT（離散時間フーリエ変換）に一致します。
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} = X(z)\big|_{z=e^{j\omega}} $$
振幅応答と位相応答
周波数応答を極形式で表すと
$$ H(e^{j\omega}) = |H(e^{j\omega})| \, e^{j\angle H(e^{j\omega})} $$
- 振幅応答（magnitude response）: $|H(e^{j\omega})|$ — 各周波数の振幅がどれだけ増幅/減衰されるか
- 位相応答（phase response）: $\angle H(e^{j\omega})$ — 各周波数の位相がどれだけ変化するか
デシベル表記の振幅応答は
$$ |H(e^{j\omega})|_{\text{dB}} = 20 \log_{10} |H(e^{j\omega})| $$
群遅延
群遅延は位相応答の微分として定義されます。
$$ \tau_g(\omega) = -\frac{d\angle H(e^{j\omega})}{d\omega} $$
群遅延は「周波数 $\omega$ の成分がフィルタを通過する際に受ける時間遅延」を表します。線形位相フィルタでは群遅延が一定（全周波数で同じ遅延）です。
零点・極の配置と周波数特性
零点と極の周波数特性への影響
周波数 $\omega$ における振幅応答は、零点と極からの距離で幾何学的に理解できます。
$$ |H(e^{j\omega})| = |b_0| \cdot \frac{\prod_{m=1}^{M} |e^{j\omega} – z_m|}{\prod_{k=1}^{N} |e^{j\omega} – p_k|} $$
ここで $|e^{j\omega} – z_m|$ は単位円上の点 $e^{j\omega}$ から零点 $z_m$ までの距離、$|e^{j\omega} – p_k|$ は極 $p_k$ までの距離です。
この幾何学的解釈から
- 零点が単位円に近い → その角周波数付近で振幅応答が小さくなる（ノッチ）
- 零点が単位円上 → その角周波数で振幅応答がゼロ（完全な遮断）
- 極が単位円に近い → その角周波数付近で振幅応答が大きくなる（ピーク）
- 極が単位円上 → その角周波数で振幅応答が無限大（不安定）
設計の指針
- 低域通過フィルタ: $\omega = \pi$（ナイキスト周波数）付近に零点を配置し、$\omega = 0$ 付近に極を配置
- 高域通過フィルタ: $\omega = 0$ 付近に零点を配置し、$\omega = \pi$ 付近に極を配置
- 帯域通過フィルタ: 通過帯域外に零点、通過帯域内に極を配置
- ノッチフィルタ: 除去したい周波数の単位円上に零点を配置
安定性条件
BIBO安定性
BIBO（Bounded-Input Bounded-Output）安定性とは、有界な入力に対して常に有界な出力が得られることです。
定義: 任意の入力 $|x[n]| \leq B_x < \infty$ に対して、$|y[n]| \leq B_y < \infty$ が成り立つとき、システムはBIBO安定である。
安定性の必要十分条件（インパルス応答）
LTIシステムがBIBO安定であるための必要十分条件は
$$ \sum_{n=-\infty}^{\infty} |h[n]| < \infty $$
すなわち、インパルス応答が絶対総和可能であることです。
証明（十分性）: $|y[n]| = \left|\sum_k h[k] x[n-k]\right| \leq \sum_k |h[k]| |x[n-k]| \leq B_x \sum_k |h[k]| < \infty$
証明（必要性）: $x[n] = \text{sgn}(h[-n])$ とおくと、$|x[n]| \leq 1$ かつ $y[0] = \sum_k |h[k]|$。$y[0]$ が有界であるためには $\sum_k |h[k]| < \infty$ が必要。$\square$
安定性の条件（極による判定）
定理: 因果的LTIシステムがBIBO安定であるための必要十分条件は、$H(z)$ の全ての極が単位円の内部にあることです。
$$ |p_k| < 1, \quad k = 1, 2, \dots, N $$
証明: 因果的システムの場合、$h[n] = 0$ ($n < 0$) です。$H(z)$ の極を $p_1, \dots, p_N$ とすると、部分分数展開により
$$ h[n] = \sum_{k=1}^{N} c_k \, p_k^n \, u[n] $$
（ここでは簡単のため、極が全て単純（1位）の場合を考えます）
$\sum_{n=0}^{\infty} |h[n]| < \infty$ のためには
$$ \sum_{n=0}^{\infty} |h[n]| \leq \sum_{k=1}^{N} |c_k| \sum_{n=0}^{\infty} |p_k|^n = \sum_{k=1}^{N} \frac{|c_k|}{1 – |p_k|} $$
この等比級数が収束するための条件は $|p_k| < 1$（全ての $k$）です。
逆に、ある $|p_k| \geq 1$ ならば $|p_k|^n \to \infty$ または $|p_k|^n = 1$ となり、$h[n]$ は絶対総和可能ではありません。$\square$
FIRフィルタの安定性
FIRフィルタの極は全て $z = 0$ にあります。$|0| = 0 < 1$ なので、FIRフィルタは常にBIBO安定です。これがFIRフィルタの重要な利点です。
線形位相条件
線形位相の定義と重要性
位相応答が周波数に対して線形であるとき、すなわち
$$ \angle H(e^{j\omega}) = -\alpha \omega + \beta $$
の形であるとき、群遅延は
$$ \tau_g(\omega) = -\frac{d(-\alpha \omega + \beta)}{d\omega} = \alpha $$
一定値になります。線形位相は、信号波形を歪ませずに一定時間だけ遅延させるため、音声・画像処理で重要です。
FIRフィルタの線形位相条件
定理: 長さ $M+1$ のFIRフィルタ $h[n]$（$n = 0, 1, \dots, M$）が線形位相を持つための必要十分条件は、次のいずれかの対称条件を満たすことです。
Type I: $h[n] = h[M-n]$（偶対称、$M$ が偶数）
Type II: $h[n] = h[M-n]$（偶対称、$M$ が奇数）
Type III: $h[n] = -h[M-n]$（奇対称、$M$ が偶数）
Type IV: $h[n] = -h[M-n]$（奇対称、$M$ が奇数）
証明（Type I の場合）
$M$ が偶数で $h[n] = h[M-n]$ の場合を証明します。
$$ \begin{align} H(e^{j\omega}) &= \sum_{n=0}^{M} h[n] e^{-j\omega n} \end{align} $$
$n \to M – n$ と変数変換します。
$$ \begin{align} H(e^{j\omega}) &= \sum_{n=0}^{M} h[M-n] e^{-j\omega(M-n)} \\ &= e^{-j\omega M} \sum_{n=0}^{M} h[M-n] e^{j\omega n} \end{align} $$
$h[n] = h[M-n]$ を用いると
$$ H(e^{j\omega}) = e^{-j\omega M} \sum_{n=0}^{M} h[n] e^{j\omega n} = e^{-j\omega M} H^*(e^{j\omega}) $$
最後の等式は $H^*(e^{j\omega}) = \sum_n h[n] e^{j\omega n}$（$h[n]$ が実数）からです。
$H(e^{j\omega}) = |H(e^{j\omega})| e^{j\phi(\omega)}$ とおくと
$$ |H(e^{j\omega})| e^{j\phi(\omega)} = e^{-j\omega M} |H(e^{j\omega})| e^{-j\phi(\omega)} $$
振幅は正なので、位相について
$$ e^{j\phi(\omega)} = e^{-j\omega M} e^{-j\phi(\omega)} $$
$$ e^{j 2\phi(\omega)} = e^{-j\omega M} $$
$$ 2\phi(\omega) = -\omega M + 2k\pi $$
$$ \phi(\omega) = -\frac{M}{2} \omega + k\pi $$
これは $\alpha = M/2$, $\beta = k\pi$ の線形位相です。群遅延は
$$ \tau_g = \frac{M}{2} \quad \text{[samples]} $$
$\square$
各タイプの制約
| タイプ | 対称性 | $M$ | $\omega = 0$ | $\omega = \pi$ | 適用可能なフィルタ |
|---|---|---|---|---|---|
| I | 偶対称 | 偶数 | 制約なし | 制約なし | 全タイプ（LPF, HPF, BPF, BSF） |
| II | 偶対称 | 奇数 | 制約なし | $H(e^{j\pi}) = 0$ | LPF, BPF |
| III | 奇対称 | 偶数 | $H(e^{j0}) = 0$ | $H(e^{j\pi}) = 0$ | BPF, 微分器 |
| IV | 奇対称 | 奇数 | $H(e^{j0}) = 0$ | 制約なし | HPF, BPF, 微分器 |
Type II は $\omega = \pi$ でゼロになるためHPFに使えません。Type III は $\omega = 0$ と $\omega = \pi$ の両方でゼロになるためLPF/HPFに使えません。
Pythonでの実装
FIRフィルタの設計と周波数応答
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== FIRフィルタの設計（窓関数法）=====
fs = 1000 # サンプリング周波数 [Hz]
fc = 100 # カットオフ周波数 [Hz]
num_taps = 51 # タップ数（奇数 → Type I）
# 正規化カットオフ周波数（ナイキスト周波数で正規化）
fc_normalized = fc / (fs / 2)
# 各種窓関数でFIRフィルタを設計
windows = ['rectangular', 'hamming', 'blackman', 'kaiser']
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
for idx, window_name in enumerate(windows):
if window_name == 'rectangular':
window = 'boxcar'
elif window_name == 'kaiser':
window = ('kaiser', 8.0)
else:
window = window_name
# FIRフィルタ設計
h = signal.firwin(num_taps, fc_normalized, window=window)
# 周波数応答
w, H = signal.freqz(h, worN=2048, fs=fs)
ax = axes[idx // 2, idx % 2]
ax.plot(w, 20 * np.log10(np.abs(H) + 1e-12), linewidth=2)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title(f'FIR LPF ({window_name}, {num_taps} taps, $f_c$={fc} Hz)')
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.set_ylim([-80, 5])
ax.axvline(x=fc, color='red', linestyle='--', alpha=0.5, label=f'$f_c$={fc} Hz')
ax.legend()
plt.tight_layout()
plt.savefig('fir_filter_design.png', dpi=150, bbox_inches='tight')
plt.show()
IIRフィルタの設計と周波数応答
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== IIRフィルタの設計（バターワース、チェビシェフ、楕円）=====
fs = 1000 # サンプリング周波数 [Hz]
fc = 100 # カットオフ周波数 [Hz]
order = 4 # フィルタ次数
# 各種IIRフィルタの設計
filters = {}
# バターワースフィルタ
b_butter, a_butter = signal.butter(order, fc, fs=fs)
filters['Butterworth'] = (b_butter, a_butter)
# チェビシェフI型（通過帯域リプル 1 dB）
b_cheby1, a_cheby1 = signal.cheby1(order, 1, fc, fs=fs)
filters['Chebyshev I (1dB)'] = (b_cheby1, a_cheby1)
# チェビシェフII型（阻止帯域減衰 40 dB）
b_cheby2, a_cheby2 = signal.cheby2(order, 40, fc, fs=fs)
filters['Chebyshev II (40dB)'] = (b_cheby2, a_cheby2)
# 楕円フィルタ（通過帯域リプル 1 dB、阻止帯域減衰 40 dB）
b_ellip, a_ellip = signal.ellip(order, 1, 40, fc, fs=fs)
filters['Elliptic (1dB/40dB)'] = (b_ellip, a_ellip)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
colors = ['blue', 'green', 'red', 'purple']
for idx, (name, (b, a)) in enumerate(filters.items()):
# 周波数応答
w, H = signal.freqz(b, a, worN=2048, fs=fs)
ax = axes[idx // 2, idx % 2]
ax.plot(w, 20 * np.log10(np.abs(H) + 1e-12), color=colors[idx], linewidth=2)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title(f'{name} (order={order}, $f_c$={fc} Hz)')
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.set_ylim([-80, 5])
ax.axvline(x=fc, color='red', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.savefig('iir_filter_design.png', dpi=150, bbox_inches='tight')
plt.show()
FIR vs IIR の比較
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== FIR vs IIR の比較（同等の遮断特性）=====
fs = 1000
fc = 100
# 4次バターワースIIR
b_iir, a_iir = signal.butter(4, fc, fs=fs)
w_iir, H_iir = signal.freqz(b_iir, a_iir, worN=4096, fs=fs)
# 同等の阻止帯域減衰を持つFIR（窓関数法）
fir_orders = [15, 31, 63, 127]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 振幅応答の比較
ax = axes[0, 0]
ax.plot(w_iir, 20*np.log10(np.abs(H_iir)+1e-12), 'k-', linewidth=2.5,
label='IIR Butterworth (N=4)')
for n_taps in fir_orders:
h_fir = signal.firwin(n_taps, fc/(fs/2), window='hamming')
w_fir, H_fir = signal.freqz(h_fir, worN=4096, fs=fs)
ax.plot(w_fir, 20*np.log10(np.abs(H_fir)+1e-12), linewidth=1.5,
label=f'FIR Hamming (N={n_taps})')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('FIR vs IIR: Magnitude Response')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.set_ylim([-80, 5])
# (2) 位相応答の比較
ax = axes[0, 1]
# IIR
phase_iir = np.unwrap(np.angle(H_iir))
ax.plot(w_iir, phase_iir * 180 / np.pi, 'k-', linewidth=2.5,
label='IIR Butterworth (N=4)')
# FIR
for n_taps in fir_orders:
h_fir = signal.firwin(n_taps, fc/(fs/2), window='hamming')
w_fir, H_fir = signal.freqz(h_fir, worN=4096, fs=fs)
phase_fir = np.unwrap(np.angle(H_fir))
ax.plot(w_fir, phase_fir * 180 / np.pi, linewidth=1.5,
label=f'FIR (N={n_taps})')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Phase [deg]')
ax.set_title('FIR vs IIR: Phase Response')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
# (3) 群遅延の比較
ax = axes[1, 0]
# IIR
w_gd, gd_iir = signal.group_delay((b_iir, a_iir), w=4096, fs=fs)
ax.plot(w_gd, gd_iir, 'k-', linewidth=2.5, label='IIR Butterworth (N=4)')
# FIR（線形位相 → 一定の群遅延）
for n_taps in fir_orders:
h_fir = signal.firwin(n_taps, fc/(fs/2), window='hamming')
w_gd_fir, gd_fir = signal.group_delay((h_fir, 1), w=4096, fs=fs)
ax.plot(w_gd_fir, gd_fir, linewidth=1.5, label=f'FIR (N={n_taps})')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Group Delay [samples]')
ax.set_title('FIR vs IIR: Group Delay')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.set_ylim([0, 70])
# (4) インパルス応答の比較
ax = axes[1, 1]
# IIR
t_iir, h_imp_iir = signal.dimpulse((b_iir, a_iir, 1/fs), n=100)
ax.stem(t_iir[0].flatten()*fs, h_imp_iir[0].flatten(),
linefmt='k-', markerfmt='ko', basefmt='k-', label='IIR (N=4)')
# FIR (31タップ)
h_fir_31 = signal.firwin(31, fc/(fs/2), window='hamming')
ax.stem(np.arange(31), h_fir_31,
linefmt='b-', markerfmt='bs', basefmt='b-', label='FIR (N=31)')
ax.set_xlabel('Sample index')
ax.set_ylabel('Amplitude')
ax.set_title('FIR vs IIR: Impulse Response')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim([0, 50])
plt.tight_layout()
plt.savefig('fir_vs_iir.png', dpi=150, bbox_inches='tight')
plt.show()
零点・極プロット
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== 零点・極プロット =====
fs = 1000
fc = 100
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
# (1) FIRフィルタ（31タップ、ハミング窓）
ax = axes[0, 0]
h_fir = signal.firwin(31, fc/(fs/2), window='hamming')
z_fir, p_fir, k_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)
ax.scatter(np.real(z_fir), np.imag(z_fir), marker='o', s=60,
facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
ax.scatter(np.real(p_fir), np.imag(p_fir), marker='x', s=60,
color='red', linewidths=2, label='Poles')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('FIR LPF (31 taps, Hamming)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim([-1.5, 1.5])
ax.set_ylim([-1.5, 1.5])
# (2) バターワースIIR（4次）
ax = axes[0, 1]
b_butter, a_butter = signal.butter(4, fc, fs=fs)
z_butter, p_butter, k_butter = signal.tf2zpk(b_butter, a_butter)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5)
ax.scatter(np.real(z_butter), np.imag(z_butter), marker='o', s=60,
facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
ax.scatter(np.real(p_butter), np.imag(p_butter), marker='x', s=60,
color='red', linewidths=2, label='Poles')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('Butterworth IIR LPF (N=4)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim([-1.5, 1.5])
ax.set_ylim([-1.5, 1.5])
# (3) チェビシェフI型IIR（4次）
ax = axes[1, 0]
b_cheby, a_cheby = signal.cheby1(4, 1, fc, fs=fs)
z_cheby, p_cheby, k_cheby = signal.tf2zpk(b_cheby, a_cheby)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5)
ax.scatter(np.real(z_cheby), np.imag(z_cheby), marker='o', s=60,
facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
ax.scatter(np.real(p_cheby), np.imag(p_cheby), marker='x', s=60,
color='red', linewidths=2, label='Poles')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('Chebyshev I IIR LPF (N=4, 1dB ripple)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim([-1.5, 1.5])
ax.set_ylim([-1.5, 1.5])
# (4) 楕円IIR（4次）
ax = axes[1, 1]
b_ellip, a_ellip = signal.ellip(4, 1, 40, fc, fs=fs)
z_ellip, p_ellip, k_ellip = signal.tf2zpk(b_ellip, a_ellip)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5)
ax.scatter(np.real(z_ellip), np.imag(z_ellip), marker='o', s=60,
facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
ax.scatter(np.real(p_ellip), np.imag(p_ellip), marker='x', s=60,
color='red', linewidths=2, label='Poles')
# 極の絶対値を表示
for p in p_ellip:
ax.annotate(f'|p|={np.abs(p):.3f}', xy=(np.real(p), np.imag(p)),
xytext=(10, 10), textcoords='offset points', fontsize=8)
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title('Elliptic IIR LPF (N=4, 1dB/40dB)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim([-1.5, 1.5])
ax.set_ylim([-1.5, 1.5])
plt.tight_layout()
plt.savefig('zpk_plots.png', dpi=150, bbox_inches='tight')
plt.show()
# 安定性の確認
print("=== 極の絶対値（安定性チェック） ===")
for name, poles in [('Butterworth', p_butter), ('Chebyshev I', p_cheby),
('Elliptic', p_ellip)]:
max_pole = np.max(np.abs(poles))
stable = "STABLE" if max_pole < 1 else "UNSTABLE"
print(f"{name}: max|p| = {max_pole:.6f} -> {stable}")
信号のフィルタリング実演
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== 実際の信号にフィルタを適用 =====
fs = 1000 # サンプリング周波数 [Hz]
duration = 0.5 # 信号長 [s]
t = np.arange(0, duration, 1/fs)
# テスト信号: 50Hz + 200Hz + ノイズ
np.random.seed(42)
x_clean_50 = np.sin(2 * np.pi * 50 * t)
x_clean_200 = 0.5 * np.sin(2 * np.pi * 200 * t)
noise = 0.3 * np.random.randn(len(t))
x = x_clean_50 + x_clean_200 + noise
# フィルタ設計: 100Hz LPF（50Hzを通過、200Hzを除去）
# FIRフィルタ（51タップ、ハミング窓）
h_fir = signal.firwin(51, 100/(fs/2), window='hamming')
# IIRフィルタ（4次バターワース）
b_iir, a_iir = signal.butter(4, 100, fs=fs)
# フィルタリング
y_fir = signal.lfilter(h_fir, 1, x)
y_iir = signal.lfilter(b_iir, a_iir, x)
fig, axes = plt.subplots(3, 2, figsize=(14, 12))
# (1) 入力信号（時間波形）
ax = axes[0, 0]
ax.plot(t * 1000, x, 'b-', alpha=0.7, linewidth=0.8)
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('Input Signal (50Hz + 200Hz + noise)')
ax.grid(True, alpha=0.3)
ax.set_xlim([0, 100])
# (2) 入力信号（スペクトル）
ax = axes[0, 1]
f_fft = np.fft.rfftfreq(len(x), 1/fs)
X_fft = np.abs(np.fft.rfft(x)) / len(x) * 2
ax.plot(f_fft, X_fft, 'b-', linewidth=1)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude')
ax.set_title('Input Spectrum')
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.axvline(x=100, color='red', linestyle='--', alpha=0.5, label='$f_c$=100Hz')
ax.legend()
# (3) FIR出力（時間波形）
ax = axes[1, 0]
ax.plot(t * 1000, x_clean_50, 'g--', alpha=0.5, linewidth=1, label='50Hz (ideal)')
ax.plot(t * 1000, y_fir, 'r-', linewidth=1, label='FIR output')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('FIR Filtered Output (51 taps, Hamming)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim([0, 100])
# (4) FIR出力（スペクトル）
ax = axes[1, 1]
Y_fir_fft = np.abs(np.fft.rfft(y_fir)) / len(y_fir) * 2
ax.plot(f_fft, Y_fir_fft, 'r-', linewidth=1)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude')
ax.set_title('FIR Output Spectrum')
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.axvline(x=100, color='red', linestyle='--', alpha=0.5)
# (5) IIR出力（時間波形）
ax = axes[2, 0]
ax.plot(t * 1000, x_clean_50, 'g--', alpha=0.5, linewidth=1, label='50Hz (ideal)')
ax.plot(t * 1000, y_iir, 'm-', linewidth=1, label='IIR output')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('IIR Filtered Output (Butterworth N=4)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim([0, 100])
# (6) IIR出力（スペクトル）
ax = axes[2, 1]
Y_iir_fft = np.abs(np.fft.rfft(y_iir)) / len(y_iir) * 2
ax.plot(f_fft, Y_iir_fft, 'm-', linewidth=1)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude')
ax.set_title('IIR Output Spectrum')
ax.grid(True, alpha=0.3)
ax.set_xlim([0, fs/2])
ax.axvline(x=100, color='red', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.savefig('filter_application.png', dpi=150, bbox_inches='tight')
plt.show()
シミュレーション結果の考察
-
窓関数の影響: 矩形窓は遷移帯域が狭いがサイドローブが大きく（ギブス現象）、ハミング窓やブラックマン窓はサイドローブを抑制する代わりに遷移帯域が広がります。カイザー窓はパラメータ $\beta$ でトレードオフを調整できます。
-
IIRフィルタの種類: バターワースは平坦な通過帯域が特徴で、チェビシェフI型は等リプルの通過帯域で急峻な遮断を実現し、楕円フィルタは通過帯域・阻止帯域ともに等リプルで最も急峻な遮断を実現します。
-
FIR vs IIR: 4次のIIRバターワースフィルタと同等の阻止帯域減衰を得るには、FIRでは63タップ以上が必要です。一方、FIRは正確な線形位相（一定の群遅延）を実現できますが、IIRの群遅延はカットオフ付近でピークを持ちます。
-
零点・極の配置: 全てのIIRフィルタで極が単位円内にあり安定です。楕円フィルタは阻止帯域に零点が単位円上に配置され、特定の周波数で完全な遮断（ノッチ）を実現しています。
-
フィルタリング結果: 50 Hz + 200 Hz + ノイズの混合信号に100 Hz LPFを適用すると、200 Hz成分とノイズの高周波成分が除去され、50 Hz成分が復元されます。FIRとIIRで概ね同等の結果ですが、過渡応答の振る舞いが異なります。
まとめ
本記事では、ディジタルフィルタの基礎（FIR/IIR）について解説しました。
- 差分方程式とz変換: フィルタは差分方程式で記述され、z変換により伝達関数 $H(z) = \frac{\sum b_m z^{-m}}{1 + \sum a_k z^{-k}}$ が導出される
- FIRフィルタ: インパルス応答が有限長。常に安定で線形位相が実現可能だが、急峻な遮断には高次数が必要
- IIRフィルタ: フィードバックにより低次で急峻な遮断が可能だが、不安定になりうる。正確な線形位相は実現不可能
- 周波数応答: $H(e^{j\omega})$ は $H(z)$ を単位円上で評価したもの。振幅応答と位相応答で特性を記述
- 零点と極: 零点は振幅を抑制し、極は振幅を増幅する。幾何学的解釈が可能
- 安定性条件: 因果的LTIシステムのBIBO安定条件は、全ての極が単位円内（$|p_k| < 1$）にあること
- 線形位相: FIRフィルタの係数が対称（$h[n] = h[M-n]$）のとき線形位相が得られ、群遅延は $M/2$ サンプルで一定
次のステップとして、以下の記事も参考にしてください。