サビツキー・ゴーレイフィルタとは
サビツキー・ゴーレイ（Savitzky-Golay）フィルタは、移動窓内で局所的に多項式を最小二乗法でフィッティングし、その多項式の値を平滑化結果として用いるフィルタです。1964年にAbraham SavitzkyとMarcel J.E. Golayによって提案されました。
移動平均フィルタやEMAフィルタが信号のピークや急峻な変化を鈍らせてしまうのに対し、サビツキー・ゴーレイフィルタは信号の形状（ピーク位置、幅、高さ）を保存しながらノイズを除去できるという大きな特徴を持ちます。
主な用途
- 分光分析データの平滑化
- クロマトグラフィーのピーク検出
- 生体信号（ECG、EEG）の前処理
- 信号の微分の推定
数学的原理
局所多項式フィッティング
窓幅 \(2m+1\)（\(m\) は片側のデータ点数）の区間内で、\(p\) 次多項式を最小二乗法でフィッティングします。
窓の中心を \(n=0\) として、データ点 \(x_{-m}, x_{-m+1}, \dots, x_m\) に対して、以下の多項式を当てはめます：
\[ \hat{x}(n) = \sum_{k=0}^{p} a_k n^k = a_0 + a_1 n + a_2 n^2 + \cdots + a_p n^p \tag{1} \]最小二乗問題は、残差の二乗和を最小化することです：
\[ \min_{\mathbf{a}} \sum_{n=-m}^{m} \left(x_n - \sum_{k=0}^{p} a_k n^k \right)^2 \tag{2} \]行列表現と正規方程式
Vandermonde行列 \(\mathbf{A}\) を定義します：
\[ \mathbf{A} = \begin{bmatrix} (-m)^0 & (-m)^1 & \cdots & (-m)^p \\ (-m+1)^0 & (-m+1)^1 & \cdots & (-m+1)^p \\ \vdots & & \ddots & \vdots \\ m^0 & m^1 & \cdots & m^p \end{bmatrix} \tag{3} \]正規方程式の解は：
\[ \mathbf{a} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{x} \tag{4} \]平滑化された中心点の値は \(\hat{x}(0) = a_0\) であり、これは \(\mathbf{x}\) に対する線形変換として表せます：
\[ \hat{x}_0 = \mathbf{c}^T \mathbf{x} \tag{5} \]ここで \(\mathbf{c}\) は \((\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T\) の最初の行（\(a_0\) に対応）です。この \(\mathbf{c}\) がフィルタの畳み込み係数に相当します。
重要な性質
- 窓幅 \(2m+1\) と多項式次数 \(p\) の関係: \(p\) が大きいほど信号の形状を保存しますが、ノイズ除去能力は低下します。\(p = 0\) のとき、サビツキー・ゴーレイフィルタは単純移動平均と一致します。
- FIRフィルタとしての性質: 畳み込み係数 \(\mathbf{c}\) は窓幅と多項式次数のみで決まり、信号に依存しません。したがって、FIRフィルタの一種です。
- 微分の推定: \(a_1, a_2, \dots\) を取り出せば、信号の1次微分、2次微分も同時に推定できます。
周波数特性
サビツキー・ゴーレイフィルタの周波数応答は、畳み込み係数 \(\mathbf{c}\) のDTFT（離散時間フーリエ変換）で得られます：
\[ H(e^{j\omega}) = \sum_{n=-m}^{m} c_n e^{-j\omega n} \tag{6} \]EMA・移動平均との比較
| フィルタ | 通過帯域 | 遮断特性 | 信号形状保存 | 群遅延 |
|---|---|---|---|---|
| サビツキー・ゴーレイ | 広い | なだらか | 優秀 | 一定（線形位相） |
| EMA | パラメータ依存 | なだらか | 中程度 | 周波数依存 |
| 移動平均 | 窓幅依存 | なだらか | 低い | 一定（線形位相） |
| バターワース | 急峻（次数依存） | 急峻 | 設計次第 | 非線形 |
サビツキー・ゴーレイフィルタは線形位相特性を持つため、信号のピーク位置がシフトしません。これは EMAフィルタ（非線形位相）との大きな違いです。
Python実装
scipy.signal.savgol_filter
scipy.signalのsavgol_filterを使えば、1行でサビツキー・ゴーレイフィルタを適用できます。
import numpy as np
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt
# テストデータ生成
np.random.seed(42)
t = np.linspace(0, 1, 500)
# ピークを含む信号 + ノイズ
signal_clean = np.sin(2 * np.pi * 3 * t) + 0.5 * np.exp(-((t - 0.3) / 0.02)**2)
noise = 0.2 * np.random.randn(len(t))
signal_noisy = signal_clean + noise
# サビツキー・ゴーレイフィルタ適用
# window_length: 窓幅（奇数）, polyorder: 多項式次数
y_sg = savgol_filter(signal_noisy, window_length=21, polyorder=3)
# 移動平均との比較
kernel = np.ones(21) / 21
y_ma = np.convolve(signal_noisy, kernel, mode='same')
# プロット
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
axes[0].plot(t, signal_clean, 'g', label='Original', linewidth=2)
axes[0].plot(t, signal_noisy, 'gray', alpha=0.5, label='Noisy')
axes[0].plot(t, y_sg, 'b', label='Savitzky-Golay (w=21, p=3)')
axes[0].plot(t, y_ma, 'r--', label='Moving Average (w=21)')
axes[0].set_title('Smoothing Comparison')
axes[0].legend()
axes[0].grid(True)
# ピーク付近の拡大
mask = (t > 0.2) & (t < 0.4)
axes[1].plot(t[mask], signal_clean[mask], 'g', label='Original', linewidth=2)
axes[1].plot(t[mask], y_sg[mask], 'b', label='Savitzky-Golay')
axes[1].plot(t[mask], y_ma[mask], 'r--', label='Moving Average')
axes[1].set_title('Peak Region (Zoomed)')
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
ピーク付近を拡大すると、移動平均がピークを大幅に鈍らせるのに対し、サビツキー・ゴーレイフィルタはピークの位置と高さをほぼ保存していることが確認できます。
パラメータの影響
窓幅と多項式次数がフィルタの特性に与える影響を可視化します。
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 窓幅の影響（多項式次数固定: p=3）
for w in [5, 11, 21, 41]:
y = savgol_filter(signal_noisy, window_length=w, polyorder=3)
axes[0].plot(t, y, label=f'w={w}')
axes[0].plot(t, signal_clean, 'k--', alpha=0.5, label='Original')
axes[0].set_title('Effect of Window Length (polyorder=3)')
axes[0].legend()
axes[0].grid(True)
# 多項式次数の影響（窓幅固定: w=21）
for p in [0, 1, 2, 3, 5]:
y = savgol_filter(signal_noisy, window_length=21, polyorder=p)
axes[1].plot(t, y, label=f'p={p}')
axes[1].plot(t, signal_clean, 'k--', alpha=0.5, label='Original')
axes[1].set_title('Effect of Polynomial Order (window=21)')
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
考察:
- 窓幅が大きいほど強い平滑化が得られますが、信号の細かい構造が失われます
- 多項式次数 \(p = 0\) は単純移動平均と等価です
- \(p\) が窓幅に近いと、ノイズ除去効果がほとんどなくなります
- 実用的には \(p = 2\) または \(p = 3\) が多くの場面で良好です
微分の推定
サビツキー・ゴーレイフィルタの重要な特徴として、平滑化と同時に微分を推定できます。
# 1次微分と2次微分の推定
dy_sg = savgol_filter(signal_noisy, window_length=21, polyorder=3, deriv=1, delta=t[1]-t[0])
d2y_sg = savgol_filter(signal_noisy, window_length=21, polyorder=4, deriv=2, delta=t[1]-t[0])
# 真の微分（解析解）
dy_true = 2 * np.pi * 3 * np.cos(2 * np.pi * 3 * t) + \
0.5 * (-2 * (t - 0.3) / 0.02**2) * np.exp(-((t - 0.3) / 0.02)**2)
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
axes[0].plot(t, dy_true, 'g', label='True derivative')
axes[0].plot(t, dy_sg, 'b', label='SG 1st derivative')
axes[0].set_title('First Derivative Estimation')
axes[0].legend()
axes[0].grid(True)
axes[1].plot(t, d2y_sg, 'b', label='SG 2nd derivative')
axes[1].set_title('Second Derivative Estimation')
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
数値微分は通常ノイズを増幅しますが、サビツキー・ゴーレイフィルタは平滑化と微分を同時に行うため、安定した微分推定が可能です。
周波数応答の可視化
from scipy.signal import savgol_coeffs, freqz
# フィルタ係数の取得
for w, p in [(11, 2), (21, 3), (41, 3)]:
coeffs = savgol_coeffs(w, p)
freq, response = freqz(coeffs, worN=2048)
plt.plot(freq / np.pi, 20 * np.log10(np.abs(response)),
label=f'w={w}, p={p}')
plt.xlabel('Normalized Frequency (×π rad/sample)')
plt.ylabel('Gain (dB)')
plt.title('Savitzky-Golay Filter Frequency Response')
plt.legend()
plt.grid(True)
plt.ylim(-60, 5)
plt.tight_layout()
plt.show()
他のフィルタとの使い分け
| 用途 | 推奨フィルタ | 理由 |
|---|---|---|
| ピークを保存した平滑化 | サビツキー・ゴーレイ | 多項式フィッティングにより形状保存 |
| リアルタイム平滑化 | EMA | \(O(1)\)の計算コスト、因果的 |
| 特定周波数帯域の除去 | バターワース / ノッチ | 急峻な遮断特性 |
| ノイズ統計が既知 | ウィーナー / カルマン | 統計的最適性 |
| 微分の推定 | サビツキー・ゴーレイ | 平滑化微分が一体 |
| 非定常信号 | ウェーブレット | 時間-周波数解析 |
関連記事
- 指数移動平均（EMA）フィルタの周波数特性 - EMAは1次IIRフィルタとして計算コストが低いですが、サビツキー・ゴーレイのような信号形状保存能力はありません。
- 移動平均フィルタの種類と比較 - サビツキー・ゴーレイフィルタの多項式次数を0にすると単純移動平均と一致します。
- 高速フーリエ変換（FFT）の仕組みとPython実装 - フィルタの周波数特性を分析するための基礎となるアルゴリズムです。
- ローパスフィルタの設計と比較 - 各種ローパスフィルタの周波数応答・群遅延を比較しています。
- FIRフィルタとIIRフィルタの比較 - サビツキー・ゴーレイフィルタはFIRフィルタの一種であり、線形位相特性を持ちます。
- バターワースフィルタの設計とPython実装 - より急峻な遮断特性が必要な場合に使用するIIRフィルタです。
- ウィーナーフィルタの理論とPython実装 - 統計的最適性に基づくフィルタ設計。ノイズの統計情報が既知の場合に有効です。
- 窓関数とパワースペクトル密度（PSD）の理論とPython実装 - サビツキー・ゴーレイの窓幅選択と関連する窓関数の理論です。
- Matplotlib実践Tips：論文品質のグラフを作る - 周波数特性や平滑化結果のグラフを高品質に仕上げるためのTipsです。
- ノッチフィルタの設計とPython実装 - 特定周波数のみを除去するフィルタ。サビツキー・ゴーレイとは異なるアプローチです。
- カルマンフィルタの理論とPython実装 - 状態空間モデルに基づく確率的フィルタリング。サビツキー・ゴーレイとは根本的に異なるパラダイムです。
参考文献
- Savitzky, A., & Golay, M. J. E. (1964). Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 36(8), 1627-1639.
- Schafer, R. W. (2011). What Is a Savitzky-Golay Filter? IEEE Signal Processing Magazine, 28(4), 111-117.
- scipy.signal.savgol_filter documentation