【信号処理の基礎】Savitzky–Golayフィルタを一般形で表す - ケィオスの時系列解析メモランダム

Savitzky–Golay（SG）フィルタは、時系列データをなめらかにしながら波形の形をできるだけ壊さない、非常に便利な平滑化手法です。生体信号解析、化学計測、音声処理、パワースペクトルの平滑化など、幅広い分野で利用されています。

　SG フィルタの特徴は、単なる平均ではなく、多項式をフィットする点にあります。窓の中のデータ点を、次数 $\displaystyle{ q }$ の多項式

$\displaystyle{ f(k)=c_0 + c_1 k + c_2 k^2 + \cdots + c_q k^q }$

で「できるだけよく近似する」ようにし、その多項式の値を平滑化値として使うのです。

　SG フィルタは、窓の中のデータに対し、多項式を最小二乗法でフィットします。最小二乗法は行列でとてもきれいに書けるため、一般次数 $\displaystyle{ q }$ の SG フィルタは、

窓幅 $\displaystyle{s}$
多項式次数 $\displaystyle{ q}$
評価点 $\displaystyle{k _ 0}$
微分次数 $\displaystyle{d}$

といったパラメタを自由に変えても同じ形式で扱えます。

　今回は以下では、SG フィルタを一般形を行列形式で見ていきます。

1. 問題設定：窓幅と多項式の次数
2. 行列で書く：設計行列の導入
3. 最小二乗法：正規方程式と擬似逆行列
4. 「どの点の値をとるか」をベクトルで表す
5. 一般次数の SG フィルタの「行列形式の一般式」
6. Rでフィルタ係数を計算
7. まとめ

1. 問題設定：窓幅と多項式の次数

まず、記号をはっきりさせておきます。

時系列： $\displaystyle{\{x _ n\}}$
フィルタをかけたい点：インデックス $\displaystyle{n}$
窓幅（サンプル数）： $\displaystyle{s = 2m+1}$ （奇数）
窓の中の相対インデックス： $\displaystyle{ k = -m, -m+1, \dots, 0, \dots, m}$
多項式の次数： $\displaystyle{q}$ （0 次〜）

SG フィルタでは、この窓の中のデータ

$\displaystyle{ x_{n-m}, \dots, x_{n}, \dots, x_{n+m} }$

を「相対インデックス $\displaystyle{k}$ 」で表して

$\displaystyle{ x_{n+k} \quad (k=-m,\dots,m) }$

と書き直します。そして、これらに対して、次数 $\displaystyle{ q }$ の多項式

$\displaystyle{ f(k) = c_0 + c_1 k + c_2 k^2 + \cdots + c_q k^q }$

を、最小二乗法でフィットし、その一部（たとえば $\displaystyle{ k=0 }$ の値）を平滑化結果として採用します。

2. 行列で書く：設計行列 $\displaystyle{A}$ の導入

　ここでは，時系列データの中から，解析したい区間を一定長の窓（ウィンドウ）として切り出し，その区間のデータ列（部分時系列）を一つの列ベクトルとして表します。

$\displaystyle{ x = \begin{bmatrix} x_{n-m}\\ x_{n-m+1}\\ \vdots\\ x_{n+m} \end{bmatrix} \in \mathbb{R}^{s} }$

　この部分時系列に当てはめる多項式

$\displaystyle{ f(k) = c_0 + c_1 k + c_2 k^2 + \cdots + c_q k^q }$

の係数も列ベクトルにします。

$\displaystyle{ c = \begin{bmatrix} c_0\\ c_1\\ \vdots\\ c_q \end{bmatrix} \in \mathbb{R}^{q+1} }$

　ここで、設計行列（デザイン行列）を

$\displaystyle{ A = \begin{bmatrix} 1 & (-m) & (-m)^2 & \cdots & (-m)^q \\ 1 & (-m+1) & (-m+1)^2 & \cdots & (-m+1)^q \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & (-1) & (-1)^2 & \cdots & (-1)^q \\ 1 & 0 & 0^2 & \cdots & 0^q \\ 1 & 1 & 1^2 & \cdots & 1^q \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & m & m^2 & \cdots & m^q \end{bmatrix} }$

と定義します。これは「各サンプル $\displaystyle{ k }$ ( $\displaystyle{ k = -m, -m+1, \cdots, m-1, m }$ ) に対して、 $\displaystyle{ 1,k,k^ 2,\dots,k^ q }$ を並べた行列」です。サイズは

