ローパスフィルタの設計と比較：移動平均・バターワース・チェビシェフ

scipy.signal.butter / scipy.signal.cheby1 / scipy.signal.filtfilt / scipy.signal.freqz でローパスフィルタをPython実装し3種比較。移動平均・バターワース・チェビシェフI型の周波数応答・群遅延・ステップ応答・実装コードと、カットオフ周波数の決め方、ノイズ除去・センサ信号平滑化への応用までまとめます。

22 min read

はじめに

ローパスフィルタ（Low-Pass Filter; LPF）は、信号処理において最も基本的かつ重要なフィルタです。カットオフ周波数以下の低周波成分を通過させ、高周波ノイズを除去する役割を持ちます。

本記事では、代表的な3種類のローパスフィルタを取り上げます。

  1. 移動平均フィルタ（FIR型）：最も単純な構造で線形位相を持つ
  2. バターワースフィルタ（IIR型）：通過域で最大限に平坦な振幅特性を持つ
  3. チェビシェフI型フィルタ（IIR型）：通過域にリプルを許容する代わりに急峻な遷移帯域を実現する

それぞれの理論的背景を示したうえで、Pythonによる設計と特性比較を行います。

本記事の位置づけ：本記事は3種のフィルタを同一の次数・カットオフ条件下で横並びに比較する入門ハブです。各フィルタの詳細な導出・エッジケース・実装上の注意点は専用記事（ 移動平均フィルタの種類と比較 、 バターワースフィルタの設計原理とPython実装 、 チェビシェフフィルタの設計原理とPython実装 ）に譲り、本記事では「同じ次数・同じカットオフを指定したときに実際は何が起こるのか」を数値とグラフで検証することに焦点を当てます。

移動平均フィルタ（Moving Average Filter）

移動平均フィルタは、直近 \(N\) 個のサンプルの算術平均を出力するFIRフィルタです。

伝達関数

入力信号 \(x[n]\) に対して、出力は次のように定義されます。

\[y[n] = \frac{1}{N}\sum_{k=0}^{N-1} x[n-k] \tag{1}\]

Z変換を適用すると、伝達関数は次のようになります。

\[H(z) = \frac{1}{N}\sum_{k=0}^{N-1} z^{-k} = \frac{1}{N} \cdot \frac{1 - z^{-N}}{1 - z^{-1}} \tag{2}\]

周波数応答

\(z = e^{j\omega}\) を代入すると、周波数応答が得られます。

\[H(e^{j\omega}) = \frac{1}{N} \cdot \frac{1 - e^{-jN\omega}}{1 - e^{-j\omega}} = \frac{1}{N} \cdot \frac{\sin(N\omega/2)}{\sin(\omega/2)} \cdot e^{-j(N-1)\omega/2} \tag{3}\]

振幅特性は \(\frac{1}{N}\left|\frac{\sin(N\omega/2)}{\sin(\omega/2)}\right|\) であり、sinc関数に似た形状を示します。この特性により、移動平均フィルタは阻止域での減衰が緩やかで、周波数選択性は高くありません。一方、位相特性は \(-(N-1)\omega/2\) であり、完全な線形位相を持つという大きな利点があります。

バターワースフィルタ（Butterworth Filter）

バターワースフィルタは、通過域において最も平坦な振幅特性（最大平坦特性）を持つIIRフィルタです。

振幅特性

\(N\) 次バターワースフィルタの振幅の2乗特性は次のように定義されます。

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

ここで、\(\omega_c\) はカットオフ周波数、\(N\) はフィルタの次数です。

この式から、以下の性質がわかります。

  • \(\omega = 0\) のとき \(|H| = 1\) （DCゲインは1）
  • \(\omega = \omega_c\) のとき \(|H| = 1/\sqrt{2}\) （\(-3\) dB）
  • 次数 \(N\) が大きいほど遷移帯域が急峻になる

バターワースフィルタは通過域・阻止域ともにリプルがなく、滑らかな特性を持つため、汎用的なフィルタとして広く使われています。

振幅特性の導出：なぜこの形に一意に定まるのか

