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

はじめに

バターワースフィルタは、通過域で最大限に平坦な振幅特性を持つ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}\]

数値検証：プリワーピングを怠ると次数を過大評価する

式 \((3)\) はアナログ角周波数の比

\[ \Omega_s/\Omega_p \]

を使う式ですが、デジタルフィルタの設計では実周波数の比をそのまま代入する誤りがよく起こります。プリワーピングの有無で必要次数がどれだけ変わるかを、具体的な仕様で確認します。

仕様: サンプリング周波数 \(f_s = 1000\) Hz、通過域端 \(f_p = 100\) Hz、阻止域端 \(f_{st} = 300\) Hz、通過域リップル \(R_p = 1\) dB、阻止域減衰量 \(R_s = 40\) dB。

import numpy as np
from scipy import signal

fs = 1000                 # サンプリング周波数 [Hz]
fp, fst = 100, 300         # 通過域端 / 阻止域端 [Hz]
rp, rs = 1, 40             # 通過域リップル [dB] / 阻止域減衰量 [dB]

eps_p = np.sqrt(10**(rp / 10) - 1)
eps_s = np.sqrt(10**(rs / 10) - 1)

# (a) プリワーピングなし：実周波数の比をそのまま式(3)に代入（誤った適用）
N_naive = np.log10(eps_s**2 / eps_p**2) / (2 * np.log10(fst / fp))

# (b) プリワーピングあり：tanで角周波数を補正してから比をとる（正しい適用）
T = 1 / fs
Wp = (2 / T) * np.tan(2 * np.pi * fp * T / 2)
Wst = (2 / T) * np.tan(2 * np.pi * fst * T / 2)
N_prewarp = np.log10(eps_s**2 / eps_p**2) / (2 * np.log10(Wst / Wp))

# (c) scipy.signal.buttord による自動設計との照合
N_buttord, wn = signal.buttord(fp / (fs / 2), fst / (fs / 2), rp, rs)

print(f"プリワーピングなし: N = {N_naive:.3f} -> 採用次数 {int(np.ceil(N_naive))}")
print(f"プリワーピングあり: N = {N_prewarp:.3f} -> 採用次数 {int(np.ceil(N_prewarp))}")
print(f"buttord:            N = {N_buttord} (Wn = {wn:.4f})")

実行結果:

プリワーピングなし: N = 4.807 -> 採用次数 5
プリワーピングあり: N = 3.658 -> 採用次数 4
buttord:            N = 4 (Wn = 0.2338)

プリワーピングを省略すると、実周波数比がそのままの比

\[ f_{st}/f_p = 3.0 \]

倍として扱われ、必要次数を 5次と過大評価します。一方、\(\tan\) によるプリワーピングを行うと阻止域端の角周波数がナイキスト周波数（\(500\) Hz）に近いぶん大きく引き伸ばされ、実効的な周波数比は

\[ W_{st}/W_p = 2752.76 / 649.84 = 4.24 \]

まで拡大します。この拡大された比を式 \((3)\) に使うと \(N=3.658\) となり、切り上げて 4次で仕様を満たすと分かります。これは scipy.signal.buttord の結果と完全に一致します。実際に4次フィルタで仕様を検証すると、\(f_p=100\) Hz で \(-1.00\) dB（要求 \(R_p=1\) dB ちょうど）、\(f_{st}=300\) Hz で \(-44.29\) dB（要求 \(R_s=40\) dB に対して \(4.29\) dB の余裕）となり、確かに仕様を満たしています。

注意点: scipy.signal.buttord はこのプリワーピングを内部で自動的に行うため、実務では式 \((3)\) を手計算する必要はほとんどありません。しかし、手計算やスプレッドシートで次数を見積もる場合、この補正を忘れると次数を1つ余分に見積もり、不要に高次（＝不要に高い計算コスト・数値誤差リスク）のフィルタを設計してしまう点に注意してください。この効果は阻止域端がナイキスト周波数に近いほど（つまりサンプリング周波数に対して比較的高い周波数を扱うほど）顕著になります。

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()