行数：窓幅 $\displaystyle{s}$
列数：多項式のパラメータ数 $\displaystyle{q+1}$

なので、 $\displaystyle{ A\in\mathbb{R}^ {s\times (q+1)} }$ です。

　このとき、多項式での近似値の列ベクトル $\displaystyle{ f }$ は

$\displaystyle{ f = \begin{bmatrix} f(-m)\\ f(-m+1)\\ \vdots\\ f(m) \end{bmatrix} = A c }$

と書けます。

3. 最小二乗法：正規方程式と擬似逆行列

　目標は、多項式 $\displaystyle{ f(k) }$ が元データ $\displaystyle{ x _ {n+k} }$ をよく近似するように、誤差の二乗和

$\displaystyle{ I(c) = \sum_{k=-m}^{m} \left\{ x_{n+k} - f(k) \right\}^2 }$

を最小にする、 $\displaystyle{ c }$ を求めることです。

　行列表現では、この誤差は

$\displaystyle{ I(c) = \|x - A c\|^2 }$

です。これを展開して整理すると、

$\displaystyle{ \begin{aligned} I(c) &= \|x - A c\|^2\\ &= (x - A c)^T (x - A c)\\ &= x^T x - x^T A c - c^T A^T x + c^T A^T A c. \end{aligned} }$

となります。

　ここでまず、行列積の転置の公式

$\displaystyle{ (XY)^T = Y^T X^T }$

を思い出してください。 $\displaystyle{ x^ T A c }$ と $\displaystyle{ c^ T A^ T x }$ は転置の関係にあり、どちらもスカラー量なので、

$\displaystyle{ x^T A c = c^T A^T x }$

が成り立ちます。したがって、

$\displaystyle{ \begin{aligned} I(c) &= x^T x - 2\,c^T A^T x + c^T A^T A c \end{aligned} }$

と整理できます。

　最小２乗法では、 $\displaystyle{I(c)}$ を最小にする $\displaystyle{c}$ を求めたいので、 $\displaystyle{I(c)}$ の勾配（ベクトル微分）が 0 になる点 $\displaystyle{c}$ を求めます。

　行列やベクトルの微分を、自信をもって計算できる人は少ないと思いますが、上の式に含まれる各項の微分は次のようになります：

$\displaystyle{\dfrac{\partial}{\partial c}\, x^ T x = 0}$ （ $\displaystyle{x}$ だけの項なので $\displaystyle{c}$ に依存しない）
$\displaystyle{\dfrac{\partial}{\partial c}\, (-2\,c^ T A^ T x) = -2 A^ T x}$
$\displaystyle{\dfrac{\partial}{\partial c}\, (c^ T A^ T A c) = 2 A^ T A c}$

ここで注意しておきたいのは、 $\displaystyle{c}$ がベクトルであるため、通常の変数で記述された関数の微分とは違うということです。とはいえ、二次形式や線形項の微分ルールはスカラーの場合とよく似た形になります。そのため、上記のような計算結果が自然に導かれます。

　結果として、

$\displaystyle{ \frac{\partial I(c)}{\partial c} = 2 A^T A c - 2 A^T x }$

が得られます。最小値ではこの値が 0 になるので，

$\displaystyle{ 2 A^T A c - 2 A^T x = 0 }$

のはずです。この式が、 $\displaystyle{c}$ を求めるための方程式です。この式を少し変形して、

$\displaystyle{ A^T A c = A^T x }$

を得ます。これが、 正規方程式（normal equation） と呼ばれる式です。

　求めた正規方程式が正しいか心配なので、成分で書き下して、確かめてみます。誤差の二乗和は

$\displaystyle{ I(c) = \sum_{k=-m}^{m} \left\{ x_{n+k} - f(k) \right\}^2 = \sum_{k=-m}^{m} \left\{ x_{n+k} - \sum_{j=0}^{q} A_{kj} c_j \right\}^2 }$

です（ここで、 $\displaystyle{ A _ {kj} }$ は $\displaystyle{A}$ の $\displaystyle{k}$ 行 $\displaystyle{j}$ 列成分，例えば $\displaystyle{ A _ {kj} = k^ j }$ です）。

　各係数 $\displaystyle{ c _ \ell }$ について偏微分して 0 とおくと，

$\displaystyle{ \frac{\partial I}{\partial c_\ell} = -2 \sum_{k=-m}^{m} A_{k\ell} \left\{ x_{n+k} - \sum_{j=0}^{q} A_{kj} c_j \right\} = 0 }$

となり，両辺を 2 で割って整理すると