式 \((4)\) は天下りに与えられることが多いですが、「最大平坦性」という要求から導出できます。実係数フィルタの振幅二乗特性 \(|H(j\omega)|^2\) は \(\omega \to -\omega\) について偶関数なので、\(\Omega = \omega/\omega_c\) を使って

\[ |H(j\omega)|^2 = \frac{1}{1 + f(\Omega)}, \qquad f(\Omega) = \sum_{m \ge 1} c_m \Omega^{2m} \]

と展開できます（奇数次の項は偶関数の条件から現れません）。「\(\Omega = 0\) の周りで最初の \(2N-1\) 階までの導関数がすべて \(0\) 」という最大平坦性の条件は、\(f\) の中に次数 \(2, 4, \ldots, 2(N-1)\) の項があるとその階数の導関数が非ゼロになってしまうため、\(c_1 = c_2 = \cdots = c_{N-1} = 0\) を要求します。残るのは最低次の項 \(f(\Omega) = c_N \Omega^{2N}\) だけであり、これが「最大平坦」を満たす唯一の（正規化を除いた）形です。最後に \(\Omega = 1\) （\(\omega = \omega_c\) ）で \(-3\,\text{dB}\) 、すなわち \(f(1) = 1\) となるよう正規化すると \(c_N = 1\) が決まり、式 \((4)\) が一意に導かれます。

極の配置（s平面）

式 \((4)\) の解析接続 \(|H(s)|^2 = H(s)H(-s)\) を考えると、極は次の方程式の解として求まります。

\[ 1 + \left(\frac{s}{j\omega_c}\right)^{2N} = 0 \quad \Longleftrightarrow \quad \left(\frac{s}{j\omega_c}\right)^{2N} = e^{j\pi(2k-1)}, \quad k \in \mathbb{Z} \]

これを \(s\) について解くと、次数 \(N\) のバターワースフィルタの極は次のように与えられます。

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

安定な（左半平面の）極だけを選ぶと、\(N\) 個の極は半径 \(\omega_c\) の円（バターワース円）上に等間隔に並びます。この「極が円周上に等間隔」という構造が、通過域の平坦性とほどよい遷移帯域の急峻さを両立させる幾何学的な理由です（\(s\) 平面の極をデジタルフィルタの\(z\) 平面に写像する双一次変換の詳細は バターワースフィルタの設計原理とPython実装 を参照）。式 \((7)\) は後述のPythonによる極配置の実測比較（scipy.signal.tf2zpk）で数値的に検証します。

設計

Pythonでは scipy.signal.butter を使って設計できます。

from scipy.signal import butter
b, a = butter(N=4, Wn=0.3)  # 4次、正規化カットオフ0.3

チェビシェフI型フィルタ（Chebyshev Type I Filter）

チェビシェフI型フィルタは、通過域にチェビシェフ多項式に基づく等リプル特性を持つIIRフィルタです。

振幅特性

\(N\) 次チェビシェフI型フィルタの振幅の2乗特性は次のとおりです。

\[|H(j\omega)|^2 = \frac{1}{1 + \varepsilon^2 T_N^2\left(\frac{\omega}{\omega_c}\right)} \tag{5}\]

ここで、\(\varepsilon\) はリプルの大きさを制御するパラメータ、\(T_N\) は \(N\) 次チェビシェフ多項式です。チェビシェフ多項式は次の漸化式で定義されます。

\[T_0(x) = 1, \quad T_1(x) = x, \quad T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x) \tag{6}\]

通過域（\(\omega \le \omega_c\) ）では \(T_N\) が \([-1, 1]\) の範囲で振動するため、振幅にリプルが生じます。その代わり、同じ次数のバターワースフィルタと比較して遷移帯域がより急峻になるという利点があります。

なぜ「等リプル」がチェビシェフ多項式から生まれるのか

