はじめに
ベッセルフィルタ（Bessel Filter）は、通過域で群遅延が最大限に平坦なIIRフィルタです。W. E. Thomsonが1952年に提案したことからThomson-Besselフィルタとも呼ばれます。
代表的なIIRフィルタの群遅延の平坦性を比較すると次のようになります（同次数の場合）：
\[\text{ベッセル} > \text{バターワース} > \text{チェビシェフ} > \text{楕円}\]ベッセルフィルタは通過帯域内で群遅延がほぼ一定であるため、異なる周波数成分の信号が同じ時間遅延でフィルタを通過します。これによりパルス波形の形状を保持できます。本記事では数学的背景から、SciPyを使った実践的な設計・実装までを解説します。
群遅延とは
群遅延（Group Delay）とは、フィルタの位相特性 \(\phi(\omega)\) の周波数に対する微分の負値です：
\[\tau(\omega) = -\frac{d\phi(\omega)}{d\omega} \tag{1}\]群遅延が一定（周波数によらず等しい）であることは、フィルタが線形位相特性を持つことを意味します。線形位相フィルタでは、すべての周波数成分が同じ時間だけ遅延するため、信号の波形がそのまま保たれます。
音声・音楽処理、パルス波形の伝送、測定機器など、波形の形状保持が重要な用途でベッセルフィルタが選ばれます。
ベッセルフィルタの数学的背景
ベッセル多項式
ベッセルフィルタの伝達関数はベッセル多項式 \(B_n(s)\) に基づいています。逆数正規化されたベッセル多項式（\(y_n(s)\) とも書く）は漸化式で定義されます：
\[B_0(s) = 1 \tag{2}\]\[B_1(s) = s + 1 \tag{3}\]\[B_n(s) = (2n-1)B_{n-1}(s) + s^2 B_{n-2}(s) \tag{4}\]低次のベッセル多項式を具体的に示します：
| 次数 \(n\) | \(B_n(s)\) |
|---|---|
| 1 | \(s + 1\) |
| 2 | \(s^2 + 3s + 3\) |
| 3 | \(s^3 + 6s^2 + 15s + 15\) |
| 4 | \(s^4 + 10s^3 + 45s^2 + 105s + 105\) |
伝達関数
\(n\) 次ベッセルローパスフィルタの伝達関数は次の形をとります：
\[H_n(s) = \frac{B_n(0)}{B_n(s/\omega_0)} = \frac{(2n-1)!!}{B_n(s/\omega_0)} \tag{5}\]ここで \(\omega_0\) は正規化角周波数、\((2n-1)!! = 1 \cdot 3 \cdot 5 \cdots (2n-1)\) は二重階乗です。
この構成により、\(s \to 0\) すなわち低周波数域での群遅延が最大限に平坦になります。
群遅延の最大平坦性
ベッセルフィルタの群遅延 \(\tau(\omega)\) を \(\omega=0\) でテイラー展開すると、\(2n-1\) 次まで導関数がゼロになります：
\[\tau(\omega) = \tau_0 \left[ 1 + O(\omega^{2n}) \right] \tag{6}\]これが「群遅延が最大平坦」の意味です。\(n\) が大きいほど平坦な領域が広がります。
他のIIRフィルタとの比較
| フィルタ | 振幅平坦性 | 遷移帯域の急峻さ | 群遅延特性 | 主な用途 |
|---|---|---|---|---|
| ベッセル | 最も劣る | 最もなだらか | 最大平坦 | 波形保持・パルス伝送 |
| バターワース | 最大平坦 | 中程度 | 比較的良好 | 汎用ローパス |
| チェビシェフI | 等リップル | バターワースより急峻 | 非線形 | 急峻な遮断が必要な場合 |
| 楕円（Cauer） | 等リップル | 最も急峻 | 最も劣る | 最小次数で仕様を満たす場合 |
ベッセルフィルタは振幅特性（遮断の鋭さ）を犠牲にして位相・群遅延特性を最大化します。遷移帯域が非常になだらかなため、同じ遮断仕様を満たすには高い次数が必要です。
SciPyによる実装
SciPyでは scipy.signal.bessel() 関数でベッセルフィルタを設計できます。デフォルトでは群遅延を \(\omega=0\) で正規化します（norm='delay'）。
基本的なローパスフィルタ設計
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== フィルタ設計パラメータ =====
fs = 1000 # サンプリング周波数 [Hz]
fc = 100 # カットオフ周波数 [Hz]
orders = [2, 4, 6]
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
for N in orders:
# ベッセルローパスフィルタ（群遅延正規化）
sos = signal.bessel(N, fc, btype='low', fs=fs, norm='delay', output='sos')
w, h = signal.sosfreqz(sos, worN=4096, fs=fs)
# 群遅延の計算
w_gd, gd = signal.group_delay((signal.bessel(N, fc, btype='low', fs=fs, norm='delay')), w=4096, fs=fs)
axes[0].plot(w, 20 * np.log10(np.abs(h) + 1e-12), label=f'N={N}')
axes[1].plot(w, np.degrees(np.unwrap(np.angle(h))), label=f'N={N}')
axes[2].plot(w_gd, gd, label=f'N={N}')
axes[0].set_xlim(0, 500)
axes[0].set_ylim(-60, 5)
axes[0].set_xlabel('周波数 [Hz]')
axes[0].set_ylabel('ゲイン [dB]')
axes[0].set_title('ベッセルフィルタ ゲイン特性')
axes[0].axvline(fc, color='gray', linestyle='--', alpha=0.5, label=f'fc={fc}Hz')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].set_xlim(0, 500)
axes[1].set_xlabel('周波数 [Hz]')
axes[1].set_ylabel('位相 [度]')
axes[1].set_title('ベッセルフィルタ 位相特性')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[2].set_xlim(0, fs / 2)
axes[2].set_ylim(0, 0.02)
axes[2].set_xlabel('周波数 [Hz]')
axes[2].set_ylabel('群遅延 [s]')
axes[2].set_title('ベッセルフィルタ 群遅延特性')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
他フィルタとの群遅延比較
ベッセルフィルタの群遅延平坦性を他フィルタと視覚的に比較します：
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000
fc = 100
N = 4
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
filters = {
'Bessel': signal.bessel(N, fc, btype='low', fs=fs, norm='delay'),
'Butterworth': signal.butter(N, fc, btype='low', fs=fs),
'Chebyshev I (rp=1dB)': signal.cheby1(N, 1.0, fc, btype='low', fs=fs),
'Elliptic (rp=1,rs=60)': signal.ellip(N, 1.0, 60.0, fc, btype='low', fs=fs),
}
colors = ['royalblue', 'tomato', 'green', 'purple']
for (name, (b, a)), color in zip(filters.items(), colors):
w, h = signal.freqz(b, a, worN=4096, fs=fs)
w_gd, gd = signal.group_delay((b, a), w=4096, fs=fs)
axes[0].plot(w, 20 * np.log10(np.abs(h) + 1e-12), label=name, color=color)
axes[1].plot(w_gd, gd, label=name, color=color)
axes[0].set_xlim(0, 500)
axes[0].set_ylim(-80, 5)
axes[0].set_xlabel('周波数 [Hz]')
axes[0].set_ylabel('ゲイン [dB]')
axes[0].set_title(f'ゲイン特性の比較（N={N}）')
axes[0].axvline(fc, color='gray', linestyle='--', alpha=0.5, label=f'fc={fc}Hz')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].set_xlim(0, 500)
axes[1].set_ylim(0, 0.025)
axes[1].set_xlabel('周波数 [Hz]')
axes[1].set_ylabel('群遅延 [s]')
axes[1].set_title(f'群遅延特性の比較（N={N}）')
axes[1].axvline(fc, color='gray', linestyle='--', alpha=0.5)
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
ハイパス・バンドパスフィルタへの拡張
btype パラメータを変更するだけでハイパス・バンドパスフィルタを設計できます：
from scipy import signal
fs = 1000
N = 4
# ===== ハイパスフィルタ =====
sos_hp = signal.bessel(N, 200, btype='high', fs=fs, norm='delay', output='sos')
print("ハイパスフィルタ: 設計完了")
# ===== バンドパスフィルタ =====
sos_bp = signal.bessel(N, [100, 300], btype='bandpass', fs=fs, norm='delay', output='sos')
print("バンドパスフィルタ: 設計完了")
# ===== バンドストップフィルタ =====
sos_bs = signal.bessel(N, [100, 300], btype='bandstop', fs=fs, norm='delay', output='sos')
print("バンドストップフィルタ: 設計完了")
パルス波形への適用例
ベッセルフィルタの最大の利点を示す例として、矩形パルス信号を各フィルタに通した場合の波形比較を示します：
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000
fc = 100
N = 4
t = np.linspace(0, 0.1, int(fs * 0.1), endpoint=False)
# 矩形パルス信号（幅5ms）
x = np.zeros(len(t))
x[(t >= 0.02) & (t < 0.025)] = 1.0
filters = {
'Bessel': signal.bessel(N, fc, btype='low', fs=fs, norm='delay', output='sos'),
'Butterworth': signal.butter(N, fc, btype='low', fs=fs, output='sos'),
'Chebyshev I': signal.cheby1(N, 1.0, fc, btype='low', fs=fs, output='sos'),
'Elliptic': signal.ellip(N, 1.0, 60.0, fc, btype='low', fs=fs, output='sos'),
}
fig, axes = plt.subplots(len(filters) + 1, 1, figsize=(10, 14), sharex=True)
axes[0].plot(t * 1000, x, 'k', label='入力パルス')
axes[0].set_ylabel('振幅')
axes[0].set_title('入力信号（矩形パルス）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
for idx, (name, sos) in enumerate(filters.items()):
y = signal.sosfilt(sos, x) # lfilter相当（因果的、リアルタイム）
axes[idx + 1].plot(t * 1000, y, label=name)
axes[idx + 1].plot(t * 1000, x, 'k--', alpha=0.3, label='元のパルス')
axes[idx + 1].set_ylabel('振幅')
axes[idx + 1].set_title(f'{name} フィルタ後（N={N}）')
axes[idx + 1].legend()
axes[idx + 1].grid(True, alpha=0.3)
axes[-1].set_xlabel('時間 [ms]')
plt.tight_layout()
plt.show()
このコードを実行すると、ベッセルフィルタが最も波形の形状を保持していること、楕円フィルタでは大きなリンギングが発生することが確認できます。
norm パラメータの意味
scipy.signal.bessel() には norm パラメータがあり、正規化方法を選べます：
norm 値 | 意味 | 用途 |
|---|---|---|
'delay' | \(\omega=0\) での群遅延が \(1/\omega_c\) になるよう正規化 | 群遅延の比較（デフォルト） |
'critical' | -3dBカットオフが指定周波数になるよう正規化 | 他フィルタとの振幅比較 |
'mag' | 通過域端のゲインを指定値に合わせる | 特定ゲイン要件がある場合 |
他のフィルタ（バターワース等）と同じカットオフ周波数で比較したい場合は norm='critical' を使用します：
from scipy import signal
fs = 1000
fc = 100
N = 4
# -3dBカットオフをfc=100Hzに設定
sos_bessel = signal.bessel(N, fc, btype='low', fs=fs, norm='critical', output='sos')
sos_butter = signal.butter(N, fc, btype='low', fs=fs, output='sos')
ベッセルフィルタを選ぶべきシナリオ
| シナリオ | 推奨フィルタ |
|---|---|
| パルス波形の形状を保持したい | ベッセルフィルタ |
| 音声信号の位相歪みを最小化したい | ベッセルフィルタ |
| 遷移帯域をできるだけ急峻にしたい | 楕円フィルタ |
| 通過域の振幅を最大限に平坦にしたい | バターワースフィルタ |
| 計算コスト（次数）を最小化したい | 楕円フィルタ |
| リップルを許容して遮断を急峻にしたい | チェビシェフフィルタ |
関連記事
- バターワースフィルタの設計原理とPython実装 - 通過域最大平坦特性を持つIIRフィルタの設計
- チェビシェフフィルタの設計原理とPython実装 - 等リップルIIRフィルタの設計
- 楕円フィルタ（Elliptic Filter）の設計原理とPython実装 - 最急峻な遷移帯域を持つIIRフィルタの設計
- ローパスフィルタの設計と比較：移動平均・バターワース・チェビシェフ - 各種ローパスフィルタの周波数応答比較
- FIRフィルタとIIRフィルタの比較 - IIRフィルタ全般の位置づけと特性比較
- 指数移動平均（EMA）フィルタの周波数特性 - シンプルな1次IIRフィルタとしてのEMAの解析
- 適応フィルタ（LMS/RLS）の理論とPython実装 - パラメータが自動調整される適応フィルタ
参考文献
- Thomson, W. E. (1952). Delay networks having maximally flat frequency characteristics. Proceedings of the IEE - Part III: Radio and Communication Engineering, 96(44), 487–490.
- Proakis, J. G., & Manolakis, D. G. (2007). Digital Signal Processing (4th ed.). Pearson.
- SciPy scipy.signal.bessel documentation
関連ツール
- DevToolBox - 開発者向け無料ツール集 - JSON整形、正規表現テスターなど85種類以上の開発者向けツール
- CalcBox - 暮らしの計算ツール - 統計計算、複利計算など61種類以上の計算ツール