$\displaystyle{ \sum_{k=-m}^{m} \sum_{j=0}^{q} A_{k\ell} A_{kj} c_j = \sum_{k=-m}^{m} A_{k\ell} x_{n+k} }$

です。左辺は $\displaystyle{ (A^ T A c) _ \ell }$ 、右辺は $\displaystyle{ (A^ T x) _ \ell }$ に対応するので，

$\displaystyle{ (A^T A c)_\ell = (A^T x)_\ell \quad (\ell = 0,\dots,q) }$

すなわちベクトル・行列表現で

$\displaystyle{ A^T A c = A^T x }$

と書き直せます。

　ということで、最小二乗法で得られる解は、正規方程式

$\displaystyle{ A^T A\, c = A^T x }$

を満たす $\displaystyle{ c }$ として与えられます。もし $\displaystyle{ A^ T A }$ が正則であれば、

$\displaystyle{ c = (A^T A)^{-1} A^T x }$

です。この

$\displaystyle{ A^\dagger := (A^T A)^{-1} A^T }$

は、線形代数でいう「擬似逆行列（Moore–Penrose の疑似逆）」に相当します。 $\displaystyle{ A^ \dagger }$ を使えば、係数はシンプルな形、

$\displaystyle{ c = A^\dagger x }$

となります。

4. 「どの点の値をとるか」をベクトルで表す

多項式

$\displaystyle{ f(k) }$

の値は、

$\displaystyle{ f(k)= \begin{bmatrix} 1 & k & k^2 & \cdots & k^q \end{bmatrix} c }$

と書けます。そこで、

$\displaystyle{ v(k) = \begin{bmatrix} 1\\ k\\ k^2\\ \vdots\\ k^q \end{bmatrix} }$

とおいて

$\displaystyle{ f(k) = v(k)^T c }$

と書くことにします。

　いま、SG 「平滑化」フィルタでは、通常 $\displaystyle{ k=0 }$ の値を採用しますので、

$\displaystyle{ y_n = f(0) = v(0)^T c }$

です。このとき、 $\displaystyle{ v(0) }$ は

$\displaystyle{ v(0) = \begin{bmatrix} 1\\ 0\\ 0\\ \vdots\\ 0 \end{bmatrix} = e_1 }$

（最初の成分だけが 1 の単位ベクトル）です。

ここで、 $\displaystyle{ c = A^ \dagger x }$ を代入すると

$\displaystyle{ y_n = e_1^T c = e_1^T A^\dagger x }$

となります。つまり、線形写像 $\displaystyle{ x\mapsto y _ n }$ の係数ベクトルは、

$\displaystyle{ h^T = e_1^T A^\dagger }$

です。この $\displaystyle{ h }$ の成分が、窓内のサンプルに対する畳み込み係数になります。

5. 一般次数 $\displaystyle{q}$ の SG フィルタの「行列形式の一般式」

　上の議論をまとめると、「窓幅 $\displaystyle{ s=2m+1 }$ 、多項式次数 $\displaystyle{ q }$ の SG 平滑化フィルタ」は次のように定義できます。

5.1 ステップ 1：設計行列 $\displaystyle{A}$

　相対インデックス $\displaystyle{ k=-m,\dots,m }$ に対して、

$\displaystyle{ A_{(k+m+1), j+1} = k^{j},\quad (j=0,1,\dots,q) }$

と定義します。これで、 $\displaystyle{ A \in \mathbb{R}^ {s\times(q+1)} }$ が定まります。

5.2 ステップ 2：擬似逆行列 $\displaystyle{A^ \dagger}$

$\displaystyle{ A^\dagger = (A^T A)^{-1} A^T }$

とすれば、パラメタ $\displaystyle{ c }$ とデータ $\displaystyle{ x }$ の関係が

$\displaystyle{ c = A^\dagger x }$

と書けます。

5.3 ステップ 3：評価点の指定ベクトル $\displaystyle{v(k _ 0)}$

　評価したい相対位置 $\displaystyle{ k _ 0 }$ （通常は 0）に対して

$\displaystyle{ v(k_0) = \begin{bmatrix} 1\\ k_0\\ k_0^2\\ \vdots\\ k_0^q \end{bmatrix} }$

とします。平滑化で中央を取りたいときは、 $\displaystyle{ k _ 0=0 }$ なので、 $\displaystyle{ v(0)=e _ 1 }$ です。

5.4 ステップ 4：フィルタ係数ベクトル $\displaystyle{h(k _ 0)}$

　出力

