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

楕円フィルタ（Elliptic / Cauer Filter）の設計理論と Python 実装を scipy.signal.ellip・scipy.signal.ellipord・scipy.signal.freqz・scipy.signal.filtfilt で解説。ヤコビ楕円関数 sn/cn/dn と等リップル特性の原理、IIR ローパス/ハイパス/バンドパス/バンドストップ設計、buttord・cheb1ord・cheb2ord による Butterworth/Chebyshev I/II との最小次数・群遅延比較、filtfilt と lfilter の使い分けまで網羅。

19 min read

はじめに

楕円フィルタ（Elliptic Filter）は、通過域と阻止域の両方に等リップルを持つIIRフィルタです。発案者の名前からCouerフィルタとも呼ばれます。

代表的なIIRフィルタの遷移帯域の急峻さを比較すると次のようになります（同次数の場合）：

\[\text{楕円} > \text{チェビシェフ（I型・II型）} > \text{バターワース}\]

楕円フィルタは同じ次数のIIRフィルタの中で最も急峻な遷移帯域を実現できますが、通過域・阻止域の両方にリップルが生じます。この記事では数学的背景から、SciPyを使った実践的な設計・実装までを解説します。

この等リップル特性は、チェビシェフの弟子であったロシアの数学者エゴール・ゾロタレフ（Yegor Zolotarev, 1847–1878）が1877年に解いた「二つの区間上での最良有理近似問題」（ゾロタレフ問題）に由来します。ゾロタレフはベルリンでワイエルシュトラスから楕円関数論を学び、それを応用してこの問題を解決しました。この解を電気フィルタの設計に初めて応用したのがヴィルヘルム・カウアー（Wilhelm Cauer）で、1933年のことです。Cauerフィルタという別名はここに由来します。

楕円フィルタの周波数特性

楕円フィルタのゲイン特性は次の式で表されます：

\[|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\) ：フィルタ次数

ヤコビ楕円有理関数：なぜ両方の帯域で等リップルになるのか

チェビシェフフィルタ のI型は、チェビシェフ多項式 \(T_n(x) = \cos(n \arccos x)\) を使って通過域の等リップルを実現します。この構成の急所は置換 \(x = \cos\theta\) にあります。\(\theta\) が実数の範囲では \(\cos\) は周期 \(2\pi\) の単一の実周期しか持たないため、\(|x| \le 1\) （\(\theta\) が実数）では等振幅で振動しますが、\(|x| > 1\) では \(\theta\) が純虚数になり \(\cos\theta \to \cosh u\) となって単調に発散します。これがチェビシェフI型の阻止域が単調減少にしかならない理由です。II型は同じ構成を逆数変数 \(\Omega_s/\Omega\) に適用することで等リップル領域を阻止域側に移しますが、通過域は単調になってしまい、依然として「等リップルにできるのは片側だけ」という制約から逃れられません。

楕円フィルタは、この制約を同時に破ります。\(|x|\le 1\) でも \(|x|\ge \xi\) でも等リップルにするには、\(\cos\theta\) のような「単一の実周期しか持たない関数」ではなく、**実周期と虚周期の両方を持つ二重周期関数（楕円関数）**による置換が必要です。ヤコビの楕円関数 \(\text{sn}(u,k)\) は実周期 \(4K(k)\) と虚周期 \(2iK'(k)\) （\(K'(k)=K(\sqrt{1-k^2})\) は補完楕円積分）を持ち、実軸方向の周期性が通過域の等リップルを、虚軸方向の周期性が阻止域の等リップルをそれぞれ生み出します。この二重周期性を利用して \(n\) 次の有理関数

\[ R_n(\xi, x) = \text{cd}\!\left(n\, \frac{K(1/\xi)}{K(k)}\, \text{cd}^{-1}(x, k),\ k\right) \]