\(x = \cos\theta\) （\(|x|\le 1\) ）と置くと、式 \((6)\) の漸化式は三角関数の倍角公式と一致し、\(T_N(\cos\theta) = \cos(N\theta)\) が成り立ちます。\(\cos(N\theta)\) は \(\theta\) が実数の範囲で常に \([-1, 1]\) の間を振動するため、通過域では \(T_N^2\) が \(0\) から \(1\) の間を繰り返し往復し、式 \((5)\) の \(|H(j\omega)|^2\) も \(1/(1+\varepsilon^2)\) から \(1\) の間で振動します。これが等リプルの直接の原因です。さらに、チェビシェフ多項式は「区間 \([-1,1]\) 上で最高次係数を固定したモニック多項式の中で、絶対値の最大値（supノルム）を最小化する」というミニマックス最適性（equioscillation theorem）を満たす一意な多項式であることが知られています。つまり、リプルの振幅を式 \((5)\) の \(\varepsilon\) で指定したとき、遷移帯域の急峻さを最大化する多項式は数学的にチェビシェフ多項式以外あり得ません（証明は チェビシェフフィルタの設計原理とPython実装 を参照）。

極の配置（s平面）：バターワースとの違い

チェビシェフI型フィルタの極は、バターワースのような円ではなく楕円上に並びます。\(\varepsilon\) をリプルパラメータとして

\[ \phi = \frac{1}{N}\sinh^{-1}\left(\frac{1}{\varepsilon}\right), \qquad \theta_k = \frac{\pi(2k-1)}{2N} \]

と定義すると、\(k\) 番目の極は

\[ s_k = -\omega_c \sinh(\phi)\sin(\theta_k) + j\,\omega_c \cosh(\phi)\cos(\theta_k), \quad k = 1, \ldots, N \tag{8} \]

で与えられます（実軸方向の半径が \(\omega_c \sinh\phi\) 、虚軸方向の半径が \(\omega_c \cosh\phi\) の楕円）。\(\sinh\phi < \cosh\phi\) なので、この楕円はバターワース円よりも虚軸（\(j\omega\) 軸）方向に扁平に潰れた形になり、極が単位円（\(z\) 平面では単位円、\(s\) 平面では虚軸）に近づきます。極が虚軸に近いほど遷移帯域は急峻になりますが、同時に減衰比（ダンピング）が下がるため、ステップ応答のオーバーシュートや群遅延のピークが大きくなります。この関係は後述のPythonによる実測比較で数値的に確認します。

エッジケース：偶数次のDC利得は0dBにならない

式 \((6)\) の漸化式から \(T_N(0)\) を計算すると、奇数次では \(T_N(0)=0\) 、偶数次では \(T_N(0) = (-1)^{N/2} = \pm 1\) となります。これを式 \((5)\) に \(\omega=0\) を代入した式へ適用すると、

\[ |H(j0)|^2 = \frac{1}{1 + \varepsilon^2 T_N^2(0)} = \begin{cases} 1 & (N \text{ が奇数}) \\ \dfrac{1}{1+\varepsilon^2} & (N \text{ が偶数}) \end{cases} \]

となり、偶数次のチェビシェフI型フィルタはDC（\(\omega=0\) ）でも \(0\,\text{dB}\) にならず、\(-R_p\,\text{dB}\) だけ減衰した状態からスタートするという見落としやすい性質があります（\(R_p\) は通過域リプル \([\text{dB}]\) ）。本記事の比較例では \(N=4\) （偶数）・\(R_p=1\,\text{dB}\) を使うため、この効果が後述のステップ応答の実測値に現れます。

設計

Pythonでは scipy.signal.cheby1 を使って設計できます。rp パラメータで通過域リプルをdB単位で指定します。

from scipy.signal import cheby1
b, a = cheby1(N=4, rp=1, Wn=0.3)  # 4次、リプル1dB、正規化カットオフ0.3

Pythonによる比較

以下のコードでは、3種類のフィルタを同一条件で設計し、周波数応答・群遅延・ステップ応答を比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import (
    butter, cheby1, freqz, group_delay, dstep, dlti, tf2zpk
)

# --- フィルタ設計 ---
N_ORDER = 4        # フィルタ次数
WN = 0.3           # 正規化カットオフ周波数 (ナイキスト周波数に対する比)
N_MA = 13           # 移動平均のタップ数 (奇数を推奨)
RP = 1.0            # チェビシェフI型の通過域リプル [dB]

