分光器で測ったスペクトルや、センサーから流れ込んでくる電圧波形には、必ずと言っていいほどノイズが乗っています。このノイズを減らしたいとき、最初に思いつくのは「近くの点を平均する」移動平均でしょう。ところが移動平均をかけると、せっかくのスペクトルの鋭いピークが低くなり、裾が広がって潰れてしまいます。「ノイズは消したいが、ピークの高さや幅は守りたい」——この一見わがままな要求に、驚くほどうまく応えてくれるのが Savitzky-Golay（サビツキー・ゴーレイ）フィルタ です。
このフィルタの正体は、「窓の中のデータに低次の多項式を最小二乗で当てはめ、その多項式の値で中央点を置き換える」という操作です。さらに同じ枠組みで、平滑化しながら同時に微分（傾きや曲率）を推定できるという、ノイズに弱い数値微分の救世主にもなります。本記事は、この手法を「直感」「最小二乗の導出」「畳み込みカーネルへの帰着」「Python実装」の4本柱で、途中式を省略せずに解説します。
応用先は広く、たとえば次のような場面で標準的に使われています。
- 分光分析・クロマトグラフィー: ノイズの多いスペクトルを平滑化しつつ、ピーク位置・高さ・面積を保ったまま解析する（もともとSavitzkyとGolayが1964年に提案した動機がこれでした）
- センサ信号の前処理と微分: 加速度・ひずみ・電流などの時系列から、ノイズに埋もれた速度や変化率を安定して推定する
左図は「窓をずらしながら多項式をあてはめる」という操作そのもので、青い窓内の点に2次放物線（青線）を当てはめ、その窓中央の値（赤点）を出力として取り出している様子がわかる。右図では、同じノイズ信号に対して移動平均（オレンジ）はピークをなだらかに潰してしまうのに対し、Savitzky-Golay（青）は真の信号（黒破線）のピーク高さをほぼ保ったまま滑らかにしている。これがフィルタの本質を一枚で示している。
本記事の内容
- 移動平均との対比でみるSavitzky-Golayフィルタの直感
- 窓内データへの多項式最小二乗フィッティングの定式化
- 正規方程式 $\bm{A}^\top \bm{A}\, \bm{c} = \bm{A}^\top \bm{y}$ の導出
- 係数行列 $(\bm{A}^\top \bm{A})^{-1}\bm{A}^\top$ の中央行が畳み込みカーネルになる仕組み
- 平滑化と微分推定の統一的な係数導出
- numpyによる係数の自前計算と
scipy.signal.savgol_filterとの一致確認 - 移動平均とのピーク保存性の比較、ノイズ信号の傾き推定の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
Savitzky-Golayフィルタは「係数が固定されたFIRフィルタの一種」であり、その係数を最小二乗法で決めるという二段構えの理論です。FIRフィルタの基礎（畳み込みと係数の意味）と、最小二乗法・線形代数の基本的な行列演算を押さえておくと、以降の導出がすっと頭に入ります。
Savitzky-Golayフィルタとは
ノイズだらけの測定値の列を、窓を1点ずつ滑らせながら眺めている状況を想像してください。素朴な移動平均は、窓の中の値を「全部足して個数で割る」だけです。これは暗黙のうちに「窓の中ではデータは平らな水平線（0次多項式 = 定数）だ」と仮定して、その水平線の高さを最小二乗で求めていることに相当します。実際、定数 $a_0$ で近似したときの二乗誤差 $\sum (y_i – a_0)^2$ を最小にする $a_0$ は、まさに窓内の平均値です。
ところが現実の信号は、ピークの近くでは平らどころか山なりに曲がっています。平らな線で無理やり近似すれば、山のてっぺんは過小評価され、麓は過大評価され、結果としてピークがなまされてしまいます。
Savitzky-Golayフィルタのアイデアはシンプルです。「窓の中を水平線で近似するのをやめ、もっと表現力のある 低次多項式（2次や4次の放物線的な曲線）で近似しよう」というものです。山なりの形を放物線で受け止められるので、ピークの高さや曲率が保たれます。そして、当てはめた多項式の窓中央での値を、その点の平滑化後の出力とします。窓を1つずらすたびに多項式を当て直し、また中央の値を取る——これを繰り返すのがSavitzky-Golayフィルタです。
ここで素晴らしいのは、後で導出するように、窓の位置をずらしても「多項式フィッティング → 中央値を取る」という操作は毎回まったく同じ重み係数の重み付き和に帰着する点です。つまり、毎回その場で最小二乗を解く必要はなく、一度きり計算した係数列との 畳み込み（FIRフィルタ）として実装できます。理論は最小二乗、実装は固定係数の畳み込み、という美しい二面性がこのフィルタの魅力です。
次節では、この「窓内多項式フィッティング」を行列で厳密に定式化し、最小二乗解を導きます。
窓内データへの多項式フィッティングの定式化
窓の幅（窓長）を奇数 $2m+1$ とし、注目している中央点を原点にとります。窓内のサンプルのインデックスを $i = -m, -m+1, \dots, -1, 0, 1, \dots, m$ とし、対応する観測値を $y_i$ とします。インデックスを原点中心に取るのがポイントで、これにより「中央での値」が極めて簡単な形になります。
この $2m+1$ 個の点に、次数 $d$ の多項式
$$ p(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_d x^d = \sum_{j=0}^{d} a_j x^j $$
を当てはめます。ここで $x$ は窓内のローカル座標（整数インデックス $i$）です。通常はサンプル間隔を $1$ とみなして $x = i$ とします。次数は窓長より小さく、$d < 2m+1$ を満たす必要があります（そうでないと係数が一意に決まりません）。
各サンプル点 $x = i$ での多項式の値と観測値 $y_i$ の食い違い（残差）の二乗和
$$ E = \sum_{i=-m}^{m} \left( p(i) – y_i \right)^2 = \sum_{i=-m}^{m} \left( \sum_{j=0}^{d} a_j\, i^j – y_i \right)^2 $$
を最小化する係数 $a_0, a_1, \dots, a_d$ を求めるのが、最小二乗フィッティングです。
これを行列で書くと見通しがよくなります。サンプル点でのべき乗を並べた 設計行列（Vandermonde行列） $\bm{A}$ を
$$ \bm{A} = \begin{pmatrix} 1 & (-m) & (-m)^2 & \cdots & (-m)^d \\ 1 & (-m+1) & (-m+1)^2 & \cdots & (-m+1)^d \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & m & m^2 & \cdots & m^d \end{pmatrix} $$
と定義します。$\bm{A}$ は $(2m+1) \times (d+1)$ 行列で、第 $i$ 行・第 $j$ 列の成分は $A_{ij} = i^{\,j}$（行は $i=-m,\dots,m$、列は $j=0,\dots,d$）です。係数ベクトルを $\bm{c} = (a_0, a_1, \dots, a_d)^\top$、観測値ベクトルを $\bm{y} = (y_{-m}, \dots, y_m)^\top$ とすると、多項式の各点での値はちょうど $\bm{A}\bm{c}$ で、残差二乗和は
$$ E = \| \bm{A}\bm{c} – \bm{y} \|^2 = (\bm{A}\bm{c} – \bm{y})^\top (\bm{A}\bm{c} – \bm{y}) $$
とコンパクトにまとまります。ここまでで問題は「$\| \bm{A}\bm{c} – \bm{y} \|^2$ を最小化する $\bm{c}$ を求めよ」という標準的な線形最小二乗問題に帰着しました。次節でこれを解きます。
左のヒートマップは窓長5・2次の設計行列 $\bm{A}$（5×3）で、各行が $(1,\ i,\ i^2)$ というVandermonde構造になっていることが色で一目でわかる。右のヒートマップは擬似逆行列 $\bm{P}=(\bm{A}^\top\bm{A})^{-1}\bm{A}^\top$（3×5）で、第0行（$a_0$行）が平滑化カーネル、第1行（$a_1$行）が1階微分カーネル、第2行（$a_2$行）が2階微分カーネルに対応する——この3行の形の違いが、平滑化と微分の統一的な枠組みを示している。
正規方程式の導出
ゴールは、$E(\bm{c}) = \| \bm{A}\bm{c} – \bm{y} \|^2$ を最小にする $\bm{c}$ を閉じた形で求めることです。多変数関数の最小値では、すべての偏微分がゼロになる（勾配がゼロベクトルになる）ことを使います。
まず $E$ を展開します。スカラーの転置はそれ自身に等しいこと（$\bm{c}^\top \bm{A}^\top \bm{y}$ はスカラーなので $\bm{y}^\top \bm{A}\bm{c}$ に等しい）を使うと、
$$ \begin{align} E(\bm{c}) &= (\bm{A}\bm{c} – \bm{y})^\top (\bm{A}\bm{c} – \bm{y}) \\ &= \bm{c}^\top \bm{A}^\top \bm{A}\bm{c} – \bm{c}^\top \bm{A}^\top \bm{y} – \bm{y}^\top \bm{A}\bm{c} + \bm{y}^\top \bm{y} \\ &= \bm{c}^\top \bm{A}^\top \bm{A}\bm{c} – 2\,\bm{c}^\top \bm{A}^\top \bm{y} + \bm{y}^\top \bm{y} \end{align} $$
となります。2行目から3行目では、$\bm{c}^\top \bm{A}^\top \bm{y}$ と $\bm{y}^\top \bm{A}\bm{c}$ が同じスカラーなので、まとめて $-2\,\bm{c}^\top \bm{A}^\top \bm{y}$ にしました。
次に $\bm{c}$ について勾配（ベクトル微分）を取ります。ベクトル微分の公式 $\dfrac{\partial}{\partial \bm{c}}(\bm{c}^\top \bm{M} \bm{c}) = 2\bm{M}\bm{c}$（$\bm{M}$ が対称のとき）と $\dfrac{\partial}{\partial \bm{c}}(\bm{c}^\top \bm{b}) = \bm{b}$ を使います。ここで $\bm{M} = \bm{A}^\top \bm{A}$ は対称行列なので、
$$ \frac{\partial E}{\partial \bm{c}} = 2\,\bm{A}^\top \bm{A}\,\bm{c} – 2\,\bm{A}^\top \bm{y} $$
が得られます。これをゼロベクトルとおくと、
$$ \bm{A}^\top \bm{A}\,\bm{c} = \bm{A}^\top \bm{y} $$
という 正規方程式（normal equations） が導かれます。$\bm{A}$ の列が線形独立（$d < 2m+1$ かつサンプル点が相異なるなら、Vandermonde行列の性質から成立）であれば、$\bm{A}^\top \bm{A}$ は $(d+1)\times(d+1)$ の正則な対称正定値行列なので逆行列を持ち、最小二乗解は
$$ \boxed{\ \bm{c} = (\bm{A}^\top \bm{A})^{-1} \bm{A}^\top\, \bm{y}\ } $$
と一意に定まります。$E$ は $\bm{c}$ について二次の凸関数（$\bm{A}^\top\bm{A}$ が正定値）なので、この停留点はまぎれもなく最小値です。
ここで現れた $(d+1) \times (2m+1)$ 行列
$$ \bm{P} = (\bm{A}^\top \bm{A})^{-1} \bm{A}^\top $$
を 擬似逆行列（Moore-Penroseの擬似逆） と呼びます。重要なのは、$\bm{P}$ が観測値 $\bm{y}$ にまったく依存せず、窓長 $2m+1$ と次数 $d$ だけで決まる点です。窓をどこに動かそうと、$\bm{A}$ も $\bm{P}$ も同じです。この「$\bm{y}$ に依らない」性質が、次節で畳み込みカーネルを生み出します。
中央行が畳み込みカーネルになる仕組み
求めた係数ベクトル $\bm{c} = \bm{P}\bm{y}$ の各成分は、当てはめた多項式 $p(x) = \sum_j a_j x^j$ の係数 $a_0, a_1, \dots, a_d$ そのものです。私たちが欲しいのは、この多項式の窓中央 $x=0$ での値でした。原点中心にインデックスを取っておいたおかげで、中央での値は
$$ p(0) = a_0 + a_1 \cdot 0 + a_2 \cdot 0^2 + \cdots = a_0 $$
と、ただの $a_0$ になります。これがインデックスを原点中心に取った最大のご利益です。$a_0$ は係数ベクトル $\bm{c}$ の第0成分なので、行列 $\bm{P}$ の 第0行（中央行） を $\bm{P}$ から取り出せば、
$$ \hat{y}_0 = a_0 = \sum_{i=-m}^{m} P_{0,i}\, y_i $$
つまり、平滑化後の中央値は、観測値 $y_i$ の 重み付き和 で書けます。この重み列
$$ g_i = P_{0,i}, \quad i = -m, \dots, m $$
こそが Savitzky-Golay係数 であり、フィルタの畳み込みカーネルです。
ここで注目すべきは、窓を信号上で1サンプルずらしても、$\bm{P}$ は不変なので 同じ重み列 $g_i$ を使い回せることです。信号 $\{s[n]\}$ の各点 $n$ に対し、その周りの窓 $s[n-m], \dots, s[n+m]$ に重み $g_i$ をかけて足すと、平滑化出力
$$ \hat{s}[n] = \sum_{i=-m}^{m} g_i\, s[n+i] $$
が得られます。これはまさにFIRフィルタの畳み込み（相互相関）の形です。「窓ごとに最小二乗を解き直す」という重そうな処理が、たった1本の固定カーネルとの畳み込みに化けたわけです。実装の軽さと理論の正しさが両立する、Savitzky-Golayフィルタの核心がここにあります。
微分推定への一般化
同じ枠組みで、平滑化だけでなく導関数も推定できます。当てはめた多項式 $p(x)=\sum_j a_j x^j$ を $k$ 回微分して $x=0$ を代入すると、テイラー展開の係数の関係から
$$ p^{(k)}(0) = k!\, a_k $$
となります。たとえば1階微分（傾き）は $p'(0) = a_1$、2階微分（曲率）は $p”(0) = 2a_2$ です。$a_k$ は $\bm{P}$ の第 $k$ 行に対応するので、$\bm{P}$ の 第 $k$ 行 を取り出して $k!$ を掛ければ、$k$ 階導関数を推定する畳み込みカーネルが得られます。
$$ \widehat{s^{(k)}}[n] = \frac{k!}{h^k} \sum_{i=-m}^{m} P_{k,i}\, s[n+i] $$
ここで $h$ はサンプル間隔です（$x=i$ を物理量に直すと $x = i\,h$ なので、$k$ 階微分には $h^{-k}$ のスケーリングが必要になります）。素朴な差分による数値微分はノイズを激しく増幅しますが、Savitzky-Golay微分は「ノイズに頑健な多項式当てはめ」を経由するため、平滑化と微分を同時にこなせます。これがセンサ信号の速度・加速度推定で重宝される理由です。
次節では、小さな具体例で実際に係数を手計算的に確かめ、その後Pythonで一般のケースを計算します。
具体例: 窓長5・2次多項式の係数
最も使われる設定の一つ、窓長 $2m+1 = 5$（$m=2$）、次数 $d=2$ のSavitzky-Golay平滑化係数を求めてみましょう。インデックスは $i = -2, -1, 0, 1, 2$ です。設計行列は
$$ \bm{A} = \begin{pmatrix} 1 & -2 & 4 \\ 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & 2 & 4 \end{pmatrix} $$
です（第 $i$ 行が $(1,\ i,\ i^2)$）。$\bm{A}^\top \bm{A}$ を計算します。各成分は $\sum_i i^p$ の形になります。対称性から奇数べきの和はゼロ（$\sum i = 0,\ \sum i^3 = 0$）になることに注意すると、
$$ \bm{A}^\top \bm{A} = \begin{pmatrix} \sum i^0 & \sum i^1 & \sum i^2 \\ \sum i^1 & \sum i^2 & \sum i^3 \\ \sum i^2 & \sum i^3 & \sum i^4 \end{pmatrix} = \begin{pmatrix} 5 & 0 & 10 \\ 0 & 10 & 0 \\ 10 & 0 & 34 \end{pmatrix} $$
ここで $\sum i^0 = 5$、$\sum i^2 = 4+1+0+1+4 = 10$、$\sum i^4 = 16+1+0+1+16 = 34$ を使いました。奇数べきの和が消えるおかげで、行列がブロック的に分かれて逆行列が求めやすくなります。
私たちが必要なのは中央値、すなわち $a_0$ に対応する $\bm{P}$ の第0行です。$\bm{c} = (\bm{A}^\top\bm{A})^{-1}\bm{A}^\top \bm{y}$ の第0成分を取り出すため、$(\bm{A}^\top\bm{A})^{-1}$ の第0行が必要です。$3\times3$ 行列の逆行列の第0行は、余因子から計算できます。行列式は
$$ \det(\bm{A}^\top\bm{A}) = 5(10\cdot34 – 0) – 0 + 10(0 – 10\cdot10) = 5\cdot340 – 10\cdot100 = 1700 – 1000 = 700 $$
第0行の各成分（逆行列の $(0,0),(0,1),(0,2)$）は、余因子を $\det$ で割って
$$ (\bm{A}^\top\bm{A})^{-1}_{0,\cdot} = \frac{1}{700}\,(\,10\cdot34-0,\ \ -(0\cdot34-0),\ \ 0-10\cdot10\,) = \frac{1}{700}(340,\ 0,\ -100) $$
となります。これを $\bm{A}^\top$ の各列にかけると、各 $i$ に対する重み $g_i$ が出ます。$\bm{A}^\top$ の第 $i$ 列は $(1,\ i,\ i^2)^\top$ なので、
$$ g_i = \frac{1}{700}\left( 340\cdot 1 + 0\cdot i + (-100)\cdot i^2 \right) = \frac{340 – 100\,i^2}{700} = \frac{17 – 5\,i^2}{35} $$
と、$i$ の関数として一発で書けます。$i = -2,-1,0,1,2$ を代入すると、
$$ g_{-2} = \frac{17-20}{35} = -\frac{3}{35},\quad g_{-1} = \frac{17-5}{35} = \frac{12}{35},\quad g_{0} = \frac{17}{35},\quad g_{1} = \frac{12}{35},\quad g_{2} = -\frac{3}{35} $$
すなわち係数列は
$$ \bm{g} = \frac{1}{35}(-3,\ 12,\ 17,\ 12,\ -3) $$
です。これは文献で「5点2次のSavitzky-Golay平滑化係数」として知られる有名な値とぴったり一致します。重みの和は $\frac{-3+12+17+12-3}{35} = \frac{35}{35} = 1$ で、定数信号をそのまま通す（DCゲイン1）ことも確認できます。
注目すべきは、両端の重み $g_{\pm 2} = -3/35$ が負であることです。単純な移動平均なら全ての重みが等しく正（$1/5 = 0.2$）ですが、Savitzky-Golayでは中央を重く（$17/35 \approx 0.486$）、両端をわずかに負にすることで、放物線的な山なり形に追従できます。この「中央を持ち上げ端を引く」重み形状こそ、ピークを潰さない秘密です。次節ではこれをnumpyで自動計算し、scipy と照合します。
左グラフは、本記事で導出した擬似逆行列の中央行から得た5点2次の係数（青）と scipy.signal.savgol_coeffs の出力（オレンジ）を並べたもので、両者は数値的に完全に一致している。右グラフは窓長11で次数を1〜5と変えたときの係数総和で、どの次数でも総和が1.0000になることを確認できる——これは「定数信号（DC成分）がそのまま通る」という直感に一致しており、フィルタが信号の基線を変えないことを保証している。
Pythonでの実装
係数の自前計算とscipyとの一致確認
まず、これまで導いた $\bm{P} = (\bm{A}^\top\bm{A})^{-1}\bm{A}^\top$ の行を取り出す手続きを、そのままnumpyで実装します。先ほど手計算した5点2次の係数を再現できるか確かめます。
import numpy as np
from math import factorial
from scipy.signal import savgol_coeffs
def sg_coeffs(window_length, polyorder, deriv=0, delta=1.0):
"""Savitzky-Golay係数を最小二乗の定義から自前計算する。
window_length: 窓長（奇数 2m+1）
polyorder: 多項式次数 d
deriv: 何階微分の係数か（0なら平滑化）
delta: サンプル間隔 h
"""
assert window_length % 2 == 1, "窓長は奇数"
m = (window_length - 1) // 2
i = np.arange(-m, m + 1) # 窓内インデックス -m..m
# 設計行列 A: 各行 (i^0, i^1, ..., i^d)
A = np.vander(i, polyorder + 1, increasing=True)
# 擬似逆 P = (A^T A)^{-1} A^T
P = np.linalg.inv(A.T @ A) @ A.T
# deriv階微分の中央値カーネル = P の第deriv行 × deriv! / h^deriv
g = P[deriv] * (factorial(deriv) / (delta ** deriv))
return g
# 5点・2次の平滑化係数
g = sg_coeffs(5, 2, deriv=0)
print("自前計算 :", g * 35, "/ 35") # 35倍して整数を見る
print("scipy :", savgol_coeffs(5, 2) * 35, "/ 35")
print("一致 :", np.allclose(g, savgol_coeffs(5, 2)))
出力では 自前計算 : [-3. 12. 17. 12. -3.] / 35 のように、手計算で得た $\frac{1}{35}(-3,12,17,12,-3)$ が再現され、scipy の savgol_coeffs(5, 2) とも True で完全一致します。これにより、「擬似逆行列の中央行がSavitzky-Golay係数である」という導出が正しいことを、コード上でも確認できました。np.vander(..., increasing=True) が設計行列そのものになっている点も、定義との対応がそのまま読み取れます。
3つのカーネルを並べると、それぞれの設計哲学の違いが視覚的に浮かび上がる。次数2（左）は中央が正の山で両端が負という基本的な「持ち上げ＋引き下げ」の形。次数4（中央）はより複雑な波形で、中央の鋭いピークをより精密に再現できる分、両端に振動が生じる。移動平均（右）はすべての重みが均一（$1/21$）で、いかなる形にも「鈍感」であることがわかる。重みの形状の違いが、そのままピーク保存性の差として現れる。
移動平均とのピーク保存性の比較
次に、このフィルタの売りである「ピークを潰さない」性質を、ガウス型のピークを持つノイズ信号で検証します。比較対象は同じ窓長の単純移動平均です。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
# ガウス型ピークを持つクリーン信号 + ノイズ
np.random.seed(0)
x = np.linspace(-10, 10, 400)
clean = np.exp(-(x**2) / 2.0) # 鋭いガウスピーク（高さ1）
noisy = clean + 0.05 * np.random.randn(x.size)
win = 31 # 共通の窓長
# Savitzky-Golay（4次）
sg = savgol_filter(noisy, win, polyorder=4)
# 単純移動平均（同じ窓長）
kernel = np.ones(win) / win
ma = np.convolve(noisy, kernel, mode="same")
plt.figure(figsize=(10, 5))
plt.plot(x, noisy, color="0.8", lw=1, label="noisy")
plt.plot(x, clean, "k--", lw=1.2, label="true peak")
plt.plot(x, ma, "tab:orange", lw=1.8, label=f"moving average (win={win})")
plt.plot(x, sg, "tab:blue", lw=1.8, label=f"Savitzky-Golay (win={win}, d=4)")
plt.xlabel("x"); plt.ylabel("amplitude")
plt.title("Peak preservation: SG vs moving average")
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print(f"true peak = {clean.max():.3f}")
print(f"MA peak = {ma.max():.3f}")
print(f"SG peak = {sg.max():.3f}")
グラフと出力から、2つのことがはっきり読み取れます。第一に、移動平均（オレンジ）はピークの高さを大きく削り、MA peak が真値1より明らかに低くなり、山が横に広がって裾が膨らみます。第二に、Savitzky-Golay（青）はノイズを同程度に抑えながらも、SG peak が真値にほぼ張り付き、山なりの形をきれいに保ちます。これは、移動平均が窓内を「平ら」と仮定するのに対し、Savitzky-Golayが4次多項式で山の曲率まで表現できるからです。前節で見た「中央を重く、端を負に」という重み形状が、まさにこのピーク保存を生んでいます。
左図では波形レベルで、Savitzky-Golay（青）が真の信号（黒破線）をほぼ忠実になぞる一方、移動平均（オレンジ）はピークを大きく削って裾を広げているのが見てとれる。右の棒グラフはピーク高さを数値で比較したもので、移動平均は真値1.000から大きく下がるのに対し、SGは1.000に極めて近い値を保持している——この差は窓内の近似モデルの「表現力」がそのまま精度に現れた結果だ。
ノイズ信号の傾き（微分）推定
最後に、Savitzky-Golayの真骨頂である「平滑化しながらの微分」を試します。素朴な差分（隣接点の引き算）と比べて、どれだけノイズに強いかを見ます。対象は傾きが既知の信号 $s(x) = \sin x$ で、その導関数 $\cos x$ をノイズ越しに復元できるか確かめます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
np.random.seed(1)
x = np.linspace(0, 4 * np.pi, 500)
h = x[1] - x[0] # サンプル間隔
clean = np.sin(x)
noisy = clean + 0.05 * np.random.randn(x.size)
true_deriv = np.cos(x) # 真の1階微分
# 素朴な数値微分（中心差分）
naive = np.gradient(noisy, h)
# Savitzky-Golay微分（deriv=1, delta=h でスケーリング込み）
sg_deriv = savgol_filter(noisy, window_length=21, polyorder=3,
deriv=1, delta=h)
plt.figure(figsize=(10, 5))
plt.plot(x, naive, color="0.7", lw=1, label="naive finite difference")
plt.plot(x, true_deriv, "k--", lw=1.5, label="true derivative cos(x)")
plt.plot(x, sg_deriv, "tab:red", lw=1.8, label="Savitzky-Golay derivative")
plt.xlabel("x"); plt.ylabel("ds/dx")
plt.title("Noisy derivative estimation")
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# RMS誤差で比較
rms_naive = np.sqrt(np.mean((naive - true_deriv) ** 2))
rms_sg = np.sqrt(np.mean((sg_deriv - true_deriv) ** 2))
print(f"RMS error (naive) = {rms_naive:.3f}")
print(f"RMS error (SG) = {rms_sg:.3f}")
このグラフから、素朴な差分（灰色）は元信号のわずかなノイズを激しく増幅し、真の $\cos x$（黒破線）がまったく見えないほどギザギザに暴れることがわかります。一方Savitzky-Golay微分（赤）は、滑らかな曲線として真の導関数をほぼ正確になぞります。出力の RMS error でも、SGの誤差は素朴差分より一桁以上小さくなります。これは、微分が「擬似逆行列 $\bm{P}$ の第1行 $\times\,1!/h$」という固定カーネルとの畳み込みであり、多項式当てはめを経由することでノイズが平均化されるためです。差分のように隣接2点だけを使うのではなく、窓内の全点を最小二乗的に使うことが、ノイズ耐性の源泉です。
上段は元の $\sin x$ 信号（ノイズあり）で、下段がその微分推定の比較だ。素朴な差分（灰色）はノイズが完全にランダムに暴れてしまい、真の $\cos x$（黒破線）を追う気配すらない。SG微分（赤）は窓内21点を多項式で平均化するため、ノイズが打ち消し合い、なめらかな $\cos x$ の曲線として復元されている。右上のRMS誤差の数値がその差を定量的に示している。
係数の形状を可視化する
最後に、平滑化カーネル・1階微分カーネルの「形」そのものを描いて、重みがどう信号を組み合わせているかを直感的に掴みます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_coeffs
win = 21
i = np.arange(-(win // 2), win // 2 + 1)
smooth = savgol_coeffs(win, polyorder=4, deriv=0) # 平滑化カーネル
deriv1 = savgol_coeffs(win, polyorder=4, deriv=1) # 1階微分カーネル
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].stem(i, smooth)
ax[0].set_title("SG smoothing kernel (win=21, d=4)")
ax[0].set_xlabel("index i"); ax[0].set_ylabel("weight g_i")
ax[0].grid(alpha=0.3)
ax[1].stem(i, deriv1)
ax[1].set_title("SG 1st-derivative kernel")
ax[1].set_xlabel("index i"); ax[1].set_ylabel("weight")
ax[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("平滑化カーネルの総和 :", smooth.sum()) # ≈ 1（DCを通す）
print("微分カーネルの総和 :", deriv1.sum()) # ≈ 0（定数の微分は0）
2つのカーネルの形から、設計思想がそのまま読み取れます。平滑化カーネル（左）は中央が正の山で対称な「丘」の形をしており、重みの総和が1になるため定数（DC成分）をそのまま通します。微分カーネル（右）は原点について反対称で、左半分が負・右半分が正——これはまさに「右の値から左の値を引く」差分の一般化であり、総和が0になるため定数信号の微分を正しく0と返します。前者がType I、後者がType IIIの線形位相FIRに対応している点は、ディジタルフィルタの基礎（FIR/IIR）を解説で扱った対称性の議論とそのままつながります。
左の平滑化カーネル（青）は完全な左右対称で「中央が最大・端が小さい」丘の形をしており、テキストボックスに示す通り総和が1.0000になっている。右の1階微分カーネル（赤）は逆に原点に関して完全に反対称（符号が反転する点対称）で、総和はほぼ0になる。この対称/反対称の違いは数学的に「偶関数を合成すると偶関数、奇関数を合成すると奇関数」という畳み込みの性質から来ており、平滑化と微分という2つの異なる操作が同じ擬似逆行列の異なる行から生まれる理由がここに凝縮されている。
周波数特性と窓長・次数のトレードオフ
Savitzky-Golayフィルタは「ローパスフィルタ」の一種です。しかし、同じ窓長であっても次数の選び方によって、どの周波数成分をどれだけ通すかが大きく変わります。このトレードオフを理解することが、実際の信号に適切なパラメータを選ぶ上で欠かせません。
周波数応答の比較
カーネルの形が決まれば、その周波数応答（フィルタがどの周波数成分をどれだけ通すか）は、カーネルのフーリエ変換で与えられます。
$$ H(f) = \sum_{i=-m}^{m} g_i\, e^{-j 2\pi f i} $$
移動平均は $\text{sinc}$ 関数状の周波数応答を持ち、ストップバンドでリップル（漏れ）が大きくなります。一方Savitzky-Golayは、次数が高いほど通過帯域での応答が平坦（ほぼ1）になり、急峻な変化もよく保存します。ただし、次数を上げるほどカットオフが緩くなり、ノイズ除去能力は下がります。
左図（リニアスケール）では、移動平均（オレンジ）の周波数応答がストップバンドで大きく振動しているのに対し、SG次数4（青）は通過帯域（低周波）でほぼ平坦1.0を保ちながらなだらかに減衰していることがわかる。右図（dBスケール）では移動平均がストップバンドで-50dB程度のリップルを示す一方、SGはより単調に減衰している。「鋭いピークを保つ」性質が周波数領域では「通過帯域の平坦性」として現れており、スペクトル解析における優位性の根拠がここにある。
窓長と次数のトレードオフ
Savitzky-Golayフィルタのパラメータは主に2つ——窓長 $2m+1$ と多項式の次数 $d$ です。この2つの組み合わせが、「滑らかさ（ノイズ除去）」と「忠実度（ピーク保存）」のバランスを決めます。
- 窓長を大きくする: より多くの点を使うため、ノイズ除去効果が上がる。ただし、急峻な変化（鋭いピーク）を鈍らせる傾向がある。
- 次数を高くする: 高次の多項式はより複雑な形状を表現できるため、ピーク保存性が上がる。ただし、窓内の自由度が増す分、ノイズ除去効果が下がる。
ガウスピークの復元実験で、この違いを視覚的に確認できます。
上段は次数3を固定し窓長を7・21・51と変えた結果だ。窓長7では信号がガタガタで平滑化が不十分、窓長51では過度に平滑化されてピークがやや低下する。窓長21がこの信号では最もバランスが良い。下段は窓長21を固定して次数を1・3・5と変えた結果で、次数1ではピークがかなり潰れるが、次数3・5ではピークをほぼ忠実に保ちながら滑らかさも確保している。実際のチューニングは「窓長を信号の変化スケールに合わせ、次数は2〜5の間で必要に応じて上げる」という方針が基本になる。
応用例：スペクトル平滑化と複数ピークの保存
最初に挙げた分光分析のユースケースを、複数の鋭いピークを持つ合成スペクトルで再現してみます。
上段のスペクトルでは、高さの異なる4本のピークが近接して並んでいる。移動平均（オレンジ）は低い・狭いピーク（$x\approx17$付近）を大きく潰してしまい、隣のピークと区別できなくなりかけている。一方SG（青）はすべてのピーク高さをほぼ保ったまま基線のノイズだけを除去している。下段の棒グラフで各ピークの高さを数値比較すると、移動平均では特に幅の狭い鋭いピークで真値からの乖離が大きいことが明確にわかる。狭いピークほど移動平均の「局所平均化」の影響を受けやすく、Savitzky-Golayの優位性が際立つ領域だ。
ここまでで、理論から実装、さらに周波数特性とトレードオフまで一通り見渡せました。次節のまとめでキーポイントを整理します。
まとめ
本記事では、Savitzky-Golayフィルタを最小二乗フィッティングの観点から導出し、Pythonで実装しました。
- 直感: 窓内を移動平均のように「平ら」と仮定するのではなく、低次多項式で近似することで、ピークの高さや曲率を保ったまま平滑化できる
- 定式化: 窓内 $2m+1$ 点への次数 $d$ の多項式当てはめは、Vandermonde設計行列 $\bm{A}$ を使った線形最小二乗問題 $\min \| \bm{A}\bm{c}-\bm{y}\|^2$ になる
- 正規方程式: 勾配をゼロとおくと $\bm{A}^\top\bm{A}\,\bm{c} = \bm{A}^\top\bm{y}$ が導かれ、解は $\bm{c} = (\bm{A}^\top\bm{A})^{-1}\bm{A}^\top\bm{y}$
- 畳み込みカーネル化: 原点中心インデックスにより中央値が $a_0$ となり、擬似逆 $\bm{P}=(\bm{A}^\top\bm{A})^{-1}\bm{A}^\top$ の 中央行がそのままSavitzky-Golay係数（FIRカーネル）になる。$\bm{P}$ は観測値に依存しないため、窓ごとに解き直さず1本の畳み込みで実装できる
- 微分推定: $\bm{P}$ の 第 $k$ 行に $k!/h^k$ を掛けたカーネルで、$k$ 階導関数をノイズに頑健に推定できる
- 検証: 5点2次で手計算した $\frac{1}{35}(-3,12,17,12,-3)$ が
scipy.signal.savgol_coeffsと一致し、移動平均よりピークを保ち、素朴な差分より微分のRMS誤差が一桁小さいことを確認した
Savitzky-Golayフィルタは「理論は最小二乗、実装は固定係数のFIR」という二面性を持つ、前処理の定番ツールです。窓長と次数の選び方（次数が高いほどピーク追従性は上がるがノイズ除去は弱まる）というトレードオフを意識すると、実データでの調整が楽になります。
次のステップとして、以下の記事も参考にしてください。