はじめに

先日離散データを微分する機会があり，差分法を採用するとデータ点が細かいほど誤差が大きくなるというジレンマに陥りました（一般にデータ点は多いほど精度の高い解析ができるはずです）．

つまり， $x$ は等間隔として，

\frac{f(x_{i+1})+\epsilon_{i+1}-(f(x_{i-1})+\epsilon_{i-1})}{x_{i+1}-x_{i-1}}=\frac{f(x_{i+1})-f(x_{i-1})}{2\Delta x}+\frac{\epsilon_{i+1}-\epsilon_{i-1}}{2\Delta x}（中心差分）

の$\Delta x$が小さいほどその逆数に比例して誤差の影響が大きくなってしまうのです．

この問題への対処法は何パターンもありそうですが，一般にどんな方法が取られるのか調べたところ，SG（Savitzky-Golay）法というものを発見しました．

SG法とは，重み付き移動平均の一種で，
一様な移動平均を各点の周囲 $2N+1$ 個のデータ点の0次近似と見做し，それをより高次に一般化したような方法です．

SG法の原理

$x$ が等間隔なとき，対応する $2N+1$ 個の点 $y_{i-N},...,y_i,...y_{i+N}$ を $m$ 次多項式で近似しすることを考えます．
そして，得られた近似曲線 $f(x)$ に $x_i$ を代入した値 $f(x_i)$ を平滑化後の $y_i$ として採用するのです．

↑ 11点のデータを3次多項式で近似し，中心の点を採用するイメージ．これを各データ点の周りで行う．

ところで，$x_i$ は元の $y_i$ に対応していた値を使うべきでしょうか？
... 答えは否で，$ \Delta x$ おきの等間隔であれば目的は達し得ます．
そこで，$x$ は $0$ を中心とした $-N \Delta x,...,0,...,N \Delta x$ で固定してしまいます．
理由はすぐ後でわかりますが，ここがポイントです．

$(-N \Delta x, y_{i-N}), ..., (N \Delta x, y_{i+N})$ の $m$ 次多項式フィッティングの係数 $\boldsymbol{A} = (a_0,...,a_m)^T$ は，最小二乗法により，

{\left(
\begin{array}\\
1 & -N \Delta x & \cdots & (-N \Delta x)^{m} \\
\vdots &\vdots &\ddots &\vdots\\
1 & 0 &\cdots & 0 \\
\vdots &\vdots &\ddots &\vdots\\
 1 & N \Delta x & \cdots & (N \Delta x)^{m}\\
\end{array}
\right)
\left(
\begin{array}\\
 a_0\\
 a_1\\
 \vdots \\
 a_{m}
\end{array}
\right)
=
\left(
\begin{array}\\
 y_{i-N} \\
 \vdots\\
 y_{i}\\
 \vdots\\
 y_{i+N}\\
\end{array}
\right)
}

\begin{align}
 \boldsymbol{X} \boldsymbol{A} &= \boldsymbol{Y}\\
\boldsymbol{A} &= (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{Y}
\end{align}

により求まり，知りたかった $y_i$ は単に $y_i = a_0$ となっていることが分かります．
これが，$y_i$ に対応する $x$ を $0$ になるよう設定したことによる恩恵です．

更に，$m$ 次多項式 $f(x)=\sum_n a_n x^n $ の $x=0$ における微分値は $\frac{df(x)}{dx}|_{x=0}=\sum_n a_nn 0^{n-1}=a_1$ であるので，上で求めた $\boldsymbol{A}$ からは微分を得ることもできます．
より高階の微分も同様です．

$(\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T$ は使用するデータ点数 $2N+1$ と次数 $m$，$x$ の間隔 $\Delta x$ のみ依存することから，ただ一度だけ計算すれば良く，
得られた $(\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T$ の１行目と $y_{i-N},...,y_i,...y_{i+N}$ の内積を計算すれば $a_0$ は求まります．
$(\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T$ の各行を個別に使うのであれば，それは最早1次元ベクトルなので，フィルターと見做して $\boldsymbol{y}$ に対して畳み込み演算を行えばデータ全体に適用できます．

（$x$ を $-N,...,0,...,N$ としても平滑化までは同じ結果になりますが，$d$ 階微分は $\left(\frac{1}{\Delta x}\right)^d$ 倍になります．）

Pythonによる実装

実装は恐ろしく簡単です．
入力値 $y$ に対し， $y_i$ の前後 $N$ 個のデータ点 $y_{i-N},...,y_i,...y_{i+N}$ を $m$ 次多項式で近似した際の $\left(\frac{d}{dx}\right)^dy$ の $1/d!$ 倍は

def SG(x, y, N, m, d=0):
    dx = x[1] - x[0]
    X = (np.c_[-N:N+1] * dx) ** np.r_[:m+1]
    C = np.linalg.pinv(X) # (X.T X)^-1 X.T
    x_ = x[N:-N]
    y_ = np.convolve(y[::-1], C[d], 'valid')[::-1]
    return x_, y_

で求まります．$x$ の範囲が両側 $N$ ずつ狭くなるので，一緒に返すようにしました．
$d=0$ は平滑化のみです．

使ってみる

試しに適当に正規分布に従うノイズの乗った関数を用意します．$\Delta x$ は $0.01$ にしました．

import numpy as np
import matplotlib.pyplot as plt

dx = 0.01
x = np.r_[-1:1+dx:dx]
coef = np.random.randn(5)
y = np.poly1d(coef)(x) + np.sin(abs(coef).sum() * x)
y_ = y + 0.1 * np.random.randn(*x.shape)

plt.plot(x, y, zorder=10, label='without noise')
plt.plot(x, y_, label='with noise')
plt.legend()
plt.show()

そのまま差分法で微分値を求めてみます．

plt.plot(x, np.gradient(y, x), zorder=10, label='without noise')
plt.plot(x, np.gradient(y_, x), label='with noise')
plt.legend()
plt.show()

ノイズは $1/\Delta x = 100$ 倍となり，ひどい有様です．

ここにSG法で平滑化したものを追加してみます．$N=5$, $m=2$ としてみました．

plt.plot(x, np.gradient(y, x), zorder=10, label='without noise')
plt.plot(x, np.gradient(y_, x), label='with noise')
plt.plot(*SG(x, y_, 10, 2, d=1), label='with noise + SG(N=5, m=2)')
plt.legend()
plt.show()

かなりノイズが抑え込めていることが分かります．

フィルターサイズ $2N+1$ を大きくするほど，次数 $m$ を小さくするほど強く平滑化されます．用途に合わせて調節しましょう．

おまけ

Whittaker Smoother も非常に簡単に実装できます．

def WS(y, n, lmd):
    m = len(y)
    D = np.diff(np.eye(m + n), n)[n:-n]
    W = np.diag(np.ones(m))
    return np.linalg.pinv(W + lmd * D.T @ D) @ W @ y

おわりに

測定データを微分することってあまりないのでしょうか？
あまりノイズが大きくなる問題に言及している記事を見かけません．
Python派の方の参考になると幸いです．

参考

Abraham. Savitzky and M. J. E. Golay, Smoothing And differentiation of Data by Simplified Least Squares Procedures. Anal. Chem. 36(8), 1627-1639 (1964).
【C++】【Eigen】Savitzky-Golay法を実装する | いざどりのtrial and error