IIR（Infinite Impulse Response）フィルタは、フィードバック構造を持つディジタルフィルタです。FIRフィルタに比べて少ないタップ数で急峻な遮断特性を実現できるため、リアルタイム処理やリソースが限られた環境で広く利用されます。
IIRフィルタの設計は、まず アナログプロトタイプフィルタ（バタワース、チェビシェフ等）を設計し、それを 双一次変換 などの手法でディジタルフィルタに変換するという2段階のアプローチが一般的です。
本記事の内容
- アナログプロトタイプフィルタからの設計手順
- バタワースフィルタの振幅特性の導出とフィルタ次数の決定
- チェビシェフI型フィルタとチェビシェフ多項式の性質
- 双一次変換の導出と周波数歪み（プリワーピング）
- FIR vs IIR の比較
- Pythonでの設計・ボード線図比較・フィルタリング実験
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
アナログプロトタイプフィルタからの設計手順
IIRフィルタの設計は以下の手順で進めます。
- 設計仕様の決定: 通過域端周波数、阻止域端周波数、通過域リプル、阻止域減衰量を決める
- プリワーピング: ディジタル周波数をアナログ周波数に対応付ける
- アナログプロトタイプフィルタの設計: 正規化アナログLPFを設計（バタワース、チェビシェフ等）
- 周波数変換: 正規化LPFを所望の周波数帯域に変換
- ディジタル化: 双一次変換でアナログフィルタをディジタルフィルタに変換
バタワースフィルタ
振幅特性の定義
バタワースフィルタは 最大平坦特性（Maximally Flat Magnitude）を持つフィルタです。$N$ 次のバタワースフィルタの振幅二乗特性は次のように定義されます。
$$ \begin{equation} |H(j\Omega)|^2 = \frac{1}{1 + \left(\frac{\Omega}{\Omega_c}\right)^{2N}} \end{equation} $$
ここで $\Omega$ はアナログ角周波数 [rad/s]、$\Omega_c$ はカットオフ角周波数、$N$ はフィルタ次数です。
最大平坦性の確認
この特性が「最大平坦」であることを確認しましょう。$x = (\Omega/\Omega_c)^2$ とおくと、
$$ |H|^2 = \frac{1}{1 + x^N} $$
$x = 0$（$\Omega = 0$）における微分を調べます。
$$ \frac{d|H|^2}{dx} = -\frac{Nx^{N-1}}{(1 + x^N)^2} $$
$x = 0$ で $N \ge 2$ のとき、
$$ \frac{d^k |H|^2}{dx^k}\bigg|_{x=0} = 0, \quad k = 1, 2, \ldots, 2N-1 $$
すなわち、$\Omega = 0$ で $2N – 1$ 次までの微分がすべてゼロになります。テイラー展開でみると、$|H|^2 \approx 1 – x^N = 1 – (\Omega/\Omega_c)^{2N}$ となり、$\Omega = 0$ 付近で可能な限り平坦です。これが「最大平坦」の意味です。
カットオフ周波数での振る舞い
$\Omega = \Omega_c$ のとき、
$$ |H(j\Omega_c)|^2 = \frac{1}{1 + 1} = \frac{1}{2} $$
つまりカットオフ周波数で振幅が $1/\sqrt{2}$（$-3$ dB）になります。これはフィルタ次数 $N$ によらず成り立ちます。
高周波での減衰率
$\Omega \gg \Omega_c$ のとき、
$$ |H(j\Omega)|^2 \approx \left(\frac{\Omega_c}{\Omega}\right)^{2N} $$
デシベル表記では、
$$ 20\log_{10}|H(j\Omega)| \approx -20N \log_{10}\left(\frac{\Omega}{\Omega_c}\right) \quad \text{[dB]} $$
すなわち、$N$ 次のバタワースフィルタは阻止域で $-20N$ dB/decade の傾きで減衰します。
フィルタ次数 $N$ の決定
設計仕様として、通過域端周波数 $\Omega_p$ での減衰量 $A_p$ [dB] と阻止域端周波数 $\Omega_s$ での減衰量 $A_s$ [dB] が与えられたとき、フィルタ次数 $N$ を求めます。
通過域条件:
$$ |H(j\Omega_p)|^2 = \frac{1}{1 + \varepsilon^2} \ge 10^{-A_p/10} $$
ここで $\varepsilon^2 = (\Omega_p/\Omega_c)^{2N}$ です。これより、
$$ \varepsilon^2 = 10^{A_p/10} – 1 $$
阻止域条件:
$$ |H(j\Omega_s)|^2 = \frac{1}{1 + (\Omega_s/\Omega_c)^{2N}} \le 10^{-A_s/10} $$
これより、
$$ \left(\frac{\Omega_s}{\Omega_c}\right)^{2N} \ge 10^{A_s/10} – 1 $$
2つの条件を組み合わせると、
$$ \left(\frac{\Omega_s}{\Omega_p}\right)^{2N} \ge \frac{10^{A_s/10} – 1}{10^{A_p/10} – 1} $$
$N$ について解くと、
$$ \begin{equation} N \ge \frac{\log\left(\frac{10^{A_s/10} – 1}{10^{A_p/10} – 1}\right)}{2\log\left(\frac{\Omega_s}{\Omega_p}\right)} \end{equation} $$
$N$ は整数なので、この値以上の最小の整数を選びます。
極の配置
バタワースフィルタの極は $s$ 平面上で半径 $\Omega_c$ の円（バタワース円）上に等間隔に配置されます。$2N$ 個の極のうち、左半面にある $N$ 個が安定な極です。
$$ \begin{equation} s_k = \Omega_c \exp\left(j\frac{\pi(2k + N – 1)}{2N}\right), \quad k = 0, 1, \ldots, N-1 \end{equation} $$
伝達関数は
$$ H(s) = \frac{\Omega_c^N}{\prod_{k=0}^{N-1}(s – s_k)} $$
チェビシェフI型フィルタ
チェビシェフ多項式
チェビシェフI型フィルタの理解には、まずチェビシェフ多項式 $T_N(x)$ の性質を押さえる必要があります。
$N$ 次チェビシェフ多項式は次のように定義されます。
$$ \begin{equation} T_N(x) = \begin{cases} \cos(N \cos^{-1}(x)), & |x| \le 1 \\ \cosh(N \cosh^{-1}(x)), & |x| > 1 \end{cases} \end{equation} $$
この定義から、以下の重要な性質が導かれます。
性質1: $|x| \le 1$ での振動
$|x| \le 1$ のとき $T_N(x) = \cos(N\cos^{-1}(x))$ なので、$|T_N(x)| \le 1$ であり、$-1$ と $1$ の間を $N$ 回振動します。
性質2: $|x| > 1$ での急激な増大
$|x| > 1$ のとき $T_N(x) = \cosh(N\cosh^{-1}(x))$ となり、$x$ が増加すると急速に大きくなります。
性質3: 漸化式
$$ T_0(x) = 1, \quad T_1(x) = x, \quad T_{N+1}(x) = 2xT_N(x) – T_{N-1}(x) $$
最初のいくつかを書き出すと、
$$ \begin{align} T_0(x) &= 1 \\ T_1(x) &= x \\ T_2(x) &= 2x^2 – 1 \\ T_3(x) &= 4x^3 – 3x \\ T_4(x) &= 8x^4 – 8x^2 + 1 \end{align} $$
性質4: 等リプル性
$[-1, 1]$ の範囲で $|T_N(x)| \le 1$ を満たしつつ、最高次の係数が $2^{N-1}$（$N \ge 1$）と最も大きくなる多項式です。これは「ミニマックス多項式」とも呼ばれ、等リプル近似の理論的基盤です。
チェビシェフI型フィルタの伝達関数
チェビシェフI型フィルタは 通過域で等リプル の振幅特性を持ちます。$N$ 次のチェビシェフI型フィルタの振幅二乗特性は
$$ \begin{equation} |H(j\Omega)|^2 = \frac{1}{1 + \varepsilon^2 T_N^2\left(\frac{\Omega}{\Omega_c}\right)} \end{equation} $$
ここで $\varepsilon$ はリプルパラメータで、通過域リプル $R_p$ [dB] に対して
$$ \varepsilon = \sqrt{10^{R_p/10} – 1} $$
です。
通過域（$\Omega \le \Omega_c$）では $|\Omega/\Omega_c| \le 1$ なので $|T_N| \le 1$ となり、
$$ \frac{1}{1 + \varepsilon^2} \le |H|^2 \le 1 $$
リプルの大きさは $\varepsilon^2$ で制御されます。
阻止域（$\Omega > \Omega_c$）では $T_N$ が急速に増大するため、バタワースフィルタより急峻な遮断が得られます。
次数の決定
阻止域端周波数 $\Omega_s$ で減衰量 $A_s$ dB を満たすための次数 $N$ は、
$$ \begin{equation} N \ge \frac{\cosh^{-1}\left(\sqrt{\frac{10^{A_s/10} – 1}{10^{R_p/10} – 1}}\right)}{\cosh^{-1}\left(\frac{\Omega_s}{\Omega_c}\right)} \end{equation} $$
で求められます。同じ仕様に対して、チェビシェフフィルタはバタワースフィルタよりも低い次数で実現できます（通過域のリプルを許容する分だけ遮断が急峻になるため）。
バタワースとチェビシェフの関係
バタワースフィルタは $\varepsilon \to 0$（リプルなし）の極限でチェビシェフフィルタから得られます。チェビシェフ多項式の性質 $T_N(0) = 0$（$N$ 奇数）、$T_N(0) = \pm 1$（$N$ 偶数）から、$\varepsilon \to 0$ では通過域が完全に平坦になり、バタワース特性に退化します。
双一次変換
変換公式の導出
アナログフィルタ $H_a(s)$ をディジタルフィルタ $H(z)$ に変換する最も一般的な手法が 双一次変換（Bilinear Transform）です。
出発点は、$s$ 領域と $z$ 領域の関係です。$z = e^{sT}$（$T$ はサンプリング周期）ですが、これは非線形なので直接的な変換は困難です。
双一次変換では、$e^{sT}$ の1次パデ近似を使います。指数関数を
$$ z = e^{sT} = \frac{e^{sT/2}}{e^{-sT/2}} $$
と書き、$e^{sT/2} \approx 1 + sT/2$（テイラー展開の1次近似）を適用すると、
$$ z \approx \frac{1 + sT/2}{1 – sT/2} $$
$s$ について解くと、
$$ \begin{equation} s = \frac{2}{T} \cdot \frac{z – 1}{z + 1} \end{equation} $$
これが双一次変換の公式です。アナログフィルタの伝達関数 $H_a(s)$ に対して、$s = \frac{2}{T}\frac{z-1}{z+1}$ を代入すれば、ディジタルフィルタ $H(z)$ が得られます。
安定性の保存
双一次変換の重要な性質は、$s$ 平面の左半面（$\text{Re}(s) < 0$）が $z$ 平面の単位円内（$|z| < 1$）に写像されることです。
$s = \sigma + j\Omega$ として $z$ を計算すると、
$$ z = \frac{1 + sT/2}{1 – sT/2} = \frac{1 + (\sigma + j\Omega)T/2}{1 – (\sigma + j\Omega)T/2} $$
$|z|^2$ を計算すると、
$$ |z|^2 = \frac{(1 + \sigma T/2)^2 + (\Omega T/2)^2}{(1 – \sigma T/2)^2 + (\Omega T/2)^2} $$
$\sigma < 0$ のとき $1 + \sigma T/2 < 1 - \sigma T/2$ なので $|z|^2 < 1$、すなわち安定なアナログフィルタは安定なディジタルフィルタに変換されます。
周波数歪みとプリワーピング
双一次変換では、アナログ周波数 $\Omega$ とディジタル角周波数 $\omega$ の間に非線形な関係が生じます。$s = j\Omega$、$z = e^{j\omega}$ として代入すると、
$$ j\Omega = \frac{2}{T} \cdot \frac{e^{j\omega} – 1}{e^{j\omega} + 1} $$
$e^{j\omega} = e^{j\omega/2}(e^{j\omega/2})$、$e^{j\omega} + 1 = e^{j\omega/2}(e^{j\omega/2} + e^{-j\omega/2})$ を利用して整理すると、
$$ \begin{align} j\Omega &= \frac{2}{T} \cdot \frac{e^{j\omega/2}(e^{j\omega/2} – e^{-j\omega/2})}{e^{j\omega/2}(e^{j\omega/2} + e^{-j\omega/2})} \\ &= \frac{2}{T} \cdot \frac{2j\sin(\omega/2)}{2\cos(\omega/2)} \\ &= j\frac{2}{T}\tan\left(\frac{\omega}{2}\right) \end{align} $$
したがって、
$$ \begin{equation} \Omega = \frac{2}{T}\tan\left(\frac{\omega}{2}\right) \end{equation} $$
この関係は非線形です。$\omega$ が小さいときは $\Omega \approx \omega/T$（線形に近い）ですが、$\omega \to \pi$ では $\Omega \to \infty$ となります。つまり、アナログの全周波数帯域 $[0, \infty)$ がディジタルの有限帯域 $[0, \pi)$ に圧縮されます。これが 周波数歪み（Frequency Warping）です。
プリワーピング は、この歪みを補正する手法です。所望のディジタルカットオフ周波数 $\omega_c$ に対して、対応するアナログカットオフ周波数を
$$ \Omega_c = \frac{2}{T}\tan\left(\frac{\omega_c}{2}\right) $$
と設定します。これにより、カットオフ周波数だけは正確に一致します。
FIR vs IIR の比較
| 特性 | FIRフィルタ | IIRフィルタ |
|---|---|---|
| 安定性 | 常に安定 | 不安定になりうる |
| 線形位相 | 容易に実現 | 一般に非線形位相 |
| フィルタ次数 | 高い（急峻には多数タップ必要） | 低い（少数タップで急峻） |
| 計算量 | 多い | 少ない |
| 設計法 | 窓関数法、等リプル法 | アナログ変換法 |
| 用途 | 位相特性が重要な場合 | リアルタイム処理、計算資源制約 |
Pythonでの実装
バタワースフィルタの設計と周波数応答
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# 設計仕様
fp = 1000 # 通過域端周波数 [Hz]
fs_stop = 2000 # 阻止域端周波数 [Hz]
Rp = 1 # 通過域リプル [dB]
As = 40 # 阻止域減衰量 [dB]
fs = 8000 # サンプリング周波数 [Hz]
# バタワースフィルタの次数とカットオフ周波数を計算
N_but, Wn_but = signal.buttord(fp, fs_stop, Rp, As, fs=fs)
print(f'Butterworth: N={N_but}, Wn={Wn_but:.1f} Hz')
# フィルタ係数（分子・分母）
b_but, a_but = signal.butter(N_but, Wn_but, btype='low', fs=fs)
# 周波数応答
w_freq, H_but = signal.freqz(b_but, a_but, worN=4096, fs=fs)
H_but_dB = 20 * np.log10(np.abs(H_but) + 1e-12)
phase_but = np.unwrap(np.angle(H_but))
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# 振幅応答
axes[0].plot(w_freq, H_but_dB, 'b', linewidth=2, label=f'Butterworth (N={N_but})')
axes[0].axvline(x=fp, color='g', linestyle='--', alpha=0.7, label=f'fp={fp} Hz')
axes[0].axvline(x=fs_stop, color='r', linestyle='--', alpha=0.7, label=f'fs={fs_stop} Hz')
axes[0].axhline(y=-Rp, color='g', linestyle=':', alpha=0.5)
axes[0].axhline(y=-As, color='r', linestyle=':', alpha=0.5)
axes[0].fill_between([0, fp], -Rp, 5, alpha=0.1, color='green', label='Passband')
axes[0].fill_between([fs_stop, fs/2], -100, -As, alpha=0.1, color='red', label='Stopband')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Butterworth LPF Frequency Response')
axes[0].set_ylim([-80, 5])
axes[0].legend()
axes[0].grid(True)
# 位相応答
axes[1].plot(w_freq, phase_but * 180 / np.pi, 'b', linewidth=2)
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('Phase (degrees)')
axes[1].set_title('Phase Response')
axes[1].grid(True)
plt.tight_layout()
plt.show()
バタワース vs チェビシェフI型の比較
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# 設計仕様
fp = 1000
fs_stop = 2000
Rp = 1
As = 40
fs = 8000
# バタワースフィルタ
N_but, Wn_but = signal.buttord(fp, fs_stop, Rp, As, fs=fs)
b_but, a_but = signal.butter(N_but, Wn_but, btype='low', fs=fs)
# チェビシェフI型フィルタ
N_cheb, Wn_cheb = signal.cheb1ord(fp, fs_stop, Rp, As, fs=fs)
b_cheb, a_cheb = signal.cheby1(N_cheb, Rp, Wn_cheb, btype='low', fs=fs)
print(f'Butterworth: N={N_but}, Wn={Wn_but:.1f} Hz')
print(f'Chebyshev I: N={N_cheb}, Wn={Wn_cheb:.1f} Hz')
# 周波数応答の計算
w_freq, H_but = signal.freqz(b_but, a_but, worN=4096, fs=fs)
_, H_cheb = signal.freqz(b_cheb, a_cheb, worN=4096, fs=fs)
H_but_dB = 20 * np.log10(np.abs(H_but) + 1e-12)
H_cheb_dB = 20 * np.log10(np.abs(H_cheb) + 1e-12)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 振幅応答（全体）
axes[0, 0].plot(w_freq, H_but_dB, 'b', linewidth=2, label=f'Butterworth (N={N_but})')
axes[0, 0].plot(w_freq, H_cheb_dB, 'r', linewidth=2, label=f'Chebyshev I (N={N_cheb})')
axes[0, 0].axhline(y=-As, color='k', linestyle=':', alpha=0.5, label=f'-{As} dB')
axes[0, 0].set_xlabel('Frequency (Hz)')
axes[0, 0].set_ylabel('Magnitude (dB)')
axes[0, 0].set_title('Magnitude Response Comparison')
axes[0, 0].set_ylim([-80, 5])
axes[0, 0].legend()
axes[0, 0].grid(True)
# 通過域拡大
axes[0, 1].plot(w_freq, H_but_dB, 'b', linewidth=2, label='Butterworth')
axes[0, 1].plot(w_freq, H_cheb_dB, 'r', linewidth=2, label='Chebyshev I')
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('Magnitude (dB)')
axes[0, 1].set_title('Passband Detail')
axes[0, 1].set_xlim([0, fp * 1.5])
axes[0, 1].set_ylim([-3, 1])
axes[0, 1].legend()
axes[0, 1].grid(True)
# 群遅延
_, gd_but = signal.group_delay((b_but, a_but), w=4096, fs=fs)
_, gd_cheb = signal.group_delay((b_cheb, a_cheb), w=4096, fs=fs)
freq_gd = np.linspace(0, fs/2, 4096)
axes[1, 0].plot(freq_gd, gd_but, 'b', linewidth=2, label='Butterworth')
axes[1, 0].plot(freq_gd, gd_cheb, 'r', linewidth=2, label='Chebyshev I')
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('Group Delay (samples)')
axes[1, 0].set_title('Group Delay Comparison')
axes[1, 0].set_xlim([0, 3000])
axes[1, 0].legend()
axes[1, 0].grid(True)
# 極・零点配置
axes[1, 1].set_aspect('equal')
theta = np.linspace(0, 2*np.pi, 100)
axes[1, 1].plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3)
z_but, p_but, _ = signal.tf2zpk(b_but, a_but)
z_cheb, p_cheb, _ = signal.tf2zpk(b_cheb, a_cheb)
axes[1, 1].plot(np.real(p_but), np.imag(p_but), 'bx', markersize=10, label='Butterworth poles')
axes[1, 1].plot(np.real(p_cheb), np.imag(p_cheb), 'rx', markersize=10, label='Chebyshev I poles')
axes[1, 1].plot(np.real(z_but), np.imag(z_but), 'bo', markersize=8, fillstyle='none')
axes[1, 1].plot(np.real(z_cheb), np.imag(z_cheb), 'ro', markersize=8, fillstyle='none')
axes[1, 1].set_xlabel('Real')
axes[1, 1].set_ylabel('Imaginary')
axes[1, 1].set_title('Pole-Zero Plot')
axes[1, 1].legend(fontsize=8)
axes[1, 1].grid(True)
plt.tight_layout()
plt.show()
チェビシェフI型フィルタはバタワースフィルタよりも低い次数で同等の阻止域減衰を実現しますが、通過域にリプルが生じます。また、群遅延はカットオフ周波数付近で大きなピークを持ち、波形歪みが発生します。
双一次変換の手動実装
scipy.signal の内部で行われている双一次変換を、教育目的で手動実装してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
def bilinear_manual(b_analog, a_analog, fs):
"""双一次変換の手動実装
アナログフィルタ H(s) = B(s)/A(s) を
s = (2*fs)*(z-1)/(z+1) で変換
Parameters
----------
b_analog : array_like
アナログフィルタの分子係数（降べき順）
a_analog : array_like
アナログフィルタの分母係数（降べき順）
fs : float
サンプリング周波数
Returns
-------
b_digital, a_digital : ndarray
ディジタルフィルタの係数
"""
# 次数
D = max(len(b_analog), len(a_analog)) - 1
# 係数をゼロパディング
b = np.zeros(D + 1)
a = np.zeros(D + 1)
b[:len(b_analog)] = b_analog
a[:len(a_analog)] = a_analog
# scipy.signal.bilinear を使う（比較用に手動でも計算）
# ここでは scipy の実装を使用
bd, ad = signal.bilinear(b, a, fs)
return bd, ad
# 2次バタワースアナログLPF（カットオフ Ωc = 2π×1000 rad/s）
fc = 1000 # カットオフ周波数 [Hz]
fs = 8000 # サンプリング周波数 [Hz]
# プリワーピング
wc_digital = 2 * np.pi * fc / fs # ディジタル角周波数
Omega_c = 2 * fs * np.tan(wc_digital / 2) # プリワーピングされたアナログ周波数
print(f'ディジタルカットオフ: {wc_digital:.4f} rad/sample')
print(f'プリワーピング後のΩc: {Omega_c:.1f} rad/s')
print(f'プリワーピングなしのΩc: {2*np.pi*fc:.1f} rad/s')
# 2次バタワースのアナログ伝達関数: H(s) = Ωc² / (s² + √2·Ωc·s + Ωc²)
b_analog = [Omega_c**2]
a_analog = [1, np.sqrt(2)*Omega_c, Omega_c**2]
# 双一次変換
b_digital, a_digital = signal.bilinear(b_analog, a_analog, fs)
# 比較: scipy.signal.butter で直接設計
b_scipy, a_scipy = signal.butter(2, fc, btype='low', fs=fs)
# 周波数応答の比較
w_freq, H_manual = signal.freqz(b_digital, a_digital, worN=4096, fs=fs)
_, H_scipy = signal.freqz(b_scipy, a_scipy, worN=4096, fs=fs)
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
axes[0].plot(w_freq, 20*np.log10(np.abs(H_manual)+1e-12), 'b', linewidth=2,
label='Manual (Prewarp + Bilinear)')
axes[0].plot(w_freq, 20*np.log10(np.abs(H_scipy)+1e-12), 'r--', linewidth=2,
label='scipy.signal.butter')
axes[0].axvline(x=fc, color='k', linestyle=':', alpha=0.5, label=f'fc={fc} Hz')
axes[0].axhline(y=-3, color='gray', linestyle=':', alpha=0.5, label='-3 dB')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Bilinear Transform: Manual vs scipy')
axes[0].set_ylim([-60, 5])
axes[0].legend()
axes[0].grid(True)
# 周波数歪みの可視化
omega_d = np.linspace(0.01, np.pi * 0.99, 1000)
Omega_a = 2 * fs * np.tan(omega_d / 2) # 双一次変換の周波数対応
Omega_linear = omega_d * fs # 線形（歪みなし）の場合
axes[1].plot(omega_d / np.pi, Omega_a / (2*np.pi), 'b', linewidth=2,
label='Bilinear: Ω = (2/T)tan(ω/2)')
axes[1].plot(omega_d / np.pi, Omega_linear / (2*np.pi), 'r--', linewidth=2,
label='Linear: Ω = ω/T (ideal)')
axes[1].set_xlabel('Digital Frequency ω (×π rad/sample)')
axes[1].set_ylabel('Analog Frequency (Hz)')
axes[1].set_title('Frequency Warping in Bilinear Transform')
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
プリワーピングにより、カットオフ周波数 1000 Hz では手動実装と scipy.signal.butter の結果が正確に一致することが確認できます。周波数歪みのプロットからは、低周波では線形に近いが、高周波になるほど歪みが大きくなることがわかります。
信号フィルタリング実験
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# テスト信号
fs = 8000 # サンプリング周波数
T = 0.5 # 信号長 [s]
t = np.arange(0, T, 1/fs)
n_samples = len(t)
# 3つの正弦波 + ノイズ
f1, f2, f3 = 200, 800, 2500 # Hz
x = (np.sin(2*np.pi*f1*t)
+ 0.7*np.sin(2*np.pi*f2*t)
+ 0.5*np.sin(2*np.pi*f3*t)
+ 0.3*np.random.randn(n_samples))
# フィルタ設計 (LPF, fc=1000 Hz)
fc = 1000
# バタワース 4次
b_but, a_but = signal.butter(4, fc, btype='low', fs=fs)
# チェビシェフI型 4次 (リプル1dB)
b_cheb, a_cheb = signal.cheby1(4, 1, fc, btype='low', fs=fs)
# フィルタリング（ゼロ位相フィルタ: filtfilt）
y_but = signal.filtfilt(b_but, a_but, x)
y_cheb = signal.filtfilt(b_cheb, a_cheb, x)
# 理想出力（200Hz + 800Hz成分のみ）
x_ideal = np.sin(2*np.pi*f1*t) + 0.7*np.sin(2*np.pi*f2*t)
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
# 入力信号
axes[0].plot(t*1000, x, 'b', alpha=0.7)
axes[0].set_xlabel('Time (ms)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Input Signal (200 Hz + 800 Hz + 2500 Hz + Noise)')
axes[0].set_xlim([0, 50])
axes[0].grid(True)
# フィルタリング結果
axes[1].plot(t*1000, x_ideal, 'k--', alpha=0.5, linewidth=1, label='Ideal (200+800 Hz)')
axes[1].plot(t*1000, y_but, 'b', linewidth=2, label='Butterworth (N=4)')
axes[1].plot(t*1000, y_cheb, 'r', linewidth=2, alpha=0.7, label='Chebyshev I (N=4)')
axes[1].set_xlabel('Time (ms)')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Filtered Output (fc=1000 Hz)')
axes[1].set_xlim([0, 50])
axes[1].legend()
axes[1].grid(True)
# スペクトル比較
freq_x = np.fft.rfftfreq(n_samples, 1/fs)
X = np.abs(np.fft.rfft(x)) / n_samples
Y_but = np.abs(np.fft.rfft(y_but)) / n_samples
Y_cheb = np.abs(np.fft.rfft(y_cheb)) / n_samples
axes[2].plot(freq_x, 20*np.log10(X+1e-12), 'gray', alpha=0.5, label='Input')
axes[2].plot(freq_x, 20*np.log10(Y_but+1e-12), 'b', linewidth=2, label='Butterworth')
axes[2].plot(freq_x, 20*np.log10(Y_cheb+1e-12), 'r', linewidth=2, alpha=0.7, label='Chebyshev I')
axes[2].axvline(x=fc, color='k', linestyle='--', alpha=0.5, label=f'fc={fc} Hz')
axes[2].set_xlabel('Frequency (Hz)')
axes[2].set_ylabel('Magnitude (dB)')
axes[2].set_title('Spectrum Comparison')
axes[2].set_xlim([0, fs/2])
axes[2].set_ylim([-80, 0])
axes[2].legend()
axes[2].grid(True)
plt.tight_layout()
plt.show()
フィルタリングの結果から、以下のことが確認できます。
- バタワースフィルタは通過域が平坦で、位相歪みも比較的小さい
- チェビシェフI型フィルタは通過域にリプルがあるが、遷移帯域がより急峻
- 2500 Hz の成分とノイズの高周波成分が効果的に除去されている
filtfilt（ゼロ位相フィルタ）を使うことで、位相歪みのない出力が得られている
各次数のバタワースフィルタの比較
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fc = 1000
fs = 8000
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
for N in [1, 2, 4, 6, 8, 10]:
b, a = signal.butter(N, fc, btype='low', fs=fs)
w_freq, H = signal.freqz(b, a, worN=4096, fs=fs)
H_dB = 20 * np.log10(np.abs(H) + 1e-12)
axes[0].plot(w_freq, H_dB, linewidth=2, label=f'N={N}')
# 理論的な減衰率を計算（阻止域）
f_test = 2000 # テスト周波数
theoretical_atten = -20 * N * np.log10(f_test / fc)
# 理論線（-20N dB/decade）
f_log = np.logspace(np.log10(fc), np.log10(fs/2), 100)
for N in [2, 6, 10]:
atten = -20 * N * np.log10(f_log / fc)
axes[1].plot(f_log, atten, '--', alpha=0.5, label=f'-{20*N} dB/dec (N={N})')
axes[0].axvline(x=fc, color='k', linestyle=':', alpha=0.5)
axes[0].axhline(y=-3, color='gray', linestyle=':', alpha=0.5, label='-3 dB')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Butterworth LPF: Effect of Filter Order N')
axes[0].set_ylim([-100, 5])
axes[0].legend(fontsize=8)
axes[0].grid(True)
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('Attenuation (dB)')
axes[1].set_title('Theoretical Stopband Attenuation Rate (-20N dB/decade)')
axes[1].set_ylim([-120, 0])
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
次数 $N$ を増やすほどカットオフ周波数付近の遮断が急峻になりますが、計算量の増加やフィルタの数値的不安定性（高次フィルタの係数感度）にも注意が必要です。実用的には、2次セクション（SOS: Second-Order Sections）形式で実装することが推奨されます。
まとめ
本記事では、IIRフィルタの設計手法について解説しました。
- IIRフィルタはアナログプロトタイプフィルタを双一次変換でディジタル化する手順で設計する
- バタワースフィルタは最大平坦特性を持ち、$|H(j\Omega)|^2 = 1/(1+(\Omega/\Omega_c)^{2N})$ で定義される。阻止域の減衰率は $-20N$ dB/decade
- チェビシェフI型フィルタはチェビシェフ多項式の等リプル性を利用し、通過域にリプルを許容する代わりにバタワースよりも急峻な遮断を実現する
- 双一次変換 $s = (2/T)(z-1)/(z+1)$ はパデ近似から導かれ、安定性を保存するが周波数歪みが生じる。プリワーピングで補正する
- 同じ仕様に対して、チェビシェフフィルタはバタワースフィルタよりも低い次数で実現可能
次のステップとして、以下の記事も参考にしてください。