# 移動平均フィルタ
b_ma = np.ones(N_MA) / N_MA
# 注意: dlti/dstep は len(a) >= len(b) を要求する（新しいSciPyでは a=[1.0] のままだと
# "Improper transfer function" で例外になる）。分母を b_ma と同じ長さまで0埋めする。
a_ma = np.concatenate(([1.0], np.zeros(N_MA - 1)))

# バターワースフィルタ
b_bw, a_bw = butter(N_ORDER, WN)

# チェビシェフI型フィルタ
b_cb, a_cb = cheby1(N_ORDER, RP, WN)

filters = [
    ("Moving Average (N=13)", b_ma, a_ma),
    ("Butterworth (N=4)", b_bw, a_bw),
    ("Chebyshev I (N=4, rp=1dB)", b_cb, a_cb),
]

# --- 1. 周波数応答 ---
fig, ax = plt.subplots(figsize=(8, 5))
for label, b, a in filters:
    w, h = freqz(b, a, worN=2048)
    freq = w / np.pi  # 正規化周波数
    mag_db = 20 * np.log10(np.abs(h) + 1e-12)
    ax.plot(freq, mag_db, label=label)

ax.set_xlabel("Normalized Frequency (×π rad/sample)")
ax.set_ylabel("Magnitude (dB)")
ax.set_title("Frequency Response Comparison")
ax.set_xlim(0, 1)
ax.set_ylim(-60, 5)
ax.axvline(WN, color="gray", linestyle="--", alpha=0.5, label=f"Cutoff = {WN}")
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

# --- 2. 群遅延 ---
fig, ax = plt.subplots(figsize=(8, 5))
for label, b, a in filters:
    w, gd = group_delay((b, a), w=2048)
    freq = w / np.pi
    ax.plot(freq, gd, label=label)

ax.set_xlabel("Normalized Frequency (×π rad/sample)")
ax.set_ylabel("Group Delay (samples)")
ax.set_title("Group Delay Comparison")
ax.set_xlim(0, 1)
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

# --- 3. ステップ応答 ---
fig, ax = plt.subplots(figsize=(8, 5))
n_steps = 60
for label, b, a in filters:
    system = dlti(b, a, dt=1)
    t, y = dstep(system, n=n_steps)
    # 注意: dstep が返す t は既に1次元配列 (n,) なので t[0] は先頭の1要素だけになってしまう
    ax.step(t, y[0].flatten(), label=label, where="post")

ax.set_xlabel("Sample")
ax.set_ylabel("Amplitude")
ax.set_title("Step Response Comparison")
ax.axhline(1.0, color="gray", linestyle="--", alpha=0.5)
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

# --- 4. 極零配置（z平面） ---
fig, axes = plt.subplots(1, 3, figsize=(12, 4.3))
theta = np.linspace(0, 2 * np.pi, 400)
for ax, (label, b, a) in zip(axes, filters):
    z, p, k = tf2zpk(b, a)
    ax.plot(np.cos(theta), np.sin(theta), color="gray", linestyle="--", linewidth=1, alpha=0.6)
    ax.scatter(z.real, z.imag, s=70, facecolors="none", edgecolors="C0", linewidths=1.8, label="zeros")
    ax.scatter(p.real, p.imag, s=70, marker="x", color="C1", linewidths=1.8, label="poles")
    ax.set_xlim(-1.4, 1.4)
    ax.set_ylim(-1.4, 1.4)
    ax.set_aspect("equal")
    ax.set_title(label, fontsize=10)
    ax.set_xlabel("Re(z)")
    ax.grid(True, alpha=0.3)
    ax.legend(loc="lower left", fontsize=7)
    max_p = np.max(np.abs(p)) if len(p) else 0.0
    ax.text(0.02, 0.98, f"max|pole|={max_p:.3f}", transform=ax.transAxes, fontsize=8, va="top")
axes[0].set_ylabel("Im(z)")
fig.suptitle("Pole-Zero Plot Comparison (z-plane)")
plt.tight_layout()
plt.show()

周波数応答の比較

群遅延の比較

ステップ応答の比較

極零配置の比較

上記のコードを実際に実行し、freqz / group_delay / dstep / tf2zpk の出力から読み取った実測値は次のとおりです（SciPy 1.18.0、NumPy 2.4.2で実行）。

