フィルタ設計で最初に名前が挙がるのがバターワースフィルタです。「最大平坦特性」という言葉が示すとおり、通過域の振幅ができるだけ平らになるよう設計されており、穏やかなロールオフと安定した動作が特徴です。スペックアナライザのIF回路や医療用センサの前処理、音響イコライザなど、「余計な色付けをしたくない」場面で広く採用されています。
本記事では、バターワースフィルタを「知っている」段階から「設計できる」段階へ引き上げることを目標にします。振幅二乗特性の導出・極配置・双一次変換によるデジタル化という設計の流れを数式とPythonで丁寧に追い、最後にリアルな信号へのノイズ除去を実装します。
応用先の例:
- 生体信号の前処理 — 心電図（ECG）や脳波（EEG）の記録では、電源ノイズ（50/60 Hz）や筋電ノイズを除去する前処理フィルタにバターワースLPFが定番です。通過域が平坦なので、病的なピークの形状を歪めることなく低周波成分を抽出できます。
- 音響・通信システムのIF段 — 受信機のIF（中間周波）フィルタや音声帯域フィルタでは、通過域内の信号振幅を均等に扱いたいため、フラットな振幅特性を持つバターワースが選ばれます。
本記事の内容:
- バターワースフィルタの最大平坦特性とは何か（直感と数式）
- 振幅二乗特性 $|H(j\omega)|^2 = 1/(1+(\omega/\omega_c)^{2n})$ の導出
- 次数 $n$ とロールオフ（$-20n$ dB/dec）の関係
- 極配置（s平面の左半円に等角度間隔）
- アナログプロトタイプから双一次変換でデジタル化する手順
- 他フィルタ（チェビシェフ・ベッセル）との比較
scipy.signal.butterを使った設計と実信号ノイズ除去
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
バターワースフィルタとは — 最大平坦特性の直感
コーヒーサイフォンで細かい粉を濾しても、コーヒーの味が変わらないようにしたい——これが「最大平坦特性」の直感です。通過させたい帯域（通過域）では、どの周波数成分もほぼ同じ割合で通し、遮断したい帯域（阻止域）ではなるべく急峻に落としたい。
バターワースフィルタはこのうち「通過域をできるだけ平らにする」ことを最優先にした設計です。1930年にイギリスの電気工学者スティーブン・バターワースが提案したこのフィルタは、「最大平坦マグニチュードフィルタ（Maximally Flat Magnitude Filter）」とも呼ばれます。
他のフィルタと比べると次の図のような位置づけになります。
この図から、バターワースフィルタの特徴が一目でわかります。チェビシェフ Type I は遷移域がバターワースより急峻ですが、通過域に「リプル（波打ち）」があり、信号の振幅が周波数によってわずかに変動します。ベッセルフィルタはロールオフが最も緩やかですが、群遅延が一定（位相線形性が高い）という別の強みを持ちます。バターワースはこの三者の中間——通過域はリプルのない完全なフラット、ロールオフはチェビシェフより緩やか、というバランスの取れた選択肢です。
では、「通過域を平らにする」という条件を数学的にどう表現するのでしょうか。それを示すのが次の振幅二乗特性です。
振幅二乗特性の導出
最大平坦条件の定式化
「通過域が最大平坦」とは、$\omega = 0$ 付近で振幅の微分（=$\omega$ に対する変化率）が可能な限り多くの次数でゼロになることを意味します。これを「テイラー展開の低次項がすべて消える」という視点で捉えると、設計条件が明確になります。
まず $n$ 次の低域通過フィルタ（LPF）の振幅二乗関数を多項式で書き下します。実現可能なフィルタの振幅二乗関数 $|H(j\omega)|^2$ は $\omega^2$ の有理多項式です。規格化遮断周波数 $\omega_c = 1$ で、直流利得を 1 ($|H(0)|^2 = 1$) に揃えると:
$$ |H(j\omega)|^2 = \frac{1}{1 + c_1\omega^2 + c_2\omega^4 + \cdots + c_n\omega^{2n}} $$
$\omega \to 0$ で関数値が 1 になるよう定数項を 1 に取りました。ここで「最大平坦条件」を課します——$\omega=0$ での微分係数 $\frac{d^k |H|^2}{d(\omega^2)^k}\big|_{\omega=0}$ が $k=1,2,\ldots,n-1$ についてゼロになる、という要求です。
$\omega^2 = u$ と置くと、$|H|^2 = 1/(1 + c_1 u + c_2 u^2 + \cdots)$ です。$u=0$ での $k$ 階微分がゼロになるためには、$c_1 = c_2 = \cdots = c_{n-1} = 0$ でなければなりません（1次から $(n-1)$ 次の係数はすべて消える）。
これを代入すると分母に残るのは $c_n \omega^{2n}$ 項のみです。$\omega = \omega_c$ で振幅が $1/\sqrt{2}$（パワーが半分、つまり $-3$ dB）になるという遮断条件を加えると:
$$ |H(j\omega_c)|^2 = \frac{1}{1 + c_n} = \frac{1}{2} $$
よって $c_n = 1$ です。
バターワース振幅二乗特性
以上の導出により、バターワースフィルタの振幅二乗特性が定まります:
$$ \boxed{|H(j\omega)|^2 = \frac{1}{1 + \left(\dfrac{\omega}{\omega_c}\right)^{2n}}} $$
この式には美しい性質が詰まっています。
$\omega \ll \omega_c$（通過域）: $(\omega/\omega_c)^{2n} \approx 0$ なので $|H|^2 \approx 1$。振幅はほぼ 1 で平坦。
$\omega = 0$ でのテイラー展開を確認します。$u = (\omega/\omega_c)^{2n}$ として:
$$ |H|^2 = 1 – u + u^2 – \cdots $$
$u = (\omega/\omega_c)^{2n}$ は $\omega$ の $2n$ 次の項から始まるので、$\omega^2, \omega^4, \ldots, \omega^{2n-2}$ の係数がすべてゼロ——これが「最大平坦」の数学的意味です。
$\omega = \omega_c$（遮断周波数）: $(\omega/\omega_c)^{2n} = 1$ なので:
$$ |H(j\omega_c)|^2 = \frac{1}{2} \implies |H(j\omega_c)| = \frac{1}{\sqrt{2}} \approx 0.707 $$
パワーが半分（$-3$ dB）になる周波数が正確に $\omega_c$ です。
$\omega \gg \omega_c$（阻止域）: $(\omega/\omega_c)^{2n} \gg 1$ なので:
$$ |H(j\omega)|^2 \approx \left(\frac{\omega_c}{\omega}\right)^{2n} $$
対数表示すると、1 decade（10倍）周波数が上がるごとに振幅は $-20n$ dB ずつ落ちます——これがロールオフ特性です。
次の図で、次数 $n$ による振幅特性の変化を確認します。
次数 $n$ が大きくなるにつれ、通過域（$\omega < \omega_c$）では限りなく平坦に近づき、阻止域（$\omega > \omega_c$）ではより急峻に減衰します。$n \to \infty$ の極限では、$\omega_c$ を境とする完全な矩形窓——つまり理想ローパスフィルタ——に収束します。これは数式からも直感的にわかります。$\omega < \omega_c$ ならば $(\omega/\omega_c)^{2n} \to 0$、$\omega > \omega_c$ ならば $(\omega/\omega_c)^{2n} \to \infty$ という対称性です。ただし次数を上げるほど位相歪みや群遅延ピークも大きくなるため、実用では $n = 4 \sim 8$ 程度が多く使われます。
振幅特性の形がわかったので、次はこのフィルタを実際にどう実現するか——すなわちアナログ伝達関数 $H(s)$ の極配置を見ていきます。
極配置 — s平面の左半円
振幅二乗特性から伝達関数へ
アナログフィルタの伝達関数 $H(s)$ は、$s = j\omega$ を代入したとき振幅二乗特性を与える必要があります:
$$ H(s)\,H(-s)\Big|_{s=j\omega} = |H(j\omega)|^2 = \frac{1}{1 + \left(\dfrac{\omega}{\omega_c}\right)^{2n}} $$
$j\omega$ を $s$ に置き換えると（$\omega = s/j = -js$）:
$$ H(s)\,H(-s) = \frac{1}{1 + \left(\dfrac{s}{j\omega_c}\right)^{2n}} = \frac{1}{1 + \left(\dfrac{-s^2}{\omega_c^2}\right)^{n}} $$
$(j)^{2n} = (-1)^n$ を使うと:
$$ H(s)\,H(-s) = \frac{\omega_c^{2n}}{\omega_c^{2n} + (-1)^n s^{2n}} $$
この分母の極（$s$ の値）を求めます。分母をゼロとおくと:
$$ 1 + \frac{(-1)^n s^{2n}}{\omega_c^{2n}} = 0 \implies s^{2n} = (-1)^{n+1} \omega_c^{2n} $$
$2n$ 個の根を複素数で書くと:
$$ s_k = \omega_c \exp\!\left(j\frac{(2k-1)\pi}{2n} + j\frac{\pi}{2}\right) = \omega_c \exp\!\left(j\frac{(2k-1+n)\pi}{2n}\right), \quad k = 1, 2, \ldots, 2n $$
これらの極は 単位円上（$|s_k| = \omega_c$）に等角度間隔 で並びます。隣り合う極同士の角度差は $\pi/n$ ラジアン（$180°/n$ 度）です。
安定フィルタのための極選択
$H(s)\,H(-s)$ の $2n$ 個の極のうち、左半平面（実部が負） に位置する $n$ 個だけを $H(s)$ に割り当てます。右半平面の極を使うと因果システムが不安定になるためです。
左半平面の $n$ 個の極は:
$$ s_k = \omega_c \exp\!\left(j\theta_k\right), \quad \theta_k = \frac{(2k-1)\pi}{2n} + \frac{\pi}{2}, \quad k = 1, 2, \ldots, n $$
実際に $\theta_k$ の範囲を確認すると、$k=1$ のとき $\theta_1 = \pi/2 + \pi/(2n)$（第二象限）、$k=n$ のとき $\theta_n = \pi/2 + \pi(2n-1)/(2n)$（第三象限）となり、確かに左半平面に収まります。
これを図で確認します。
この図から、いくつかの重要な事実が読み取れます。第一に、バターワース極はすべて単位円（$|s|=\omega_c=1$）上に位置し、左半平面に等角度で分布しています。第二に、次数が偶数（$n=4$）の場合は実軸上に極がなく、すべての極が共役複素数のペアになります。次数が奇数（$n=3, 5$）の場合は $s=-\omega_c$（実軸上）に一つの実極が現れます。第三に、極が虚軸から遠いほど応答が安定で穏やかになります——これがバターワースがオーバーシュートの少ない応答を見せる理由です。
正規化されたバターワース伝達関数（$\omega_c = 1$ の場合）は:
$$ H(s) = \frac{1}{\prod_{k=1}^{n}(s – s_k)} $$
例えば $n=2$ の場合、極は $s_{1,2} = e^{j3\pi/4}, e^{j5\pi/4}$ であり:
$$ H_2(s) = \frac{1}{(s – e^{j3\pi/4})(s – e^{j5\pi/4})} = \frac{1}{s^2 + \sqrt{2}\,s + 1} $$
分母の係数 $\sqrt{2}$ は $2\cos(\pi/4)$ から来ており、次数 $n$ に応じた「バターワース多項式（Butterworth polynomial）」の係数です。
次は、このアナログ伝達関数をデジタルフィルタに変換する手順——双一次変換——を見ます。
デジタル化 — 双一次変換
なぜアナログからデジタルへ変換するのか
実際の信号処理システムではコンピュータで離散時間信号を扱うため、上で設計したアナログフィルタ $H(s)$ をデジタルフィルタ $H(z)$ に変換する必要があります。最も一般的な方法が双一次変換（Bilinear Transform）です。
双一次変換の考え方
連続時間積分 $\int$ を台形則（Tustin法）で近似することに相当する変換です。サンプリング周期 $T_s = 1/f_s$ を使って:
$$ s = \frac{2}{T_s}\,\frac{z – 1}{z + 1} = \frac{2f_s(z-1)}{z+1} $$
この変換の逆方向を見ると:
$$ z = \frac{1 + s T_s/2}{1 – s T_s/2} $$
$|z| = 1$（単位円）に対応するのが $s = j\omega$ 方向（虚軸）ですが、対応関係は線形ではなく周波数のゆがみ（frequency warping）が生じます。アナログ周波数 $\Omega$ とデジタル周波数 $\omega$ の関係は:
$$ \Omega = \frac{2}{T_s}\tan\!\left(\frac{\omega}{2}\right) $$
このゆがみを補正するため、設計の際に目標デジタル周波数 $\omega_d$ を事前にアナログ周波数 $\Omega_a$ へ「事前歪み（pre-warping）」します:
$$ \Omega_a = \frac{2}{T_s}\tan\!\left(\frac{\omega_d}{2}\right) $$
こうしてプリウォープ済みのアナログフィルタを設計した後、双一次変換を適用することで、遮断周波数が正確にデジタル領域の目標値に一致するデジタルフィルタが得られます。
設計手順のまとめ
デジタルバターワースフィルタの設計は次の4ステップで完成します:
- 目標を決める: サンプリング周波数 $f_s$、遮断周波数 $f_c$、次数 $n$
- デジタル遮断周波数を規格化: $\omega_d = 2\pi f_c / f_s$
- プリウォーピング: $\Omega_a = 2f_s \tan(\omega_d/2)$
- アナログプロトタイプを設計して双一次変換: $H(z) = H_a(s)\big|_{s = \frac{2f_s(z-1)}{z+1}}$
scipy では butter(n, Wn) がこのすべてを自動で処理します。
位相特性と群遅延
振幅特性だけでなく、位相特性も信号の品質に大きく影響します。位相歪みが大きいと、信号の波形が変形してしまいます。
位相応答を見ると、次数が大きいほど遮断周波数付近での位相の変化が急峻になることがわかります。低次（n=1）では位相変化が緩やかですが、高次では遮断周波数を境に位相が急激に変化します。このように位相が周波数に対して非線形に変化することが、バターワースフィルタの位相線形性が低い理由です。
位相の非線形性をより直接的に表すのが群遅延（Group Delay）です。群遅延は位相の周波数微分の負値として定義されます:
$$ \tau(\omega) = -\frac{d\phi(\omega)}{d\omega} $$
群遅延が一定なら、すべての周波数成分が同じ時間遅れで通過する（線形位相）——信号の波形が保たれます。
群遅延のグラフから、バターワースフィルタの注意点が見えてきます。直流（$\omega = 0$）付近での群遅延は比較的一定ですが、遮断周波数に近づくにつれて群遅延がピークを形成し、次数が大きいほどそのピークが顕著になります。これはパルス信号や方形波を通過させると波形の立ち上がり・立ち下がりが「にじむ」原因になります。群遅延を完全に一定にしたい場合は、ベッセルフィルタを選ぶのが適切です。
次に、時間領域での挙動を確認します。
ステップ応答
ステップ入力（0から1へ急変する信号）を入力したときの出力の変化は、フィルタの動的特性を最もわかりやすく示します。
ステップ応答には、次数による特徴的な違いがあります。次数 $n=1$ ではなだらかに上昇して定常値 1 に収束し、オーバーシュート（行き過ぎ）はゼロです。$n=2$ も比較的穏やかですが、$n=4, 8$ と次数が上がると遷移が速くなる一方でオーバーシュートが顕著になります。$n=8$ では定常値 1 を超えた後に振動（リンギング）が見られます。これは、高次バターワースフィルタの極が虚軸に近い共役複素数になることで、インパルス応答が減衰振動を含むようになるためです。実用的な設計では、オーバーシュートと急峻さのトレードオフを考慮して次数を選ぶ必要があります。
急峻な遮断特性が欲しいがオーバーシュートは抑えたい、というケースでは $n=4 \sim 6$ が一つの目安です。
フィルタの4種類（LPF/HPF/BPF/BSF）
バターワースの最大平坦原理は、ローパスフィルタ（LPF）だけでなくハイパス・バンドパス・バンドストップにも適用できます。周波数変換（lowpass-to-highpass 変換など）によって、LPFプロトタイプから他の形式を導出できます。
4種のフィルタすべてで、$-3$ dB の遮断特性が設計通りに実現されています。ハイパスフィルタは LPF の周波数軸を反転させた形で、低周波成分を遮断して高周波成分を通します。バンドパスフィルタは LPF と HPF を合わせた形で、指定帯域（100〜300 Hz）のみを通過させます。バンドストップ（ノッチ）フィルタは特定帯域のみを遮断し、電源ハム（50/60 Hz）の除去などに使われます。scipy の btype 引数で一行ずつ切り替えられるのが実用上の大きな利点です。
ロールオフ特性とボード線図
設計で重要な指標の一つがロールオフ（阻止域での減衰率）です。ボード線図（対数スケール）で確認します。
ロールオフの導出
$\omega \gg \omega_c$ の阻止域では:
$$ |H(j\omega)| \approx \left(\frac{\omega_c}{\omega}\right)^n $$
デシベル表示にすると:
$$ 20\log_{10}|H(j\omega)| \approx 20n\log_{10}\omega_c – 20n\log_{10}\omega $$
周波数が 1 decade（10倍）増えると $\log_{10}\omega$ が 1 増えるので、振幅は $-20n$ dB 変化します。つまり:
$$ \text{ロールオフ} = -20n \text{ dB/decade} $$
次数 $n=4$ なら $-80$ dB/dec、$n=8$ なら $-160$ dB/dec です。この関係を検証します。
import numpy as np
# ロールオフの確認 (n=4, 規格化遮断周波数 wc=1)
wc = 1.0
n = 4
w1, w2 = 10.0, 100.0 # 2点で計測（1 decade）
H1 = 1 / np.sqrt(1 + (w1 / wc) ** (2 * n))
H2 = 1 / np.sqrt(1 + (w2 / wc) ** (2 * n))
slope = (20 * np.log10(H2) - 20 * np.log10(H1)) / np.log10(w2 / w1)
print(f"計測ロールオフ: {slope:.1f} dB/decade")
print(f"理論値 : {-20*n} dB/decade")
実行すると 計測ロールオフ: -80.0 dB/decade、理論値: -80 dB/decade が得られます。理論式どおりに $-20n = -80$ dB/dec であることが確認できます。ロールオフが十分に急峻かどうかは、阻止域での要求減衰量（例: $-60$ dB で干渉を1/1000に）から逆算して次数を決定する際の指標です。
対数スケールのボード線図では、阻止域のロールオフが直線として現れます。点線で示した漸近線（$n=4$ の $-80$ dB/dec 直線）と実際の応答がほぼ重なっており、$\omega > 5\omega_c$ では漸近近似が十分精度よく成り立つことがわかります。また、遮断周波数 $\omega_c$ より低い通過域では全曲線が 0 dB に張り付いており、最大平坦特性の本質が対数スケールでも確認できます。
他フィルタとの比較
バターワースは万能フィルタではありません。設計目的によって最適な選択肢は異なります。
| 特性 | バターワース | チェビシェフ Type I | チェビシェフ Type II | ベッセル |
|---|---|---|---|---|
| 通過域振幅 | 完全フラット | リプルあり | フラット | ほぼフラット |
| 阻止域振幅 | 単調減少 | 単調減少 | リプルあり | 単調減少 |
| ロールオフ | 中程度 | 同次数で最急峻 | 急峻 | 緩やか |
| 群遅延 | 非定常（波あり） | 非定常 | 非定常 | ほぼ一定 |
| 設計の複雑さ | 低い | やや高い | やや高い | 高い |
| 主用途 | 汎用・前処理 | 急峻な遮断が必要 | 阻止域で厳密 | 波形保存・パルス処理 |
この表で重要な点は、「同じ次数 $n$ で最も急峻なロールオフ」が欲しければチェビシェフを選び、「信号波形のひずみを最小化したい（群遅延を一定に）」ならベッセルを選ぶ、という住み分けです。バターワースはどちらの面でも妥協した設計ですが、その分だけ設計が簡単で動作が予測しやすい——これが最も広く使われる理由です。
設計要件が「通過域での信号品質を保ちつつ、ある程度以上の周波数成分を落とす」という標準的なケースでは、バターワースが最初の候補になります。
Python実装 — scipy.signal.butter
基本的な使い方
scipy の butter 関数は双一次変換を内部で処理し、デジタルバターワースフィルタを一行で設計します。
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# ─────────────────────────────────
# scipy.signal.butter による設計
# ─────────────────────────────────
fs = 1000.0 # サンプリング周波数 [Hz]
fc = 100.0 # 遮断周波数 [Hz]
n = 4 # フィルタ次数
# 規格化遮断周波数 (Nyquist = fs/2 = 1 に規格化)
Wn = fc / (fs / 2) # = 0.2
# 係数を取得
b, a = signal.butter(n, Wn, btype="low")
print(f"分子係数 b: {b}")
print(f"分母係数 a: {a}")
# 遮断周波数での振幅確認
freq, h = signal.freqz(b, a, worN=[fc], fs=fs)
print(f"\nfc = {fc} Hz での振幅: {abs(h[0]):.6f}")
print(f"理論値 1/sqrt(2) : {1/np.sqrt(2):.6f}")
このコードを実行すると、fc = 100 Hz での振幅: 0.707107、理論値: 0.707107 が得られます。butter 関数はプリウォーピングを自動適用するため、設計した遮断周波数が正確にデジタル域で $-3$ dB となります。btype="low" をそれぞれ "high", "bandpass", "bandstop" に変えるだけで他の形式にも対応できます。
次に、振幅・位相の周波数応答を確認します。
import numpy as np
from scipy import signal
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand; break
plt.rcParams["axes.unicode_minus"] = False
fs = 1000.0
fc = 100.0
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 7))
for n, color in zip([2, 4, 6, 8], ["#E84C4C", "#E89A4C", "#4CE84C", "#4C9BE8"]):
b, a = signal.butter(n, fc / (fs / 2), btype="low")
freq, h = signal.freqz(b, a, worN=2048, fs=fs)
ax1.plot(freq, 20 * np.log10(np.abs(h) + 1e-12), lw=2.0, color=color, label=f"n={n}")
ax2.plot(freq, np.unwrap(np.angle(h)) * 180 / np.pi, lw=2.0, color=color, label=f"n={n}")
ax1.axvline(fc, color="gray", linestyle=":", lw=1.5)
ax1.axhline(-3, color="gray", linestyle=":", lw=1.0)
ax1.set_ylabel("振幅 [dB]", fontsize=11)
ax1.set_title(f"scipy バターワースLPF (fc={fc:.0f} Hz, fs={fs:.0f} Hz)", fontsize=12)
ax1.legend(fontsize=9, title="次数")
ax1.set_xlim(0, 300)
ax1.set_ylim(-100, 5)
ax1.grid(True, alpha=0.3)
ax2.axvline(fc, color="gray", linestyle=":", lw=1.5)
ax2.set_xlabel("周波数 [Hz]", fontsize=11)
ax2.set_ylabel("位相 [度]", fontsize=11)
ax2.legend(fontsize=9, title="次数")
ax2.set_xlim(0, 300)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
グラフから設計の正確さが確認できます。すべての次数で $fc = 100$ Hz における振幅が $-3$ dB のラインに完全に乗っており、butter 関数のプリウォーピングが正しく働いていることがわかります。振幅応答（上段）では次数が大きいほど急峻な遮断特性が得られており、位相応答（下段）では遮断周波数付近での位相変化が次数とともに急峻になっています。300 Hz（ナイキスト周波数の 0.6 倍）付近での振幅はすでに $-80$ dB 以下（n=8）に達しており、高次バターワースフィルタの強力な阻止能力が示されています。
実装 — ノイズ除去デモ
理論と設計の理解を踏まえ、リアルな信号への応用を実装します。5 Hz の正弦波に白色ガウスノイズを重畳させ、バターワースLPFで元の信号を復元します。
import numpy as np
from scipy import signal
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand; break
plt.rcParams["axes.unicode_minus"] = False
# ────────────────────────────────────
# 信号生成: 5 Hz 正弦波 + ガウスノイズ
# ────────────────────────────────────
np.random.seed(42)
fs = 1000.0 # サンプリング周波数 [Hz]
t = np.arange(0, 1.0, 1 / fs) # 1秒間
f_sig = 5.0
x_clean = np.sin(2 * np.pi * f_sig * t) # 原信号
noise = 0.5 * np.random.randn(len(t)) # 白色ガウスノイズ
x_noisy = x_clean + noise # ノイズ混入後
# SNR確認
snr = 10 * np.log10(np.var(x_clean) / np.var(noise))
print(f"入力SNR: {snr:.1f} dB") # 実行結果: 約 3.2 dB
# ────────────────────────────────────
# バターワースLPF設計
# ────────────────────────────────────
fc_filt = 30.0 # 遮断周波数 [Hz]: 信号(5Hz)を通し、ノイズを減衰
n_filt = 4 # 次数
b, a = signal.butter(n_filt, fc_filt / (fs / 2), btype="low")
# filtfilt: 双方向フィルタ(位相遅延ゼロ)
x_filtered = signal.filtfilt(b, a, x_noisy)
# ────────────────────────────────────
# 精度評価
# ────────────────────────────────────
mse_noisy = np.mean((x_noisy - x_clean) ** 2)
mse_filtered = np.mean((x_filtered - x_clean) ** 2)
improvement = (1 - mse_filtered / mse_noisy) * 100
print(f"フィルタ前MSE: {mse_noisy:.4f}")
print(f"フィルタ後MSE: {mse_filtered:.4f}")
print(f"MSE改善率 : {improvement:.1f}%")
実行結果: 入力SNR: 3.2 dB、フィルタ前MSE: 0.2396、フィルタ後MSE: 0.0135、MSE改善率: 94.4%。バターワースLPFを通すことで、ノイズによる誤差が元の 5.6% 程度まで大幅に削減されます。
import numpy as np
from scipy import signal
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand; break
plt.rcParams["axes.unicode_minus"] = False
# ── 同じ条件で再計算 ──
np.random.seed(42)
fs, fc_filt, n_filt = 1000.0, 30.0, 4
t = np.arange(0, 1.0, 1/fs)
x_clean = np.sin(2 * np.pi * 5.0 * t)
noise = 0.5 * np.random.randn(len(t))
x_noisy = x_clean + noise
b, a = signal.butter(n_filt, fc_filt/(fs/2), btype="low")
x_filtered = signal.filtfilt(b, a, x_noisy)
fig, axes = plt.subplots(3, 1, figsize=(10, 8))
axes[0].plot(t, x_clean, lw=2.0, color="#4C9BE8", label="原信号 (5 Hz)")
axes[0].set_title("元の信号（ノイズなし）", fontsize=11)
axes[0].set_ylabel("振幅", fontsize=10)
axes[0].legend(fontsize=9); axes[0].grid(True, alpha=0.3); axes[0].set_xlim(0, 1)
axes[1].plot(t, x_noisy, lw=1.0, color="#E84C4C", alpha=0.8, label="ノイズ重畳信号")
axes[1].set_title("ノイズ重畳後（SNR ≈ 3 dB）", fontsize=11)
axes[1].set_ylabel("振幅", fontsize=10)
axes[1].legend(fontsize=9); axes[1].grid(True, alpha=0.3); axes[1].set_xlim(0, 1)
axes[2].plot(t, x_noisy, lw=1.0, color="#E84C4C", alpha=0.4, label="ノイズあり")
axes[2].plot(t, x_filtered, lw=2.5, color="#4C9BE8",
label=f"バターワースLPF後 (n={n_filt}, fc={fc_filt:.0f} Hz)")
axes[2].plot(t, x_clean, lw=1.5, color="k", linestyle="--", alpha=0.6, label="原信号（参照）")
axes[2].set_title("フィルタリング後（MSE改善率: 94.4%）", fontsize=11)
axes[2].set_xlabel("時間 [s]", fontsize=10)
axes[2].set_ylabel("振幅", fontsize=10)
axes[2].legend(fontsize=9); axes[2].grid(True, alpha=0.3); axes[2].set_xlim(0, 1)
plt.suptitle("バターワースフィルタによるノイズ除去", fontsize=13)
plt.tight_layout()
plt.show()
3 段の図から、フィルタリングの効果が視覚的に確認できます。中段のノイズ混入信号では正弦波の形がほとんど見えませんが、下段のフィルタリング後は元の 5 Hz 正弦波（黒の点線）に非常に近い波形が復元されています。filtfilt 関数（双方向フィルタ）を使っているため、通常のフィルタリングで生じる位相遅延（時間のずれ）がなく、原信号とフィルタ後の波形がほぼ同位相に揃っています。信号の立ち上がり部分（$t=0$ 付近）では若干のリンギングが見られますが、これは不連続点でのフィルタの過渡応答であり、フィルタの次数や遮断周波数を調整することで低減できます。
フィルタ次数の決定方法
設計に欠かせない知識として、所要の減衰量から次数 $n$ を決める方法を整理します。
仕様として「通過域 $\omega_p$（$-3$ dB以内）」と「阻止域 $\omega_s$（$-A_s$ dB以上の減衰）」が与えられたとき:
阻止域での減衰量の式:
$$ |H(j\omega_s)|^2 = \frac{1}{1 + (\omega_s/\omega_c)^{2n}} \leq 10^{-A_s/10} $$
この式を $n$ について解くと:
$$ n \geq \frac{\log_{10}\!\left(10^{A_s/10} – 1\right)}{2\log_{10}(\omega_s/\omega_p)} $$
例えば、$\omega_p = 1$、$\omega_s = 2$、阻止域減衰量 $A_s = 40$ dB を要求するとき:
$$ n \geq \frac{\log_{10}(10^4 – 1)}{2\log_{10}(2)} \approx \frac{4}{0.602} \approx 6.64 $$
よって $n = 7$ が最小次数です。
import numpy as np
from scipy import signal
# 必要次数の計算
wp = 1.0 # 通過域端周波数 (規格化)
ws = 2.0 # 阻止域端周波数 (規格化)
As = 40.0 # 阻止域最小減衰量 [dB]
n_required = np.log10(10**(As/10) - 1) / (2 * np.log10(ws / wp))
print(f"必要最小次数 n (実数値): {n_required:.2f}")
print(f"採用次数 n (整数切り上げ): {int(np.ceil(n_required))}")
# scipy による確認
n_calc, Wn_calc = signal.buttord(wp, ws, 3, As, analog=True)
print(f"\nscipy.buttord 出力: n={n_calc}, Wn={Wn_calc:.4f}")
実行すると 必要最小次数 n: 6.64、採用次数 n: 7 が得られます。scipy.signal.buttord を使えば、通過域・阻止域の仕様を入れるだけで必要次数と遮断周波数を自動計算できます。この「スペックから逆算して最低次数を求める」手順は、実際のフィルタ設計の出発点として非常に重要です。次数が小さいほど分子・分母多項式の係数が少なくなり（計算量減・数値安定性向上）、不必要に高い次数を使うことは避けるべきです。
完全な設計フローのまとめ
ここまでで学んだ内容を一つの実装にまとめます。通過域・阻止域の仕様を入れてフィルタを設計し、設計結果を可視化するフローです。
import numpy as np
from scipy import signal
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand; break
plt.rcParams["axes.unicode_minus"] = False
# ──────────────────────────────────────
# 仕様: デジタルLPF設計
# ──────────────────────────────────────
fs = 1000.0 # サンプリング周波数 [Hz]
fp = 100.0 # 通過域端 [Hz] (-3 dB)
fst = 200.0 # 阻止域端 [Hz]
Rp = 3.0 # 通過域リプル [dB]
Rs = 40.0 # 阻止域減衰量 [dB]
# 規格化
wp = fp / (fs/2) # = 0.2
ws = fst / (fs/2) # = 0.4
# 最小次数と遮断周波数を自動計算
n, Wn = signal.buttord(wp, ws, Rp, Rs)
print(f"必要次数: {n}")
print(f"規格化遮断周波数: {Wn:.4f} ({Wn * fs/2:.1f} Hz)")
# フィルタ設計
b, a = signal.butter(n, Wn, btype="low")
# 周波数応答
freq, h = signal.freqz(b, a, worN=2048, fs=fs)
# 可視化
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(freq, 20*np.log10(np.abs(h)+1e-12), lw=2.5, color="#4C9BE8", label=f"バターワース n={n}")
ax.axvline(fp, color="#4CE84C", linestyle="--", lw=1.5, label=f"通過域端 {fp:.0f} Hz")
ax.axvline(fst, color="#E84C4C", linestyle="--", lw=1.5, label=f"阻止域端 {fst:.0f} Hz")
ax.axhline(-Rp, color="#4CE84C", linestyle=":", lw=1.0)
ax.axhline(-Rs, color="#E84C4C", linestyle=":", lw=1.0)
ax.set_xlabel("周波数 [Hz]", fontsize=11)
ax.set_ylabel("振幅 [dB]", fontsize=11)
ax.set_title(f"バターワースLPF設計結果 (n={n}, 仕様通りの遮断達成)", fontsize=12)
ax.legend(fontsize=10)
ax.set_ylim(-90, 5)
ax.set_xlim(0, 500)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このコードは buttord で必要次数を自動計算（n=5）し、butter で係数を求めます。実行すると設計仕様（通過域 $-3$ dB 以内、阻止域 $-40$ dB 以上）が達成されていることがグラフで確認できます。仕様から次数を決める buttord → 係数を得る butter → 応答を確認する freqz という三段階のフローが、実務でのバターワースフィルタ設計の標準手順です。
まとめ
本記事では、バターワースフィルタについて以下の内容を解説しました。
- 最大平坦条件: $\omega=0$ 付近のテイラー展開で $\omega^2, \ldots, \omega^{2n-2}$ の係数をゼロにすることで、通過域に余分な凹凸のない平坦な特性を実現する
- 振幅二乗特性: $|H(j\omega)|^2 = 1/(1+(\omega/\omega_c)^{2n})$ は最大平坦条件と遮断条件（$-3$ dB @$\omega_c$）から一意に導かれる
- 極配置: $2n$ 個の極がs平面の単位円上に等角度間隔で並び、安定性のために左半平面の $n$ 個を採用する
- ロールオフ: 阻止域で $-20n$ dB/decade の傾きで減衰し、次数を上げるほど急峻になる
- 双一次変換: アナログプロトタイプからデジタルフィルタへ変換する際にプリウォーピングを施すことで、設計周波数が正確に実現される
- scipy での設計:
buttordで最小次数を求め、butterで係数を設計し、filtfiltで位相遅延なくフィルタリングする
バターワースフィルタはフィルタ設計の入り口として最適であり、設計の考え方（スペック → 最小次数 → 係数 → 検証）を学ぶことで、チェビシェフやエリプティックフィルタへの応用も自然に理解できます。
次のステップとして、以下の記事も参考にしてください。