【C++】【Eigen】Savitzky-Golay法を実装する

2021-03-082024-04-16C++,Eigen,スペクトル,信号処理,数学

スペクトル等のノイズを含むデータをスムージングするために広く使われる手法の1つにSavitzky-Golay法があります。このSavitzky-Golay法をC++とEigenを使って実装してみたので紹介します。

2021.03.08　初版公開
2021.10.10　リライト
2024.01.25　記事に間違いがありましたので修正しました

1. Savitzky-Golay法とは
2. C++とEigenによる実装
3. スムージングの例
4. まとめ

Savitzky-Golay法とは

1964年にSavitzkyとGolayが発表したデータのスムージング手法[1]です。隣接する$2N+1$点のデータに対して、$m$次多項式による近似を行い、$2N+1$点の中央の値を近似した多項式の値に置き換えていくことで、ノイズを低減させていきます。
この手法は、実質的に重み付き移動平均であり、[1]の論文にスムージングのために必要な係数が与えられています。計算量も少なく手軽に使えるので便利な手法です（ただしオリジナルの係数表に一部間違いがあります[2]）。

分析化学では、スペクトル等の解析においてノイズを含む生データのスムージングで広く使われているので、名前を聞いたことのある方も多いと思います。

なお、Savitzky-Golay法と比べあまり知られていない、Whittaker Smootherなるスムージング手法もあります。以下のエントリにて紹介していますので、興味のある方はこちらもご覧ください。

【C++】【Eigen】Whittaker Smootherを実装する

スペクトル等のノイズを含むデータをスムージングするために広く使われる手法の1つにSavitzky-Golay法 ...

https://izadori.net/math-whittaker/

特徴

多項式近似にもかかわらず、近似に使うデータの点数と近似する多項式の次数だけで係数が決まる
スムージングだけでなく、平滑化微分も計算できる（3次以上の多項式で近似する場合は2次微分以上も可能）
スムージングの手法にもかかわらず、ピークが鈍りにくい
実質は重み付き移動平均なので、処理が高速

原理

最初に、等間隔に並んだデータについて、$2N+1$点を$m$次多項式$f(x)=a_{0}+a_{1}x+\cdots+a_{m}x^{m}$で近似することを考えます。

このときデータは$\Delta{x}$で等間隔に並んでいるので、$x$の値は$-N\Delta{x}$から$N\Delta{x}$と置き換えることができます。

すると、$(x, y)$のデータの組$(-N\Delta{x}, y_{-N}), \cdots, (0, y_{0}), \cdots, (N\Delta{x}, y_{N})$を$m$次多項式で近似すると考えればよく、その時の最小二乗法、

$$
\displaystyle
S = \sum_{i=-N}^{N}(y_{i} – f(i\Delta{x}))^{2}
$$

の解は、

$$
\displaystyle
\begin{pmatrix}
1 & -N & \cdots & (-N)^{m} \\
\vdots & \vdots & & \vdots \\
1 & 0 & \ddots & 0 \\
\vdots & \vdots & & \vdots \\
1 & N & \cdots & N^{m}
\end{pmatrix}
\begin{pmatrix}
a_{0} \\
\vdots \\
a_{m}\Delta{x}^{m}
\end{pmatrix} =
\begin{pmatrix}
y_{-N} \\
\vdots \\
y_{0} \\
\vdots \\
y_{N} \\
\end{pmatrix}
$$

を解くことで得られることになります。左の行列を$X$と置くと、近似多項式の係数$(a_{0}, \cdots, a_{m}\Delta{x}^{m})$は、

$$
\displaystyle
\begin{pmatrix}
a_{0} \\
\vdots \\
a_{m}\Delta{x}^{m}
\end{pmatrix} =
(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}
\begin{pmatrix}
y_{-N} \\
\vdots \\
y_{0} \\
\vdots \\
y_{N} \\
\end{pmatrix}
$$

で与えられます。ここからがSavitzky-Golay法のミソで、実際の処理で近似多項式の値に置き換えられるのは$2N+1$点の中央の値、すなわち$y_{0}$であり、これに対応する$f(x)$の値は$x=0$を代入した$f(0)=a_{0}$だという点です。そして、この$a_{0}$は、$(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}$の最初の行とベクトル$(y_{-N}, \cdots, y_{N})$の内積で計算できます。実はこの連立方程式を解く必要はまったくなかったりします。