周波数応答の読み取り

  • 移動平均フィルタ（\(N=13\) ）は、公称カットオフ \(\text{Wn}=0.3\) の時点ですでに \(-31.6\,\text{dB}\) まで減衰しているものの、阻止域の第一サイドローブ（\(0.2205\pi\) 付近）でのピークは \(-13.09\,\text{dB}\) にとどまり、以降のサイドローブも \(-17\) 〜\(-22\,\text{dB}\) 程度でしか減衰しない。「阻止域の減衰が浅い」という定性的な記述は、この \(-13\,\text{dB}\) という実測値で裏付けられる
  • バターワースフィルタ（\(N=4\) ）は \(\text{Wn}=0.3\) ちょうどで \(-3.014\,\text{dB}\) 、実測の \(-3\,\text{dB}\) クロス周波数は \(0.2999\pi\) と、butter の Wn パラメータが文字通り \(-3\,\text{dB}\) 点を指すことが確認できる
  • チェビシェフI型フィルタ（\(N=4\) 、\(R_p=1\,\text{dB}\) ）は \(\text{Wn}=0.3\) で \(-1.005\,\text{dB}\) （\(\approx -R_p\) ）にすぎず、実測の \(-3\,\text{dB}\) クロス周波数は \(0.3135\pi\) とバターワースより広い。つまり同じ Wn=0.3 を指定しても、cheby1 の Wn は \(-3\,\text{dB}\) 点ではなく通過域リプルの端を指しており、実際の \(-3\,\text{dB}\) 帯域幅はバターワースより約4.5%広い（詳細は後述のエッジケースを参照）。周波数 \(0.5\pi\) における減衰量で比較すると、バターワースが \(-23.45\,\text{dB}\) に対しチェビシェフI型は \(-33.11\,\text{dB}\) と、同次数でも約\(10\,\text{dB}\) 深く減衰しており、遷移帯域の急峻さの違いが定量的に確認できる

群遅延の読み取り

  • 移動平均フィルタは理論値どおり、低周波（\(0.01\pi\) ）・カットオフ（\(0.3\pi\) ）・\(0.5\pi\) のいずれでも群遅延が \(6.000\) サンプル（\(=(N-1)/2=(13-1)/2\) ）で完全に一定。線形位相であることが数値的に確認できる
  • バターワースフィルタの群遅延は低周波で \(2.566\) サンプル、カットオフ付近（\(0.3\pi\) ）で \(4.567\) サンプルとピークを迎え、\(0.5\pi\) では \(1.528\) サンプルまで下がる
  • チェビシェフI型フィルタの群遅延はカットオフ付近（\(0.3\pi\) ）で \(9.870\) サンプルに達し、同じ条件のバターワースの \(4.567\) サンプルの2倍以上。等リプルによる急峻な遷移帯域と引き換えに、カットオフ近傍での位相非線形性がバターワースよりも顕著に大きくなることが実測値からわかる

ステップ応答の読み取り

  • 移動平均フィルタはオーバーシュートなし（\(0.00\%\) ）で定常値 \(1.0000\) に到達し、2%整定時間は12サンプル
  • バターワースフィルタは定常値 \(1.0000\) に対しピーク \(1.1322\) 、オーバーシュート \(13.22\%\) 、2%整定時間は11サンプル
  • チェビシェフI型フィルタはピーク \(1.1065\) に達するが、定常値そのものが \(1.0000\) ではなく \(0.8915\) （\(\approx 10^{-1/20}=0.8913\) 、偶数次のDC利得が\(-R_p\,\text{dB}\) になるという前述のエッジケースの帰結）。この定常値を基準にすると見かけ上のオーバーシュートは \(24.12\%\) に膨らみ、2%整定時間も21サンプルと3種の中で最も長い

極配置から見た特性の違い