次数ごとのゲイン・ロールオフの数値検証

「カットオフ周波数で-3dB」「ロールオフは\(-20N\) dB/decade」という性質を、実際に scipy.signal で数値検証します。\(f_c = 100\) Hz、\(N \in \{2, 4, 6, 8\}\) で設計し、カットオフの2倍・10倍・100倍の周波数でのゲインを計算します。

import numpy as np
from scipy import signal

fs = 100_000   # 高域まで確認するためサンプリング周波数を高く取る
fc = 100       # カットオフ周波数 [Hz]
orders = [2, 4, 6, 8]

print(f"{'N':<4}{'gain@fc[dB]':>14}{'gain@2fc[dB]':>15}{'gain@10fc[dB]':>16}{'gain@100fc[dB]':>17}")
for N in orders:
    sos = signal.butter(N, fc, fs=fs, output='sos')
    freqs = np.array([fc, 2 * fc, 10 * fc, 100 * fc])
    _, h = signal.sosfreqz(sos, worN=freqs, fs=fs)
    gdb = 20 * np.log10(np.abs(h))
    print(f"{N:<4}{gdb[0]:>14.4f}{gdb[1]:>15.4f}{gdb[2]:>16.4f}{gdb[3]:>17.4f}")

実行結果:

N      gain@fc[dB]   gain@2fc[dB]  gain@10fc[dB]  gain@100fc[dB]
2          -3.0103       -12.3047       -40.0061        -80.5850
4          -3.0103       -24.0997       -80.0113       -161.1700
6          -3.0103       -36.1252      -120.0170       -241.7550
8          -3.0103       -48.1656      -160.0226       -322.3400

\(f_c\) における減衰量はどの次数でも \(-3.0103\) dB とまったく同一であることが確認できます（式 \((1)\) から \(|H(j\Omega_c)|^2 = 1/2\) となるのは \(N\) に依存しないため）。周波数特性を次数別に重ねてプロットすると、ロールオフの傾きの違いが視覚的にも明確です。

バターワースフィルタの次数別周波数特性（fc=100Hz）

次に、\(10f_c\) から \(100f_c\) の区間（1デカード）でのゲイン変化量から実測ロールオフを求め、理論値 \(-20N\) dB/decade と比較します。

import numpy as np
from scipy import signal

fs = 100_000
fc = 100
orders = [2, 4, 6, 8]

theoretical = [-20 * N for N in orders]
measured = []
for N in orders:
    sos = signal.butter(N, fc, fs=fs, output='sos')
    freqs = np.array([10 * fc, 100 * fc])
    _, h = signal.sosfreqz(sos, worN=freqs, fs=fs)
    gdb = 20 * np.log10(np.abs(h))
    slope = (gdb[1] - gdb[0]) / (np.log10(100 * fc) - np.log10(10 * fc))
    measured.append(slope)

for N, th, ms in zip(orders, theoretical, measured):
    print(f"N={N}: 理論値 {th} dB/decade, 実測値 {ms:.3f} dB/decade, 差分 {ms - th:+.3f} dB")

実行結果:

N=2: 理論値 -40 dB/decade, 実測値 -40.579 dB/decade, 差分 -0.579 dB
N=4: 理論値 -80 dB/decade, 実測値 -81.159 dB/decade, 差分 -1.159 dB
N=6: 理論値 -120 dB/decade, 実測値 -121.738 dB/decade, 差分 -1.738 dB
N=8: 理論値 -160 dB/decade, 実測値 -162.317 dB/decade, 差分 -2.317 dB

次数ごとのロールオフ勾配：理論値と実測値の比較