しかも、$(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}$がデータの点数$N$と近似する多項式の次数$m$だけで構成されているので、元データの値に依存せず、データ点数と多項式の次数が同じならば毎回同じ値が利用できます。
$a_{0}$と同様に$a_{1}, a_{2}, \cdots$も計算できますが、$f'(0)=a_{1}\Delta{x}, f"(0)=2a_{2}\Delta{x}^{2}, \cdots$であることがわかれば、微分したデータも簡単に求めることができます。

ちなみにデータの間隔が不定の場合の計算方法もあるようです[3]。しかしながら、通常スペクトル等のデータは等間隔で並んでいるはずなので、あまり出番はなさそうな気がします（論文未見）。

参考文献

C++とEigenによる実装

原理がわかれば実装は比較的簡単です。行列$X$を用意した後に、$(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}$を作り、単なる平滑化データであれば1行目を取り出すことで、係数が得られます。$k$次微分した平滑化データを求める場合は、$k+1$行目を取り出して、$k!$をかけ、さらに$\Delta{x}^{k}$で割ります。

実装例

係数を計算する部分を関数にしたものを以下に示します。エラーチェックはしていないのでご注意を。コードをGitHubにも公開しましたので、興味のある方はご覧ください。

#include <Eigen/Dense>

// Savitzky-Golay法の係数を求める
//
// window   : データ点数（2N+1点）
// polyorder: 近似多項式の次数
// derive   : 微分の次数
// delta    : データの間隔
const Eigen::RowVectorXd CalcSavGolCoeff(
    const size_t window, const unsigned int polyorder, const unsigned derive, const double delta
)
{
    int n = (int)window / 2;
    double p = 1;

    // 行列Xの生成
    Eigen::VectorXd v = Eigen::VectorXd::LinSpaced(window, -n, n);
    Eigen::MatrixXd x = Eigen::MatrixXd::Ones(window, polyorder + 1);

    for(unsigned int i = 1; i <= polyorder; i++){
      x.col(i) = (x.col(i - 1).array() * v.array()).matrix();
    }

    // (X^T * X)^-1 * X^Tを計算する
    Eigen::MatrixXd coeff_mat = (x.transpose() * x).inverse() * x.transpose();

    // 階乗の計算
    if(derive > 0){
      for(unsigned int i = 1; i <= derive; i++){
        p *= i;
      }
      p /= std::pow(delta, derive)
    }

    return coeff_mat.row(derive) * p;
}

使用例

9点のデータに対して5次関数で平滑化するための係数を算出した例を示します。

int main()
{
    double delta = 0.01;
    auto a0 = CalcSavGolCoeff(9, 5, 0, delta);
    auto a1 = CalcSavGolCoeff(9, 5, 1, delta);
    auto a2 = CalcSavGolCoeff(9, 5, 2, delta);

    std::cout << "CalcSavGolCoeff(9, 5, 0, 0.01) = " << std::endl;
    std::cout << a0 << std::endl << std::endl;

    std::cout << "CalcSavGolCoeff(9, 5, 1, 0.01) = " << std::endl;
    std::cout << a1 << std::endl << std::endl;

    std::cout << "CalcSavGolCoeff(9, 5, 2, 0.01) = " << std::endl;
    std::cout << a2 << std::endl << std::endl;

    return 0;
}

出力結果

CalcSavGolCoeff(9, 5, 0, 0.01) = 
 0.034965 -0.128205 0.0699301  0.314685  0.417249  0.314685 0.0699301 -0.128205  0.034965

CalcSavGolCoeff(9, 5, 1, 0.01) = 
-2.96037  16.0956 -26.4452 -33.5548        0  33.5548  26.4452 -16.0956  2.96037

CalcSavGolCoeff(9, 5, 2, 0.01) = 
-734.266     2162  879.953  -1229.6 -2156.18  -1229.6  879.953     2162 -734.266

動作を確認した環境

OS
1. Windows 10
  1. Microsoft C/C++ Optimizing Compiler 19.28.29337 for x64
  2. MinGW-W64 GCC version 8.1.0
2. Linux
  1. GCC version 9.3.0
Eigen version 3.3.9