ケィオスの時系列解析メモランダム

時系列解析，生体情報学，数学・物理などの解説です．

【Rで時系列解析】Savitzky-Golayフィルタで平滑化

時系列解析 R言語

今回はギザギザ，凸凹したデータの平滑化法の一つであるSavitzky-Golayフィルタの紹介です．

Savitzky-Golayフィルタの考え方 (2次でスケール5の例)：部分区間に多項式をフィットして中央の値を平滑化された値として採用する．実用上は，畳み込み和を使って実装されることがほとんどで，最小2乗多項式を計算する必要はない．

この前，単純な移動平均を使った平滑化を拡張する方法として，修正ダニエルスムーザとその畳込み（複数回の修正ダニエルスムーザ）を紹介しました．
chaos-kiyono.hatenablog.com

　この前は，移動平均を使った平滑化の名前を「ダニエルスムーザ」と紹介しました．同じものでも，呼び方は業界によって異なます．今回は，移動平均を使った平滑化の名前を「0次Savitzky-Golayフィルタ (Savitzky-Golayスムーザ)」と呼びます．"Savitzky-Golay"のカタカナ読みについては，自信がありませんが，サヴィツキー-ゴレイと私は読んでいます．

　式で説明する前に，実際に平滑化している様子を図に示します．注意してほしいことは，実用上は，畳込み演算 (重み付き和)でフィルタは計算され，曲線をフィットしたりしません．

Savitzky-Golayフィルタを使った平滑化

1. Savitzky-Golayフィルタの考え方

　下の図は，2次のSavitzky-Golayフィルタで，スケール (幅の点数) 5の例を使って手順 (考え方)を描いています．

2次Savitzky-Golayフィルタの適用例 (スケール5)．

$q$ 次Savitzky-Golayフィルタの考え方は以下です．

奇数の長さのスケール $s$ (幅の点数) を設定し，その区間に $q$ 次多項式を最小2乗法でフィットします．
フィットした多項式の区間における中央の点を，平滑化された値として採用します．
多項式をフィットする区間を一つずつ右にシフトして，上の手順を繰り返す．

2. Savitzky-Golayフィルタの畳込み表現

　時系列を $y[i$ ]としたとき，スケール $s=2m+1$ のSavitzky-Golayフィルタは，

　 $\displaystyle \hat{y}[i] = \sum_{j=-m}^{m} w_{j}\, y[i-j]$

の形で表すことができます．興味があれば， $w_{j}$ を導いてみてください．下の図は， $w_{j}$ を描いたものです．

Savitzky-Golayフィルタの重み $w_j$ ．

w_j — Savitzky-Golayフィルタの重み $w_j$ ．

　さらに，高速な計算アルゴリズムも存在します．
Phys. Rev. E 93, 053304 (2016) - Fast algorithm for scaling analysis with higher-order detrending moving average method

3. Savitzky-Golayフィルタの特徴

　Savitzky-Golayフィルタの特徴として私が強調したいことは，

$q$ を奇数として， $q$ 次Savitzky-Golayフィルタは， $(q+1)$ 次多項式で記述される曲線を再現できる．次数が1大きいぞ！

ということです．偶数次のSavitzky-Golayフィルタは，1だけ小さい奇数次のSavitzky-Golayフィルタと同じ結果になります．

　他に，FIRフィルタなので，位相遅れがないとか，周波数特性とかありますが，私がSavitzky-Golayフィルタを使う主な理由はこれです．

　ここでは，数式の証明は打つのが面倒なので与えませんが，数値実験の結果を下に載せておきます．

0次Savitzky-Golayフィルタは，直線を再現できる．

0次Savitzky-Golayフィルタを直線に適用した結果．

0次Savitzky-Golayフィルタは，2次関数を再現できない．

0次Savitzky-Golayフィルタを2次関数に適用した結果．ずれてます．

2次Savitzky-Golayフィルタは，直線を再現できる．

2次Savitzky-Golayフィルタを直線に適用した結果．

2次Savitzky-Golayフィルタは，2次関数を再現できる．

2次Savitzky-Golayフィルタを2次関数に適用した結果．

2次Savitzky-Golayフィルタは，3次関数を再現できる．

2次Savitzky-Golayフィルタを3次関数に適用した結果．

2次Savitzky-Golayフィルタは，4次関数を再現できない．

2次Savitzky-Golayフィルタを4次関数に適用した結果．わずかですがずれてます．

4. Rで実行

　Rでは，signalパッケージを導入すれば，Savitzky-Golayフィルタをかけることができます．コマンドは，

sgolayfilt(x, p, n, m = 0, ts = 1)

です．詳しい使い方の説明は次回にします．