上図の極零プロット（tf2zpk の実測値）は、これまでの数値の背景にある幾何学的な理由を示しています。

  • 移動平均フィルタ：分母が定数（\(a=[1,0,\ldots,0]\) ）のFIRフィルタなので、12個の極はすべて原点 \(z=0\) （max|pole|=0.000）に重なる。極が原点にしかないことは「群遅延が周波数に依存しない＝線形位相」であることの必要条件そのもの
  • バターワースフィルタ：4個の極は原点からの最大距離（最大極半径）\(0.726\) の位置にあり、式\((7)\) の円配置がz平面（双一次変換後）でも概ね保たれる
  • チェビシェフI型フィルタ：4個の極の最大半径は \(0.893\) とバターワースより単位円に近い。式\((8)\) で導いた楕円配置により極が単位円に近づくほど、周波数応答の遷移は急峻になる代わりに減衰比が下がり、ステップ応答のオーバーシュート・整定時間の悪化、カットオフ付近の群遅延の増大という形で現れる——これはまさに上記で実測した数値どおりの結果です

フィルタ選択ガイド

特性移動平均バターワースチェビシェフI型
通過域平坦最大平坦等リプル
遷移帯域緩やか中程度急峻
阻止域減衰浅い中程度深い
位相特性線形非線形非線形
群遅延一定周波数依存周波数依存
計算コスト低い中程度中程度
主な用途ノイズ平滑化汎用フィルタリング急峻な遮断が必要な場合

用途に応じた選び方の目安は次のとおりです。

  • 位相歪みを避けたい場合（波形の形状が重要）→ 移動平均フィルタ
  • 通過域の平坦性が重要な場合（計測信号の忠実な再現）→ バターワースフィルタ
  • 遷移帯域をできるだけ狭くしたい場合（隣接する周波数成分の分離）→ チェビシェフI型フィルタ

ゼロ位相フィルタリング（filtfilt）による群遅延の除去

バターワースやチェビシェフI型のようなIIRフィルタは非線形位相を持ちますが、オフライン処理（記録済みデータの後処理）であれば scipy.signal.filtfilt で信号を順方向と逆方向の両方から通すことで、見かけ上の群遅延をゼロにできます。順方向・逆方向のフィルタ処理は同一フィルタを2回通すことと等価なので、周波数応答は2乗されます。

import numpy as np
from scipy.signal import butter, cheby1, freqz, filtfilt

b_bw, a_bw = butter(4, 0.3)
b_cb, a_cb = cheby1(4, 1.0, 0.3)

w, h_bw = freqz(b_bw, a_bw, worN=8192)
w, h_cb = freqz(b_cb, a_cb, worN=8192)
freq = w / np.pi
idx_wn = np.argmin(np.abs(freq - 0.3))

# 単発フィルタの利得 |H|^2 [dB] = filtfilt 適用後の実効利得
mag_bw_single = 20 * np.log10(np.abs(h_bw[idx_wn]))
mag_bw_filtfilt = 20 * np.log10(np.abs(h_bw[idx_wn]) ** 2)
mag_cb_single = 20 * np.log10(np.abs(h_cb[idx_wn]))
mag_cb_filtfilt = 20 * np.log10(np.abs(h_cb[idx_wn]) ** 2)
print(f"Butterworth: 単発 {mag_bw_single:.3f} dB → filtfilt後 {mag_bw_filtfilt:.3f} dB")
print(f"Chebyshev I: 単発 {mag_cb_single:.3f} dB → filtfilt後 {mag_cb_filtfilt:.3f} dB")

# 群遅延がゼロになることをインパルス応答のピーク位置で確認
x = np.zeros(200)
x[100] = 1.0
y = filtfilt(b_bw, a_bw, x)
print("filtfilt後のインパルス応答ピーク位置:", np.argmax(np.abs(y)), "（入力ピーク位置100と一致すれば遅延ゼロ）")

実行結果:

Butterworth: 単発 -3.014 dB → filtfilt後 -6.027 dB
Chebyshev I: 単発 -1.005 dB → filtfilt後 -2.011 dB
filtfilt後のインパルス応答ピーク位置: 100 （入力ピーク位置100と一致すれば遅延ゼロ）