エッジケース: 次数が上がるほど理論値からのずれが大きくなっています。これはアナログプロトタイプが正確に \(-20N\) dB/decade に従うのに対し、双一次変換によるデジタル化では高周波側の周波数軸が対数的に圧縮される（ワーピングの影響）ためです。実際、同じフィルタを signal.butter(..., analog=True) でアナログ設計のまま評価すると、\(10f_c\) 〜\(100f_c\) の勾配は \(N=2,4,6,8\) に対してそれぞれ厳密に \(-40.000, -80.000, -120.000, -160.000\) dB/decade となり、理論値と完全に一致します。つまりここで見えているズレはバターワースの理論の誤差ではなく、ナイキスト周波数に近い帯域を扱う際にデジタルIIRフィルタが本質的に持つ非線形な周波数軸の歪みであり、高次フィルタやサンプリング周波数に対して比較的高いカットオフ周波数を設計する際に注意すべき点です。

ノイズ除去の応用例

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）	等リップル	等リップル	最も狭い	最も急峻

同じ次数では楕円フィルタが最も急峻な遷移帯域を持ちますが、通過域にリップルが生じます。リップルが許容できない用途ではバターワースが最適です。他フィルタとの次数・群遅延の定量比較（buttord・cheb1ord・cheb2ord・ellipord による最小次数の横並び比較を含む）は楕円フィルタ（Elliptic Filter）の理論とPython実装にまとめているので、本記事ではバターワース単体の設計原理に絞って解説しています。

近年の研究動向

バターワースフィルタは1930年発表の古典的な設計ですが、現在も次数決定アルゴリズムや実装形態を巡る研究が続いています。Dey, Roy & Sarkar (2026) は、通常の整数次数バターワースフィルタを一般化した分数次（fractional-order）バターワースフィルタの設計に、カオス写像と準対立学習（quasi-oppositional learning）を組み合わせた新しいメタヒューリスティック最適化アルゴリズム「Chaotic Quasi-Oppositional Bald Eagle Search（CQOBES）」を適用し、ローパスだけでなくハイパス・バンドパス型の分数次バターワースフィルタでも目標特性への誤差を最小化できることを示しました（Circuits, Systems, and Signal Processing, Springer, 2026）。分数次フィルタは通過域・阻止域の間の遷移特性を整数次数よりも細かく制御できるため、通常のバターワース設計では次数を1つ上げるかどうかの二択しかない場面で、より仕様に近い特性を実現できる点が実用上のメリットです。

参考文献

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.
Dey, S., Roy, P. K., & Sarkar, A. (2026). “Accurate Design of Digital Fractional-Order Butterworth Filters Using a Novel Chaotic Quasi-Oppositional Bald Eagle Search Algorithm”. Circuits, Systems, and Signal Processing. https://link.springer.com/article/10.1007/s00034-026-03623-1
SciPy scipy.signal.butter documentation

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

はじめに

振幅特性

振幅二乗関数

最大平坦特性

周波数による振る舞い

極の配置

設計手順

1. フィルタ次数の決定

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

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

数値検証：プリワーピングを怠ると次数を過大評価する

Python実装

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

次数ごとのゲイン・ロールオフの数値検証

ノイズ除去の応用例

他のIIRフィルタとの比較

近年の研究動向

おすすめ書籍

関連記事

参考文献

はじめに

振幅特性

振幅二乗関数

最大平坦特性

周波数による振る舞い

極の配置

設計手順

1. フィルタ次数の決定

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

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

数値検証：プリワーピングを怠ると次数を過大評価する

Python実装

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

次数ごとのゲイン・ロールオフの数値検証

ノイズ除去の応用例

他のIIRフィルタとの比較

近年の研究動向

おすすめ書籍

関連記事

参考文献

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

関連記事

非線形カルマンスムーザ：Extended RTS SmootherとUnscented RTS Smootherの理論とPython実装

テーガー・カイザーエネルギー演算子（TKEO）の理論とPython実装：ヒルベルト変換との比較・軸受診断への応用

カルマンスムーザ3種比較（RTS固定区間・固定ラグ・固定点）：理論・遅延と計算量のトレードオフ・Python実装