楕円フィルタ(Elliptic Filter)の設計原理とPython実装

楕円フィルタ(Cauer Filter)の数学的基礎から、等リップル特性の導出、バターワース・チェビシェフとの比較、SciPyを使ったローパス/ハイパス/バンドパスフィルタのPython実装までを解説します。

はじめに

楕円フィルタ(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 ベッセル
計算コスト(次数)を最小化したい楕円フィルタ
リアルタイム処理で低遅延が必要バターワース(低次数)

関連記事