filtfilt はフィルタを2回通すため、同じ Wn を指定しても実効的な減衰量はちょうど2倍（dB表示で）になる——バターワースはカットオフでの減衰が \(-3.014\,\text{dB}\) から \(-6.027\,\text{dB}\) へ、チェビシェフI型の通過域リプルも \(-1.005\,\text{dB}\) から \(-2.011\,\text{dB}\) へと拡大します。つまり filtfilt で使うフィルタは、単発で使う場合よりもやや緩め（次数を落とすかカットオフを高めに取る）に設計しないと、意図より狭い実効帯域になってしまいます。一方でインパルス応答のピーク位置は入力と同じ100番目のサンプルに一致しており、位相遅延が完全にゼロになっていることも確認できます。

実務上のエッジケース・注意点

これまでの導出・実測から明らかになった、実装時に見落としやすい注意点をまとめます。

  1. Wn パラメータの意味はフィルタ関数ごとに異なる：scipy.signal.butter の Wn は \(-3\,\text{dB}\) 点を指しますが、scipy.signal.cheby1 の Wn は通過域リプルの端（\(-R_p\,\text{dB}\) 点）を指します。前述のとおり、同じ Wn=0.3 でも実測の \(-3\,\text{dB}\) 点はバターワースが \(0.2999\pi\) 、チェビシェフI型が \(0.3135\pi\) とずれるため、「次数とカットオフを揃えたから同じ帯域幅のはず」という前提で比較すると誤った結論を導きます

  2. 移動平均フィルタの阻止域減衰は本質的に浅い：第一サイドローブでも \(-13.09\,\text{dB}\) にしかならないため、単一の妨害周波数を狙い撃ちで除去したい場合はサイドローブの零点にちょうど合わせる（\(N\) の調整）か、 サビツキー・ゴーレイフィルタ やIIRフィルタへの切り替えを検討すべきです

  3. 偶数次チェビシェフI型のDC利得は \(0\,\text{dB}\) ではない：本記事の実測（\(N=4\) 、\(R_p=1\,\text{dB}\) ）では定常値が \(0.8915\) にしかならず、直流成分の忠実な伝達が必要な用途では見込み違いの原因になります

  4. 高次IIRフィルタは伝達関数形式（b, a）のままだと数値的に不安定になりうる：\(N=4\) 程度では問題になりませんが、次数を上げると危険です。実際に butter(N, 0.3) を次数 \(N=20,30,40,50\) で設計し、本来完全に平坦なはずの通過域（最初の400ビン）でのゲインのばらつきを b,a 形式と sos（2次セクション）形式で比較すると次のようになります。

    N=20: max|a|=3.258e+02  b,a形式の通過域ばらつき=2.103e-10 dB   sos形式=2.604e-14 dB
N=30: max|a|=1.186e+04  b,a形式の通過域ばらつき=3.248e-07 dB   sos形式=3.857e-14 dB
N=40: max|a|=4.462e+05  b,a形式の通過域ばらつき=1.295e-04 dB   sos形式=4.436e-14 dB
N=50: max|a|=1.745e+07  b,a形式の通過域ばらつき=2.441e-01 dB   sos形式=4.436e-14 dB

    \(N=50\) では分母係数の絶対値が最大 \(1.7\times10^7\) にまで達し、b,a 形式の通過域ゲインが本来 \(0\,\text{dB}\) で平坦であるべきところ最大 \(0.24\,\text{dB}\) もばらつきます。sos 形式は同じ条件で常に \(10^{-14}\,\text{dB}\) 未満の誤差に収まるため、次数が2桁に達するような設計では output='sos' を使うべきです（本記事の \(N=4\) の比較では実害はありません）

  5. scipy.signal.dlti / dstep にまつわる2つの罠：(a) dlti / dstep は len(a) >= len(b) を要求するため、移動平均フィルタのような純粋なFIR（\(a=[1.0]\) ）をそのまま渡すと、新しいSciPy（本記事の検証環境では1.18.0）では ValueError: Improper transfer function になります。分母を分子と同じ長さまで0埋めしてから渡す必要があります。(b) dstep が返す時刻配列 t はすでに1次元配列 (n,) であり、t[0] はその先頭の1要素（スカラー）にすぎません。単入力単出力系のプロットでは t[0].flatten() ではなく t をそのまま使う必要があります（y の方は出力ごとのタプルなので y[0].flatten() で正しい）。この2点はいずれも本記事のコードで対応済みです