$\displaystyle{ y_n = f(k_0) = v(k_0)^T c }$

に、 $\displaystyle{ c = A^ \dagger x }$ を代入すると

$\displaystyle{ y_n = v(k_0)^T A^\dagger x }$

です。

　したがって、窓内のデータ $\displaystyle{ x _ {n-m},\dots,x _ {n+m} }$ に対する畳み込み係数は

$\displaystyle{ h(k_0)^T = v(k_0)^T A^\dagger }$

となります。

　成分で書くと、

$\displaystyle{ y_n = \sum_{r=-m}^{m} h_r(k_0)\, x_{n+r} }$

ここで、ベクトル $\displaystyle{ h(k _ 0)\in\mathbb{R}^ {s} }$ の各要素 $\displaystyle{ h _ r(k _ 0) }$ が、「相対インデックス $\displaystyle{ r }$ のサンプルに対する SG フィルタ係数」ということになります。

これが、次数 $\displaystyle{ q }$ の SG フィルタの一般行列表現です。

6. Rでフィルタ係数を計算

　上の 5 で説明した流れにそってフィルタ係数を求める関数のRスクリプトをのせておきます。

###############################################
# Savitzky–Golay フィルタ係数 h(k0) を求める関数
#  - m  : 窓の片側長 (s = 2m+1)
#  - q  : 多項式次数
#  - k0 : 評価点（相対インデックス，通常は 0）
###############################################

sg_coeff <- function(m, q, k0 = 0) {
  #-------------------------------------------
  # Step 0: 基本量の定義
  #   k = -m, ..., m  （相対インデックス）
  #   s = 2m+1        （窓幅）
  #-------------------------------------------
  k <- -m:m
  s <- length(k)  # = 2m+1
  
  #-------------------------------------------
  # Step 1: 設計行列 A （s × (q+1)）
  #   A_{(k+m+1), j+1} = k^j,  j = 0,...,q
  #   outer() を使うと簡潔に書ける
  #-------------------------------------------
  A <- outer(k, 0:q, function(kk, j) kk^j)
  # 行: k=-m,...,m, 列: 1, k, k^2, ..., k^q
  # この A が本文の A ∈ R^{s×(q+1)} に対応
  
  #-------------------------------------------
  # Step 2: 擬似逆行列 A^dagger = (A^T A)^{-1} A^T
  #-------------------------------------------
  AtA    <- t(A) %*% A
  A_pinv <- solve(AtA) %*% t(A)  # A^dagger
  
  #-------------------------------------------
  # Step 3: 評価点ベクトル v(k0)
  #   v(k0) = [1, k0, k0^2, ..., k0^q]^T
  #-------------------------------------------
  v <- k0^(0:q)           # 長さ q+1 のベクトル
  v <- matrix(v, ncol = 1)  # 列ベクトルとして扱う
  
  #-------------------------------------------
  # Step 4: フィルタ係数ベクトル h(k0)
  #   h(k0)^T = v(k0)^T A^dagger
  #   → h は長さ s (=2m+1) のベクトル
  #-------------------------------------------
  h_row <- t(v) %*% A_pinv     # 1 × s 行ベクトル
  h     <- as.numeric(h_row)   # 普通の数値ベクトルに変換
  
  # 相対インデックス r = -m,...,m を名前として付与（任意）
  names(h) <- paste0("k=", k)
  
  return(h)
}

　フィルタ係数を求める場合は、以下の命令を実行してください。

# 窓幅 s = 2m+1 = 5, 多項式次数 q = 2, 中央値 (k0=0) の SG フィルタ係数の例
sg_coeff(m = 2, q = 2, k0 = 0)

7. まとめ

　SG フィルタの一般形を行列で見直すことは、単に一つのフィルタを理解するだけでなく、線形代数そのものの威力を実感する良い機会だと思います。設計行列 $\displaystyle{ A }$ 、擬似逆行列 $\displaystyle{ A^ \dagger }$ 、直交射影 $\displaystyle{ P = AA^ \dagger }$ といった道具立ては、今回の SG フィルタだけでなく、最小二乗回帰、主成分分析、一般的な FIR フィルタ設計、さらには機械学習の多くのアルゴリズムにも共通して現れるものです。

　こうした線形代数の視点をしっかり身につけておくことで、「なぜこのフィルタがうまく働くのか」、「どこまで一般化できるのか」を、公式暗記ではなく構造として理解できるようになります。SG フィルタを入り口にして、線形代数・最小二乗・射影といった概念を、学んでみてください。

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。