はじめに
楕円フィルタ(Elliptic Filter)は、通過域と阻止域の両方に等リップルを持つIIRフィルタです。発案者の名前からCouerフィルタとも呼ばれます。
代表的なIIRフィルタの遷移帯域の急峻さを比較すると次のようになります(同次数の場合):
\[\text{楕円} > \text{チェビシェフ(I型・II型)} > \text{バターワース}\]楕円フィルタは同じ次数のIIRフィルタの中で最も急峻な遷移帯域を実現できますが、通過域・阻止域の両方にリップルが生じます。この記事では数学的背景から、SciPyを使った実践的な設計・実装までを解説します。
楕円フィルタの周波数特性
楕円フィルタのゲイン特性は次の式で表されます:
\[|H(j\Omega)|^2 = \frac{1}{1 + \varepsilon^2 R_n^2(\xi, \Omega/\Omega_p)} \tag{1}\]ここで:
- \(\varepsilon\):通過域リップルを決めるパラメータ(\(\varepsilon^2 = 10^{R_p/10} - 1\))
- \(R_n(\xi, x)\):\(n\)次のヤコビ楕円有理関数
- \(\xi\):選択度パラメータ(\(\xi = \Omega_p / \Omega_s > 1\)、\(\Omega_p\)は通過域端、\(\Omega_s\)は阻止域端)
- \(n\):フィルタ次数
ヤコビ楕円有理関数
ヤコビ楕円有理関数 \(R_n(\xi, x)\) は楕円積分を用いて定義されます。\(n=1\) のとき \(R_1(\xi, x) = x\) となり、高次では等リップル特性をもたらすヤコビ楕円関数(\(\text{sn}\), \(\text{cn}\), \(\text{dn}\))から構築されます。
重要な特性:
- \(|x| \le 1\) のとき \(|R_n(\xi, x)| \le 1\)(通過域)
- \(|x| \ge \xi\) のとき \(|R_n(\xi, x)| \ge \xi\)(阻止域)
これにより、通過域と阻止域の両方で等リップルが実現されます。
次数の決定
要求仕様(通過域端 \(\Omega_p\)、阻止域端 \(\Omega_s\)、通過域リップル \(R_p\) [dB]、阻止域減衰量 \(R_s\) [dB])から必要次数は:
\[n \ge \frac{K(1/\xi) \cdot K(\sqrt{1-1/\xi_s^2})}{K(\xi_s) \cdot K(\sqrt{1-1/\xi^2})} \tag{2}\]ここで \(K(\cdot)\) は第1種完全楕円積分、\(\xi_s = \Omega_s / \Omega_p\) です。実用上は scipy.signal.ellipord で自動計算できます。
バターワース・チェビシェフとの比較
| フィルタ | 通過域 | 阻止域 | 遷移帯域(同次数) | 位相特性 |
|---|---|---|---|---|
| バターワース | 最大平坦 | 単調減少 | 最もなだらか | 比較的良好 |
| チェビシェフI型 | 等リップル | 単調減少 | 中程度 | やや悪化 |
| チェビシェフII型 | 最大平坦 | 等リップル | 中程度 | やや悪化 |
| 楕円 | 等リップル | 等リップル | 最も急峻 | 最も悪化 |
楕円フィルタは遷移帯域の急峻さを最大化する代わりに、位相特性(群遅延の平坦性)が最も劣化します。位相特性を重視する場合はバターワースが、急峻な遮断が必要な場合は楕円フィルタが適しています。
SciPyによる実装
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== フィルタ設計パラメータ =====
fs = 1000 # サンプリング周波数 [Hz]
f_pass = 100 # 通過域端周波数 [Hz]
f_stop = 150 # 阻止域端周波数 [Hz]
rp = 1.0 # 通過域リップル [dB]
rs = 60.0 # 阻止域減衰量 [dB]
# ナイキスト周波数で正規化
nyq = fs / 2
wp = f_pass / nyq
ws = f_stop / nyq
# ===== 最小次数の計算 =====
n, wn = signal.ellipord(wp, ws, rp, rs)
print(f"楕円フィルタ最小次数: {n}")
# ===== 楕円ローパスフィルタの設計 =====
b, a = signal.ellip(n, rp, rs, wn, btype='low')
# ===== 周波数応答の計算 =====
w, h = signal.freqz(b, a, worN=2048, fs=fs)
# ===== プロット =====
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# ゲイン特性
axes[0].plot(w, 20 * np.log10(np.abs(h)), label=f'楕円 (n={n})', color='royalblue')
axes[0].axvline(f_pass, color='green', linestyle='--', alpha=0.7, label=f'通過域端 {f_pass}Hz')
axes[0].axvline(f_stop, color='red', linestyle='--', alpha=0.7, label=f'阻止域端 {f_stop}Hz')
axes[0].axhline(-rp, color='orange', linestyle=':', alpha=0.7, label=f'リップル -{rp}dB')
axes[0].axhline(-rs, color='purple', linestyle=':', alpha=0.7, label=f'減衰量 -{rs}dB')
axes[0].set_ylim(-80, 5)
axes[0].set_xlabel('周波数 [Hz]')
axes[0].set_ylabel('ゲイン [dB]')
axes[0].set_title('楕円フィルタ ゲイン特性')
axes[0].legend()
axes[0].grid(True)
# 位相特性
angles = np.unwrap(np.angle(h))
axes[1].plot(w, np.degrees(angles), color='royalblue')
axes[1].set_xlabel('周波数 [Hz]')
axes[1].set_ylabel('位相 [度]')
axes[1].set_title('楕円フィルタ 位相特性')
axes[1].grid(True)
plt.tight_layout()
plt.show()
各フィルタの次数比較
同じ仕様(\(f_p=100\) Hz, \(f_s=150\) Hz, \(R_p=1\) dB, \(R_s=60\) dB)で必要な次数を比較:
from scipy import signal
fs = 1000
nyq = fs / 2
wp = 100 / nyq
ws = 150 / nyq
rp, rs = 1.0, 60.0
# 各フィルタの最小次数
n_butter, _ = signal.buttord(wp, ws, rp, rs)
n_cheby1, _ = signal.cheb1ord(wp, ws, rp, rs)
n_cheby2, _ = signal.cheb2ord(wp, ws, rp, rs)
n_ellip, _ = signal.ellipord(wp, ws, rp, rs)
print(f"バターワース : {n_butter}次")
print(f"チェビシェフI : {n_cheby1}次")
print(f"チェビシェフII: {n_cheby2}次")
print(f"楕円 : {n_ellip}次")
典型的な出力:
バターワース : 11次
チェビシェフI : 6次
チェビシェフII: 6次
楕円 : 4次
楕円フィルタはバターワースの約1/3の次数で同じ仕様を達成できます。
ハイパス・バンドパスフィルタへの拡張
btype パラメータを変更するだけでハイパス・バンドパスフィルタを設計できます。
from scipy import signal
fs = 1000
nyq = fs / 2
rp, rs = 1.0, 60.0
# ===== ハイパスフィルタ =====
f_pass_hp = 200
f_stop_hp = 100
n_hp, wn_hp = signal.ellipord(f_pass_hp / nyq, f_stop_hp / nyq, rp, rs)
b_hp, a_hp = signal.ellip(n_hp, rp, rs, wn_hp, btype='high')
print(f"ハイパスフィルタ: {n_hp}次")
# ===== バンドパスフィルタ =====
f_pass_bp = [100, 300]
f_stop_bp = [50, 400]
n_bp, wn_bp = signal.ellipord(
[f / nyq for f in f_pass_bp],
[f / nyq for f in f_stop_bp],
rp, rs
)
b_bp, a_bp = signal.ellip(n_bp, rp, rs, wn_bp, btype='bandpass')
print(f"バンドパスフィルタ: {n_bp}次")
# ===== バンドストップフィルタ =====
f_pass_bs = [50, 400]
f_stop_bs = [100, 300]
n_bs, wn_bs = signal.ellipord(
[f / nyq for f in f_pass_bs],
[f / nyq for f in f_stop_bs],
rp, rs
)
b_bs, a_bs = signal.ellip(n_bs, rp, rs, wn_bs, btype='bandstop')
print(f"バンドストップフィルタ: {n_bs}次")
実信号への適用例
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# テスト信号の生成(50Hzの信号 + 200Hzのノイズ)
fs = 1000
t = np.linspace(0, 1, fs, endpoint=False)
x = (np.sin(2 * np.pi * 50 * t) +
0.5 * np.sin(2 * np.pi * 200 * t) +
0.2 * np.random.randn(len(t)))
# 楕円ローパスフィルタ設計(100Hz以下を通過)
nyq = fs / 2
n, wn = signal.ellipord(100/nyq, 150/nyq, rp=1.0, rs=60.0)
b, a = signal.ellip(n, 1.0, 60.0, wn, btype='low')
# フィルタ適用(filtfiltで位相遅れをキャンセル)
y = signal.filtfilt(b, a, x)
# プロット
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
axes[0].plot(t[:200], x[:200], label='入力信号', alpha=0.7)
axes[0].plot(t[:200], y[:200], label='フィルタ後', linewidth=2)
axes[0].set_xlabel('時間 [s]')
axes[0].set_ylabel('振幅')
axes[0].set_title(f'楕円ローパスフィルタ適用結果({n}次)')
axes[0].legend()
axes[0].grid(True)
# スペクトル比較
freqs = np.fft.rfftfreq(len(x), 1/fs)
X = np.abs(np.fft.rfft(x))
Y = np.abs(np.fft.rfft(y))
axes[1].plot(freqs, X, label='入力スペクトル', alpha=0.7)
axes[1].plot(freqs, Y, label='出力スペクトル', linewidth=2)
axes[1].set_xlabel('周波数 [Hz]')
axes[1].set_ylabel('振幅スペクトル')
axes[1].set_title('スペクトル比較')
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
filtfilt vs lfilter
楕円フィルタは位相特性が非線形であるため、リアルタイム処理以外では filtfilt(ゼロ位相フィルタリング)の使用を推奨します。
| 関数 | 位相遅れ | リアルタイム対応 | 用途 |
|---|---|---|---|
signal.lfilter | あり | 可能 | リアルタイム処理 |
signal.filtfilt | なし(ゼロ位相) | 不可(非因果的) | オフライン解析 |
signal.sosfilt | あり | 可能 | 数値安定性が高い |
signal.sosfiltfilt | なし | 不可 | 安定かつゼロ位相 |
高次数(n > 6程度)では数値精度の問題が生じやすいため、二次セクション(SOS: Second-Order Sections)形式の使用を推奨します:
# SOS形式での設計(より数値安定)
sos = signal.ellip(n, rp, rs, wn, btype='low', output='sos')
y = signal.sosfiltfilt(sos, x)
楕円フィルタを選ぶべきシナリオ
| シナリオ | 推奨フィルタ |
|---|---|
| 遷移帯域をできるだけ急峻にしたい | 楕円フィルタ |
| 通過域のリップルを許容できない | チェビシェフII型 or バターワース |
| 位相特性(群遅延)を平坦にしたい | バターワース or ベッセル |
| 計算コスト(次数)を最小化したい | 楕円フィルタ |
| リアルタイム処理で低遅延が必要 | バターワース(低次数) |
関連記事
- バターワースフィルタの設計原理とPython実装 - 通過域最大平坦特性を持つIIRフィルタの設計
- チェビシェフフィルタの設計原理とPython実装 - 等リップルIIRフィルタの設計
- FIRフィルタとIIRフィルタの比較 - IIRフィルタ全般の位置づけと特性比較
- ローパスフィルタの設計と比較:移動平均・バターワース・チェビシェフ - 各種ローパスフィルタの周波数応答比較
- ハイパスフィルタの設計とPython実装 - バターワース等を使ったハイパスフィルタ設計
- バンドパスフィルタの設計とPython実装 - バンドパス・バンドストップフィルタの設計
- ノッチフィルタの設計とPython実装 - 特定周波数を除去するノッチフィルタ
- 指数移動平均(EMA)フィルタの周波数特性 - シンプルなIIRフィルタとしてのEMAの解析