近年の研究動向：ヒューリスティック・FIR・IIRの実データ横断比較

本記事のような「異なるフィルタ族を同一条件で比較する」というアプローチは、実データを使った検証研究でも行われています。Raju, Friedman, Bouman, Komogortsev（2023）は、視線追跡装置EyeLink 1000から得られる眼球運動データに対し、ヒューリスティックな移動平均系フィルタ（STD/EXTRA）、サビツキー・ゴーレイフィルタ、IIR（バターワース）フィルタ、FIRフィルタという異なる系統のフィルタを横並びに適用し、ノイズ除去性能を定量比較しました（Journal of Eye Movement Research, 2023）。

この研究では、\(-30\,\text{dB}\) 減衰に達する周波数がFIRフィルタで\(93\,\text{Hz}\) 、IIR（バターワース）フィルタで\(102\,\text{Hz}\) 、サビツキー・ゴーレイフィルタでは\(204\,\text{Hz}\) という結果が報告されており、遷移帯域の鋭さでFIR・IIRがヒューリスティック系フィルタを大きく上回ることが示されています。また速度信号のノイズ（標準偏差）は、フィルタ適用前の\(120.19\,\text{deg/s}\) からFIR・IIRの両方で約\(23\,\text{deg/s}\) まで低減されました。さらに、FIRとIIR（バターワース）はいずれも本記事の「ゼロ位相フィルタリング」の節で扱ったfiltfilt相当の双方向フィルタリングとして実装され、位相遅延を除去している点も本記事の内容と直接対応します。一方で、全てのデジタルフィルタが信号のラグ1自己相関を未フィルタ時の\(0.58\) から約\(0.97\) まで押し上げるという副作用も報告されており、ノイズ除去と引き換えに時系列の統計的性質（自己相関構造）が変化する点は、フィルタ選択時に見落とされがちな注意点として参考になります。

この結果は、本記事で比較した「移動平均（FIRの特殊形）・バターワース（IIR）・チェビシェフI型（IIR）」という枠組みが、理論だけでなく実データに基づく検証でも同様の傾向（IIR・FIR系はヒューリスティックな平滑化より遷移帯域が鋭く、ノイズ抑制性能も高い）を示すことを裏付けています。

まとめ

本記事では、移動平均フィルタ、バターワースフィルタ、チェビシェフI型フィルタの3種のローパスフィルタについて、振幅特性の導出（最大平坦性・等リプル性）と極配置を示したうえで、同一の次数・カットオフ条件下でPythonにより周波数応答・群遅延・ステップ応答・極零配置を実測比較しました。

実測でわかった要点は次のとおりです。

  • 同じ Wn を指定しても、butter は \(-3\,\text{dB}\) 点、cheby1 は通過域リプル端を意味するため、実際の \(-3\,\text{dB}\) 帯域幅はチェビシェフI型の方が約4.5%広くなる（実測: バターワース\(0.2999\pi\) 、チェビシェフI型\(0.3135\pi\) ）
  • 移動平均フィルタの阻止域第一サイドローブは\(-13.09\,\text{dB}\) にとどまり、阻止域減衰は本質的に浅い
  • チェビシェフI型（偶数次）のステップ応答は定常値が\(1.0\) ではなく\(0.8915\) になり、見かけ上のオーバーシュートが\(24.12\%\) に達する
  • 極零プロットで見ると、チェビシェフI型の最大極半径（\(0.893\) ）はバターワース（\(0.726\) ）より単位円に近く、これがカットオフ付近の群遅延（\(9.870\) vs \(4.567\) サンプル）やオーバーシュートの差を生む幾何学的な理由になっている

各フィルタにはそれぞれ長所と短所があり、アプリケーションの要件に応じて適切なフィルタを選択することが重要です。また、EMAフィルタ（1次IIRフィルタ）との関連については 指数移動平均（EMA）フィルタの周波数特性 も併せてご参照ください。状態空間モデルに基づくフィルタリング手法（カルマンフィルタ、拡張カルマンフィルタ、UKF、粒子フィルタ）については 信号処理におけるフィルタリング手法の基礎 で体系的に解説しています。

参考

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

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

関連記事