（\(\text{cd} = \text{cn}/\text{dn}\) はヤコビ楕円関数の一種、\(k\) は \(\xi\) から定まるモジュラス）を構成すると、\(R_n\) は次の性質を満たします：

  • \(|x| \le 1\) のとき \(|R_n(\xi, x)| \le 1\) 、かつ区間内で \(n\) 回等振幅に振動する（通過域の等リップル）
  • \(|x| \ge \xi\) のとき \(|R_n(\xi, x)| \ge \xi\) 、かつ区間内で等振幅に振動する（阻止域の等リップル）

この \(R_n(\xi,x)\) はゾロタレフ有理関数（Zolotarev rational function）と呼ばれ、チェビシェフ多項式の「多項式版・単一周期」の等リップル理論を「有理関数版・二重周期」に一般化したものです（Orfanidis, 2006）。式(1)の \(R_n\) にこれを代入すると、通過域・阻止域の両方が等リップルになる楕円フィルタの振幅特性が得られます。

次数の決定

要求仕様（通過域端 \(\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 で自動計算できます。

式(2)の導出：楕円ノームによる次数決定則

式(2)がなぜ成り立つのかを、楕円関数の**ノーム（nome）**という量を使って導きます。次のパラメータを定義します。

  • 選択度: \(k = \Omega_p / \Omega_s\) （\(0 < k < 1\) 、式(2)の \(1/\xi\) に相当）
  • 弁別度:
\[ k_1 = \frac{\varepsilon}{\varepsilon_1}, \qquad \varepsilon = \sqrt{10^{R_p/10}-1}, \qquad \varepsilon_1 = \sqrt{10^{R_s/10}-1} \]
  • 補モジュラス \(k' = \sqrt{1-k^2}\) 、\(K'(k) = K(k')\)
  • 楕円ノーム: \(q(k) = \exp(-\pi K'(k)/K(k))\)

前節のゾロタレフ有理関数 \(R_n(\xi,x)\) は、モジュラス \(k\) を持つ「入力側」の楕円関数配置を、\(n\) 重の合成によってモジュラス \(k_1\) を持つ「出力側」の配置へ写す構成になっています。この \(n\) 重合成は、三角関数における倍角公式の反復 \(\cos(n\theta) = T_n(\cos\theta)\) の楕円関数版（ランデン変換の \(n\) 回反復）に対応し、ノームに対しては単純な冪乗則

\[q(k_1) = \big[q(k)\big]^{\,n}\]

として作用します（Orfanidis, 2006）。この式を \(n\) について解くと、

\[ n = \frac{\ln q(k_1)}{\ln q(k)} = \frac{K(k)\, K'(k_1)}{K'(k)\, K(k_1)} \]

となり、右辺は \(\xi = 1/k\) 、\(\xi_s = 1/k_1\) と置き換えれば式(2)に一致します。等号は要求仕様をちょうど満たす（実数の）次数を与えるため、実際に使用する次数は \(n\) 以上の最小の整数（\(\lceil n \rceil\) ）です。

数値検証：\(f_p=100\) Hz, \(f_s=150\) Hz, \(R_p=1\) dB, \(R_s=60\) dB のとき、\(k = 100/150 = 0.6667\) 、\(\varepsilon = 0.5088\) 、\(\varepsilon_1 = 999.9995\) より \(k_1 = \varepsilon/\varepsilon_1 = 5.088 \times 10^{-4}\) です。実際に scipy.special.ellipk で \(K, K'\) を計算すると

import numpy as np
from scipy.special import ellipk
from scipy import signal

rp, rs = 1.0, 60.0
fp, fs_stop = 100, 150

eps  = np.sqrt(10**(rp/10) - 1)
eps1 = np.sqrt(10**(rs/10) - 1)
k  = fp / fs_stop
k1 = eps / eps1

def K(m):
    return ellipk(m**2)  # scipy.special.ellipk は m = k^2 を引数に取る

Kk, Kpk   = K(k),  K(np.sqrt(1 - k**2))
Kk1, Kpk1 = K(k1), K(np.sqrt(1 - k1**2))

n_formula = (Kk * Kpk1) / (Kpk * Kk1)
print(f"k={k:.4f}  k1={k1:.6e}")
print(f"n (実数解) = {n_formula:.4f}  -> 採用する次数 = {int(np.ceil(n_formula))}")

n_ellipord, _ = signal.ellipord(fp/500, fs_stop/500, rp, rs)
print(f"scipy.signal.ellipord による次数 = {n_ellipord}")

を実行すると、

k=0.6667  k1=5.088473e-04
n (実数解) = 5.4267  -> 採用する次数 = 6
scipy.signal.ellipord による次数 = 6

となり、式(2)（＝ノームの冪乗則から導いた式）の理論値 \(n=5.4267\) を切り上げた \(6\) 次と、scipy.signal.ellipord が実際に返す次数 \(6\) 次が完全に一致することが確認できます。

極・零点配置：有限伝達零点

急峻な遷移帯域は、伝達関数 \(H(s)\) の零点配置によって実現されます。ゾロタレフ有理関数 \(R_n(\xi,x)\) は \(x=\xi_1, \xi_2, \dots\) （\(\xi\) 以上の有限の点）で \(0\) になるため、対応する \(H(s)\) は虚軸上の有限の点 \(s=j\Omega_1, j\Omega_2,\dots\) に零点を持ちます。これを**有限伝達零点（finite transmission zeros）**と呼びます。有限伝達零点の位置ちょうどでゲインが理論上 \(-\infty\) dB になり、この「深い谷」が阻止域の等リップルの谷を作っています。

ここで、Butterworth・チェビシェフI型・チェビシェフII型・楕円フィルタの零点配置は同じではないという点に注意が必要です。実際に scipy.signal のアナログプロトタイプ関数で確認します。

from scipy import signal

N = 6
z_butter, _, _ = signal.buttap(N)
z_cheby1, _, _ = signal.cheb1ap(N, 1.0)
z_cheby2, _, _ = signal.cheb2ap(N, 60.0)
z_ellip,  _, _ = signal.ellipap(N, 1.0, 60.0)

print("バターワース  の有限零点数:", len(z_butter))
print("チェビシェフI の有限零点数:", len(z_cheby1))
print("チェビシェフII の有限零点数:", len(z_cheby2))
print("楕円          の有限零点数:", len(z_ellip))

実行結果：

バターワース  の有限零点数: 0
チェビシェフI の有限零点数: 0
チェビシェフII の有限零点数: 6
楕円          の有限零点数: 6

バターワースとチェビシェフI型は多項式（\(T_n\) そのもの）で振幅特性を構成するため有限零点を持たず、全ての零点は \(s=\infty\) にあります（阻止域は単調にしか減衰しない一因です）。一方チェビシェフII型は、I型の構成を逆数変数に適用した有理関数であるため、実は楕円フィルタと同様に阻止域に有限零点を持ちます。つまり「有限伝達零点を持つかどうか」の境界は「Butterworth/チェビシェフ vs 楕円」ではなく、「多項式近似（Butterworth・チェビシェフI型）vs 有理関数近似（チェビシェフII型・楕円）」にあります。楕円フィルタが特別なのは、この有理関数近似を通過域側にも同時に適用し、両側を等リップルにしている点です。

有限伝達零点の個数にも次数の偶奇でエッジケースがあります。\(n\) 次楕円フィルタの有限零点対の数は、\(n\) が偶数なら \(n/2\) 対（全ての極零点ペアが有限）、\(n\) が奇数なら \((n-1)/2\) 対（1つの極が対応する零点を持たず \(s=\infty\) 側に「余る」）になります。

from scipy import signal
for N in [5, 6, 7]:
    z, p, k = signal.ellipap(N, 1.0, 60.0)
    print(f"n={N}: 極数={len(p)}  有限零点数={len(z)}")

n=5: 極数=5  有限零点数=4
n=6: 極数=6  有限零点数=6
n=7: 極数=7  有限零点数=6

この非対称性は実用上重要な帰結を持ちます。偶数次では零点対の数と極対の数が一致するため、\(\Omega\to\infty\) でゲインは阻止域リップルの床（\(-R_s\) dB）に漸近し、それ以上下がりません。奇数次では極が1つ余るため、最後の有限零点を超えた先でゲインはさらに単調に減衰します。

実際に \(n=6\) （scipy.signal.ellipord(100, 150, 1, 60, fs=1000) が返す次数）で有限伝達零点の周波数を計算し、可視化します。

import numpy as np
from scipy import signal

fs = 1000
rp, rs = 1.0, 60.0
n, wn = signal.ellipord(100, 150, rp, rs, fs=fs)  # Hz単位で統一（fs=を渡すwnもHz）

z, p, k = signal.ellip(n, rp, rs, wn, btype='low', fs=fs, output='zpk')
b, a = signal.zpk2tf(z, p, k)

print("零点の絶対値（単位円上か確認）:", np.round(np.abs(z), 6))

# 零点自身の角周波数でH(z)を厳密に評価し、深さを確認
uniq_ang = sorted(set(round(abs(np.angle(zz)), 10) for zz in z))
_, h_exact = signal.freqz(b, a, worN=np.array(uniq_ang))
for wa, hh in zip(uniq_ang, h_exact):
    f_hz = wa * fs / (2 * np.pi)
    print(f"f={f_hz:6.2f} Hz   |H|={abs(hh):.2e}   ({20*np.log10(abs(hh)+1e-300):.1f} dB)")

実行結果：

零点の絶対値（単位円上か確認）: [1. 1. 1. 1. 1. 1.]
f=133.73 Hz   |H|=1.31e-12   (-237.6 dB)
f=163.41 Hz   |H|=3.02e-13   (-250.4 dB)
f=303.88 Hz   |H|=1.85e-14   (-274.6 dB)

全ての零点が単位円上（\(|z|=1\) 、つまり解析的には虚軸上）にあり、通過域端（100 Hz）と阻止域端（150 Hz）のすぐ外側から阻止域全体にかけて3対の有限零点が配置されていることが分かります。理論上のゲインは \(-\infty\) dB ですが、浮動小数点演算の丸め誤差により実測では \(-238\sim-275\) dB という、通常の用途では実質無限大とみなせる深さになっています。興味深いのは最初の零点（133.73 Hz）が、指定した阻止域端 150 Hz より内側（遷移帯域の途中）に位置することです。ellipord が保証するのは「阻止域端 150 Hz 以遠で \(R_s\) dB 以上の減衰」であって、それより手前で減衰が一時的にさらに深くなることは妨げられないためです。

次の図は、この楕円フィルタ（\(n=6\) ）の全体特性（上）と通過域の等リップル構造の拡大（下）を示したものです。上図の赤破線が有限伝達零点の位置で、いずれも阻止域リップルの谷の底と一致しています。

楕円ローパスフィルタの有限伝達零点と等リップル構造

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

フィルタ通過域阻止域遷移帯域（同次数）位相特性
バターワース最大平坦単調減少最もなだらか比較的良好
チェビシェフI型等リップル単調減少中程度やや悪化
チェビシェフII型最大平坦等リップル中程度やや悪化
楕円等リップル等リップル最も急峻最も悪化

楕円フィルタは遷移帯域の急峻さを最大化する代わりに、位相特性（群遅延の平坦性）が最も劣化します。位相特性を重視する場合はバターワースが、急峻な遮断が必要な場合は楕円フィルタが適しています。

群遅延特性とパルス応答のトレードオフ

急峻な振幅特性の代償を定量的に確認します。群遅延は位相特性の周波数微分の符号を反転したもの

\[\tau(\omega) = -\frac{d\phi(\omega)}{d\omega}\]

で定義され、群遅延が周波数によらず一定であれば、フィルタは線形位相特性を持ち、パルス波形の形状が保たれます。 ベッセルフィルタ は通過域でこの群遅延を最大限平坦にする設計思想であり、楕円フィルタとは正反対の設計指向を持ちます。楕円フィルタは通過域端付近に有理関数由来の急激な位相回転を持ち込むため、そこで群遅延が大きく変動します。

同一次数（\(n=4\) ）のバターワース・チェビシェフI型・楕円・ベッセルの4フィルタについて、通過域内（1〜90 Hz）の群遅延の変動幅と、単位ステップ応答のオーバーシュート率を実測します。

import numpy as np
from scipy import signal

fs = 1000
fc = 100
N = 4
rp, rs = 1.0, 60.0

filters = {
    'Butterworth': signal.butter(N, fc, btype='low', fs=fs),
    'Chebyshev I': signal.cheby1(N, rp, fc, btype='low', fs=fs),
    'Elliptic':    signal.ellip(N, rp, rs, fc, btype='low', fs=fs),
    'Bessel':      signal.bessel(N, fc, btype='low', fs=fs, norm='delay'),
}

print("=== 通過域(1-90Hz)における群遅延の変動 ===")
for name, ba in filters.items():
    w, gd_samples = signal.group_delay(ba, w=4096, fs=fs)
    gd_ms = gd_samples / fs * 1000.0          # サンプル数 -> ミリ秒
    mask = (w >= 1) & (w <= 90)
    gd = gd_ms[mask]
    print(f"{name:12s}: 平均={gd.mean():.3f}ms  "
          f"最小={gd.min():.3f}ms  最大={gd.max():.3f}ms  "
          f"変動幅={(gd.max()-gd.min()):.3f}ms")

print()
print("=== 単位ステップ応答のオーバーシュート ===")
for name, ba in filters.items():
    b, a = ba
    t = np.arange(0, 0.3, 1/fs)
    y = signal.lfilter(b, a, np.ones(len(t)))
    steady = y[-50:].mean()
    overshoot = (y.max() - steady) / steady * 100
    print(f"{name:12s}: 定常値={steady:.4f}  ピーク={y.max():.4f}  オーバーシュート={overshoot:.2f}%")

実行結果：

=== 通過域(1-90Hz)における群遅延の変動 ===
Butterworth : 平均=4.848ms  最小=4.021ms  最大=6.522ms  変動幅=2.501ms
Chebyshev I : 平均=5.697ms  最小=4.148ms  最大=10.450ms  変動幅=6.302ms
Elliptic    : 平均=5.539ms  最小=3.903ms  最大=10.598ms  変動幅=6.694ms
Bessel      : 平均=1.582ms  最小=1.539ms  最大=1.669ms  変動幅=0.130ms

=== 単位ステップ応答のオーバーシュート ===
Butterworth : 定常値=1.0000  ピーク=1.1191  オーバーシュート=11.91%
Chebyshev I : 定常値=0.8913  ピーク=1.0923  オーバーシュート=22.56%
Elliptic    : 定常値=0.8913  ピーク=1.1017  オーバーシュート=23.61%
Bessel      : 定常値=1.0000  ピーク=1.0807  オーバーシュート=8.07%

同じ4次でも、通過域内の群遅延変動幅はベッセルの0.130 msに対し楕円は6.694 msと、実に50倍以上に達します。これは表中で「位相特性：最も悪化」としていた定性的な評価を、実測値で裏付けるものです。ステップ応答のオーバーシュートもベッセルの8.07%に対し楕円は23.61%と最大で、急峻な振幅特性と引き換えに矩形的な入力に対して強いリンギング（振動的なオーバーシュート）が生じることが分かります。なお、チェビシェフI型（通過域のみ等リップル）と楕円（両側等リップル）のオーバーシュート・群遅延変動幅はこの条件では近い値になっており、通過域の非平坦性それ自体が線形位相からのずれの主要因であることも読み取れます。パルス波形の保持が最優先される用途（測定器・パルス伝送など）では楕円フィルタは避け、ベッセルフィルタを選ぶべきです。

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}次")

実際に実行した出力（SciPy 1.18.0）：

バターワース  : 17次
チェビシェフI :  9次
チェビシェフII:  9次
楕円          :  6次

楕円フィルタはバターワースの1/3以下の次数で同じ仕様を達成できます。次数は伝達関数の分母多項式の次数、すなわち必要な演算量・状態変数（メモリ）の数にほぼ直結するため、組み込み機器などリソースが限られる環境では次数の差がそのまま実装コストの差になります。この6次という値は、前節でノームの冪乗則から導いた式(2)の理論値 \(n=5.4267\) を切り上げた値と完全に一致しています。

ハイパス・バンドパスフィルタへの拡張

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

近年の研究動向

楕円フィルタに代表される古典的アナログプロトタイプ法（Butterworth・Chebyshev・楕円・ベッセル）は、通過域リップル・阻止域減衰・遷移帯域幅という3つの仕様を閉形式の式で満たす設計として、現在も信号処理・音響機器・通信システムの標準的な選択肢であり続けています。SciPyやMATLABが提供するのもこの古典的設計です。

一方で2023年以降、IIRフィルタの基本構造（特に2次セクション＝biquad）を深層学習パイプラインに組み込む研究が進んでいます。Lutati, Zimerman, Wolf (EMNLP 2023) の “Focus Your Attention (with Adaptive IIR Filters)” は、入力に応じて係数が変化する2次の適応IIRフィルタをTransformerの前処理層として使い、長系列タスクでAttentionの効率を高める手法を提案しています。また Malek, Schulz, Wuebbelmann (DAFx25, 2025) は、Kolmogorov-Arnold Networks（KAN）を使ってbiquadフィルタの係数そのものを学習によって最適化する手法を報告しています。

これらの研究は、本記事で扱ったような「与えられた通過域・阻止域仕様から最小次数の閉形式解を導く」古典的設計を置き換えるものではなく、IIRフィルタという構造そのものを、固定仕様に縛られない学習可能なビルディングブロックとしてニューラルネットワークに組み込む方向性である点に注意が必要です。古典的アナログプロトタイプ法とニューラルベースの手法は、現時点では代替関係というより、適用場面（固定された周波数仕様 vs 入力依存・タスク依存の適応的な特性）が異なる補完的な技術として併存していると捉えるのが実情に近いといえます。

関連記事

参考文献

  • Zolotarev, E. I. (1877). Application of elliptic functions to questions of functions deviating least and most from zero. Zapiski Imperatorskoi Akademii Nauk, 30(5), 1-71.
  • Cauer, W. (1933). Ein Interpolationsproblem mit Funktionen mit positivem Realteil. Mathematische Zeitschrift, 38, 1-44.
  • Orfanidis, S. J. (2006). Lecture Notes on Elliptic Filter Design. Rutgers University, Department of Electrical & Computer Engineering.
  • Proakis, J. G., & Manolakis, D. G. (2007). Digital Signal Processing (4th ed.). Pearson.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • Lutati, S., Zimerman, I., & Wolf, L. (2023). Focus Your Attention (with Adaptive IIR Filters). Proceedings of EMNLP 2023.
  • Malek, A., Schulz, D., & Wuebbelmann, F. (2025). Biquad Coefficients Optimization via Kolmogorov-Arnold Networks. Proceedings of the 28th International Conference on Digital Audio Effects (DAFx25).
  • SciPy scipy.signal.ellip documentation
  • SciPy scipy.signal.ellipord documentation

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox Frequency Converter Hz/kHz/MHz/rad/s 単位を相互変換 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション

関連記事