- Savitzky-Golayの係数を計算するsavgol_coeffsの内容に関するメモ．
Referencesの"Ying, and Jing Bai. 2005. Savitzky-Golay smoothing and differentiation filter for even number data. Signal Process. 85, 7 (July 2005), 1429-1434."と，その中で引かれている"S.J. Orfanidis, Introduction to Signal Processing, Prentice- Hall, EnglewoodCliffs, NJ, 1996 (Chapter8)." を参考にした．
【関連記事】
- [A] 線形回帰 - Notes_JP
Savitzky-Golayでやること
窓範囲に含まれる1点を，窓範囲に含まれるデータを使ってフィルタリングする．この記事では- 窓の長さ：$M=\text{window\_length}$
- フィルタリング結果として返す窓の座標：$0 \leq p = \text{pos} < M$
また，着目している窓範囲のデータは
\{(x_{i}, y_{i})\}_{i = 0}^{M - 1}
\end{aligned}
savgol_coeffsでやっていること
後で定義する行列$X$に対して，$X^{T} \bm{c} = k! \bm{e}_{k}$の解$c$を係数として返している．ただし，\bm{e}_{k}
=
\left(
\begin{array}{c}
0 \\
\vdots \\
1\\
\vdots \\
0
\end{array}
\right)
(k ~~~
\end{aligned}
というわけで，やっていることは超単純．
以下，なぜこれで係数が計算できるのかを確認する．
導出
データ$\{(x_{i}, y_{i})\}_{i = 0}^{M - 1}$を$N$次多項式y = \sum_{j = 0}^{N} \beta_{j} x^{j}
\end{aligned}
こうすると，フィルタリング結果
\hat{y}_{p}
\end{aligned}
\frac{\mathrm{d}^{k}\hat{y}_{p}}{\mathrm{d}x^{k}}
\end{aligned}
スムージング
後で見る$k$回微分の結果の$k = 0$の場合として以下の結果も含まれるが，まず微分しない場合を理解しておくとラク．データはまとめて
\underbrace{
\begin{pmatrix}
y_{0} \\
\vdots \\
y_{M - 1} \\
\end{pmatrix}}_{= \bm{y}}
=
\underbrace{
\begin{pmatrix}
(-p)^{0} & (-p)^{1} & \cdots & (-p)^{N} \\
(1-p)^{0} & (1-p)^{1} & \cdots & (1-p)^{N} \\
\vdots & \vdots & \ddots & \vdots \\
(M - 1-p)^{0} & (M - 1-p)^{1} & \cdots & (M - 1-p)^{N} \\
\end{pmatrix}}_{= X(-p)}
\underbrace{
\begin{pmatrix}
\beta_{0} \\
\vdots \\
\beta_{N} \\
\end{pmatrix}}_{= \bm{\beta}}
\end{aligned}
パラメータ$\boldsymbol{\beta}$の推定結果は，
\hat{\boldsymbol{\beta}}
&= \underbrace{(X^{T} X)^{-1} X^{T}}_{=G^{T}} \boldsymbol{y}
\end{aligned}
ここで，$G = (\bm{g}_{0}, \ldots, \bm{g}_{M - 1})$とすると
\hat{\boldsymbol{\beta}}
&=
\begin{pmatrix}
\bm{g}_{0}^{T} \\
\vdots \\
\bm{g}_{M - 1}^{T}
\end{pmatrix} \boldsymbol{y}
=
\begin{pmatrix}
\bm{g}_{0}^{T} \boldsymbol{y} \\
\vdots \\
\bm{g}_{M - 1}^{T} \boldsymbol{y}
\end{pmatrix}
\end{aligned}
データ$\boldsymbol{y}$のフィルタリング結果は$x = 0$における値なので
\hat{y}_{p}
&= \hat{\beta}_{0} \\
&=\bm{g}_{0}^{T} \boldsymbol{y}
= \bm{g}_{0} \cdot \boldsymbol{y}
\end{aligned}
したがって，$X$だけで決まる（データ$\boldsymbol{y}$に依存しない！）$\bm{g}_{0}$がわかれば，各データ$\boldsymbol{y}$に$\bm{g}_{0}$をかけるだけでフィルタリング結果が得られる．以下，$\bm{g}_{0}$を求める方法を考える．
$G = X(X^{T} X)^{-1} $に対して
X^{T} G = I
\end{aligned}
& (X^{T}\bm{g}_{0}, \ldots, X^{T}\bm{g}_{M - 1})
= (\bm{e}_{0}, \ldots, \bm{e}_{M - 1}) \\
&\Rightarrow X^{T} \bm{g}_{0} = \bm{e}_{0}
\end{aligned}
$k$回微分
以下の内容は，$k = 0$とすれば上の内容も含まれる．モデルは
y = \sum_{j = 0}^{N} \beta_{j} x^{j}
\end{aligned}
\frac{\mathrm{d}^{k} y}{\mathrm{d} x^{k}} (x)
&= \sum_{j = 0}^{N} \beta_{j} \gamma_{j} x^{j - k} \\
\gamma_{j}
&=
\begin{cases}
\, 0 & (j < k) \\
\, k! & (j = k) \\
\, \frac{j!}{(j - k)!} & (j > k)
\end{cases}
\end{aligned}
よって，パラメータ$\boldsymbol{\beta}$として上で導いた$\hat{\boldsymbol{\beta}}$を使って導関数を推定することにすれば，フィルタリング結果は$x = 0$における値なので$x^{0} = 1$の項だけが残り
\frac{\mathrm{d}^{k} \hat{y}_{p}}{\mathrm{d} x^{k}}
&= k! \beta_{k} \\
&= k! \bm{g}_{k}^{T} \bm{y} \\
&=k! \bm{g}_{k} \cdot \bm{y}
\end{aligned}
したがって，$X$だけで決まる（データ$\boldsymbol{y}$に依存しない！）$\bm{g}_{k}$がわかれば，各データ$\boldsymbol{y}$に$k! \bm{g}_{k}$をかけるだけでフィルタリング結果が得られる．以下，$\bm{g}_{k}$（$k! \bm{g}_{k}$）を求める方法を考える．
$G = X(X^{T} X)^{-1} $に対して
X^{T} G = I
\end{aligned}
& (X^{T}\bm{g}_{0}, \ldots, X^{T}\bm{g}_{M - 1})
= (\bm{e}_{0}, \ldots, \bm{e}_{M - 1}) \\
&\Leftrightarrow X^{T}\bm{g}_{k} = \bm{e}_{k} \quad (k = 0,1,...,M - 1) \\
&\Leftrightarrow X^{T}(k! \bm{g}_{k}) = k! \bm{e}_{k} \quad (k = 0,1,...,M - 1)
\end{aligned}
