バターワースフィルタの設計原理とPython実装

バターワースフィルタの最大平坦特性の数学的導出から、アナログ・デジタル変換、SciPyによる設計・適用までを解説します。

3 min read

はじめに

バターワースフィルタは、通過域で最大限に平坦な振幅特性を持つIIRフィルタです。1930年にStephen Butterworthが発表したこのフィルタは、リップル（波打ち）のない滑らかな周波数応答を特徴とし、信号処理で最も広く使われるフィルタ設計の一つです。

本記事では、バターワースフィルタの数学的基礎から、SciPyを使った実践的な設計・実装までを解説します。

振幅特性

振幅二乗関数

\(N\) 次バターワースローパスフィルタの振幅二乗関数は次のように定義されます。

\[|H(j\Omega)|^2 = \frac{1}{1 + \left(\frac{\Omega}{\Omega_c}\right)^{2N}} \tag{1}\]

ここで \(\Omega_c\) はカットオフ角周波数（-3dBの周波数）、\(N\) はフィルタ次数です。

最大平坦特性

式 \((1)\) を \(\Omega = 0\) の周りでテイラー展開すると、\(|H(j\Omega)|^2\) の最初の \(2N-1\) 個の導関数がすべてゼロになることが示されます。これが「最大平坦（maximally flat）」と呼ばれる理由です。

周波数による振る舞い

  • \(\Omega = 0\): \(|H| = 1\)（完全通過）
  • \(\Omega = \Omega_c\): \(|H| = 1/\sqrt{2} \approx -3\,\text{dB}\)
  • \(\Omega \gg \Omega_c\): \(|H| \approx (\Omega_c/\Omega)^N\)（ロールオフ率 \(-20N\) dB/decade）

極の配置

バターワースフィルタの極は、s平面上でバターワース円と呼ばれる半径 \(\Omega_c\) の円上に等間隔に配置されます。\(N\) 次フィルタの \(k\) 番目の極は次のとおりです。

\[s_k = \Omega_c \exp\left(j\frac{\pi(2k + N - 1)}{2N}\right), \quad k = 1, 2, \ldots, N \tag{2}\]

安定なフィルタを得るため、左半平面（実部が負）にある極のみを使用します。

設計手順

1. フィルタ次数の決定

通過域・阻止域の仕様から必要な次数を求めます。

  • 通過域: \(|H(j\Omega_p)| \geq 1/\sqrt{1+\varepsilon_p^2}\)（通過域リップル）
  • 阻止域: \(|H(j\Omega_s)| \leq 1/\sqrt{1+\varepsilon_s^2}\)（阻止域減衰）

最小次数は次の式で求まります。

\[N \geq \frac{\log(\varepsilon_s^2/\varepsilon_p^2)}{2\log(\Omega_s/\Omega_p)} \tag{3}\]

2. アナログプロトタイプの設計

正規化カットオフ周波数 \(\Omega_c = 1\) のアナログバターワースフィルタを設計します。

3. デジタルフィルタへの変換

**双一次変換（bilinear transform）**を用いてアナログフィルタをデジタルフィルタに変換します。

\[s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}} \tag{4}\]

双一次変換では周波数軸の歪み（ワーピング）が生じるため、事前にアナログ周波数を補正（プリワーピング）する必要があります。

\[\Omega_a = \frac{2}{T}\tan\left(\frac{\omega_d T}{2}\right) \tag{5}\]

Python実装

SciPyによるフィルタ設計と適用

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# --- フィルタ設計 ---
fs = 1000           # サンプリング周波数 [Hz]
fc = 100            # カットオフ周波数 [Hz]
orders = [2, 4, 8]  # フィルタ次数

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

for N in orders:
    sos = signal.butter(N, fc, fs=fs, output='sos')
    w, h = signal.sosfreqz(sos, worN=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[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Butterworth Filter - Magnitude Response')
axes[0].set_xlim(0, 500)
axes[0].set_ylim(-80, 5)
axes[0].axvline(fc, color='gray', linestyle='--', alpha=0.5)
axes[0].axhline(-3, color='gray', linestyle=':', alpha=0.5, label='-3 dB')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 位相特性のプロット設定
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [degrees]')
axes[1].set_title('Butterworth Filter - Phase Response')
axes[1].set_xlim(0, 500)
axes[1].legend()
axes[1].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
t = np.arange(0, 1, 1/fs)
clean = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 30 * t)
noisy = clean + 0.8 * np.random.randn(len(t))

# --- バターワースフィルタの設計と適用 ---
sos = signal.butter(4, 50, fs=fs, output='sos')
filtered_causal = signal.sosfilt(sos, noisy)       # 因果フィルタ
filtered_zero = signal.sosfiltfilt(sos, noisy)      # ゼロ位相フィルタ

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)

axes[0].plot(t, noisy, alpha=0.5, label='Noisy')
axes[0].plot(t, clean, 'k', linewidth=1.5, label='Original')
axes[0].set_title('Input Signal')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, filtered_causal, label='sosfilt (causal)')
axes[1].plot(t, clean, 'k', linewidth=1.5, label='Original')
axes[1].set_title('Causal Filtering (sosfilt)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t, filtered_zero, label='sosfiltfilt (zero-phase)')
axes[2].plot(t, clean, 'k', linewidth=1.5, label='Original')
axes[2].set_title('Zero-Phase Filtering (sosfiltfilt)')
axes[2].set_xlabel('Time [s]')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

sosfilt は因果フィルタ（位相遅れあり）、sosfiltfilt は前方・後方の2回フィルタリングによりゼロ位相を実現しますが、リアルタイム処理には使えません。

他のIIRフィルタとの比較

フィルタ通過域阻止域遷移帯域位相特性
バターワース最大平坦単調減少広い比較的滑らか
チェビシェフI型等リップル単調減少狭い急峻な変化
チェビシェフII型単調減少等リップル狭い比較的滑らか
楕円（Cauer）等リップル等リップル最も狭い最も急峻

同じ次数では楕円フィルタが最も急峻な遷移帯域を持ちますが、通過域にリップルが生じます。リップルが許容できない用途ではバターワースが最適です。

関連記事

参考文献

  • Butterworth, S. (1930). “On the theory of filter amplifiers”. Wireless Engineer, 7(6), 536-541.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • SciPy scipy.signal.butter documentation

関連記事