【信号処理の基礎】Savitzky-Golay平滑化・微分フィルタの周波数特性

今回は，Savitzky-Golayフィルタの畳込み表現の係数 (下図のピンク図形)と周波数特性の話です．

Savitzky-Golay平滑化フィルタ．元の時系列（上段灰色）は，下段の破線にノイズを加えたもの．下段の灰色破線はノイズを加える前の時系列．赤実線（上下両方）は，元時系列（上段灰色）にSavitzky-Golay平滑化フィルタをかけた結果．ピンクの図形は畳込み係数を表す．

　Savitzky-Golayフィルタは，中央移動平均を一般化したものの一つです．中央移動平均の場合，奇数の長さ $s$ の部分区間に，定数関数

　 $f_n = c_0$

をフィットし，中央の値を平滑化された値として採用します．部分区間は1点ずつシフトしていきます．Savitzky-Golayフィルタでは，当てはめる関数を $q$ 次多項式

　 $f_n = c_0+c_1 n + c_2 n^2 + \cdots + c_m n^q$

に変えてフィットし，中央の値を平滑化された値として採用します．以前にも，この考え方は説明しました．
chaos-kiyono.hatenablog.com

　実際にRを使ってフィルタをかけるとき，私はsignalパッケージのsgolayfiltを使います．

　今回は，Savitzky-Golayフィルタについての理解を深めるために，Savitzky-Golayフィルタを畳込み

　 $y_n = \displaystyle \sum_{k=-\infty}^{\infty} a_k\, x_{n-k}$

で表現したときの係数を与えます．さらに，Savitzky-Golay微分フィルタを紹介します．

Savitzky-Golay平滑化フィルタの係数
Savitzky-Golay微分フィルタの係数
Rスクリプト

Savitzky-Golay平滑化フィルタの係数

　計算を簡単にするために，多項式をフィットする部分区間を $\displaystyle \left[ -\frac{s-1}{2}, \frac{s-1}{2} \right]$ にします．部分区間の場所が変っても畳込みの係数は変りません．

　 $\displaystyle m = \frac{s-1}{2}$ とおくと，考える部分時系列は $\left\{ x_{-m}, x_{-m+1}, \cdots, 0, x_{m-1}, x_{m} \right\}$ です．

　今回は２次のSavitzky-Golayフィルタを考えます．部分時系列 $\left\{ x_{-m}, x_{-m+1}, \cdots, 0, x_{m-1}, x_{m} \right\}$ に最小２乗法で２次関数を当てはめます．

　 $\displaystyle I \left( \left\{c_i \right\} \right) = \sum_{k = -m}^{m} \left\{ x_{k} - \left(c_0 + c_1 k + c_2 k^2 \right) \right\}^2$

を最小にするのは，

　 $\displaystyle \frac{\partial I \left( \left\{c_i \right\} \right)}{\partial c_0} = 0, \quad \frac{\partial I \left( \left\{c_i \right\} \right)}{\partial c_1} = 0, \quad \frac{\partial I \left( \left\{c_i \right\} \right)}{\partial c_2} = 0$

を満たす $\{c_0, c_1, c_2\}$ です．上の条件を計算して，行列にまとめれば，

　 $\begin{eqnarray} \begin{bmatrix} \displaystyle \sum_{k=-m} ^{m} 1 & 0 & \displaystyle \sum_{k=-m} ^{m} k^2 \\ 0 &\displaystyle \sum_{k=-m} ^{m} k^2 & 0 \\ \displaystyle\sum_{k=-m} ^{m} k^2 & 0 & \displaystyle \sum_{k=-m} ^{m} k^4 \\ \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \end{bmatrix} &=& \begin{bmatrix} \displaystyle \sum_{k=-m} ^{m} x_k \\ \displaystyle \sum_{k=-m} ^{m} k \, x_k \\ \displaystyle \sum_{k=-m} ^{m} k^2 \, x_k \\ \end{bmatrix} \\ \begin{bmatrix} s & 0 & \frac{s^3-s}{12} \\ 0 & \frac{s^3-s}{12} & 0 \\ \frac{s^3-s}{12} & 0 & \frac{3 s^5-10 s^3+7 s}{240} \\ \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \end{bmatrix} &=& \begin{bmatrix} \displaystyle \sum_{k=-m} ^{m} x_k \\ \displaystyle \sum_{k=-m} ^{m} k \, x_k \\ \displaystyle \sum_{k=-m} ^{m} k^2 \, x_k \\ \end{bmatrix} \end{eqnarray}$

これを， $\{c_0, c_1, c_2\}$ について解けば，

　 $\begin{eqnarray} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \end{bmatrix} &=& \begin{bmatrix} \frac{21-9 s^2}{16 s-4 s^3} & 0 & \frac{15}{4 s-s^3} \\ 0 & -\frac{12}{s-s^3} & 0 \\ \frac{15}{4 s-s^3} & 0 & \frac{180}{s^5-5 s^3+4 s} \\ \end{bmatrix}\begin{bmatrix} \displaystyle \sum_{k=-m} ^{m} x_k \\ \displaystyle \sum_{k=-m} ^{m} k \, x_k \\ \displaystyle \sum_{k=-m} ^{m}k^2 \, x_k \\ \end{bmatrix}\\ \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \end{bmatrix} &=& \begin{bmatrix} \displaystyle \frac{21-9 s^2}{16 s-4 s^3} \sum_{k=-m} ^{m} x_n + \frac{15}{4 s-s^3} \sum_{k=-m} ^{m} k^2 \, x_k \\ \displaystyle -\frac{12}{s-s^3} \sum_{k=-m} ^{m} n \, x_k \\ \displaystyle \frac{15}{4 s-s^3} \sum_{k=-m} ^{m} x_k + \frac{180}{s^5-5 s^3+4 s} \sum_{k=-m} ^{m} k^2 \, x_k \\ \end{bmatrix}\\ \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \end{bmatrix} &=& \begin{bmatrix} \displaystyle \sum_{k=-m} ^{m} \left( \frac{9 s^2 - 21}{4 s^3 - 16 s} - \frac{15 k^2}{s^3-4 s} \right) x_k \\ \displaystyle \sum_{k=-m} ^{m} \frac{12 k}{s^3-s}\, x_n \\ \displaystyle \sum_{k=-m} ^{m} \left( \frac{180 k^2}{s^5-5 s^3+4 s} - \frac{15}{s^3-4 s}\right) x_k \\ \end{bmatrix} \end{eqnarray}$

　部分区間 $\displaystyle \left[ -\frac{s-1}{2}, \frac{s-1}{2} \right]$ の中央点 $n=0$ のみを採用するので，

　 $\begin{eqnarray}y_0 &=& c_0 + c_1 \cdot 0 + c_2 \cdot 0^2 = c_0 \\ &=& \sum_{k=-m} ^{m} \left( \frac{9 s^2 - 21}{4 s^3 - 16 s} - \frac{15 k^2}{s^3-4 s} \right) x_k \\ &=& \sum_{k=-m} ^{m} \frac{9 s^2 - 21 - 60 k^2}{4 s^3 - 16 s} x_k \end{eqnarray}$

になります．

　ということで，畳込み用に係数の順番をひっくり返しても同じなので，

幅 $s$ の２次Savitzky-Golay平滑化フィルタの係数は，
　 $\displaystyle a_k = \frac{9 s^2 - 21 - 60 k^2}{4 s^3 - 16 s}$

です．

Savitzky-Golay微分フィルタの係数

　Savitzky-Golayフィルタでは，平滑化した時系列の微分係数を求めることも可能です．

　２次Savitzky-Golayフィルタで部分区間にフィットした関数

　 $f_k = c_0 + c_1 k + c_2 k^2$

を微分し， $k=0$ を代入すると，

　 $\displaystyle \left.\frac{d f_k}{d k} \right|_{k=0} = c_1 = \displaystyle \sum_{k=-m} ^{m} \frac{12 n}{s^3-s}\, x_n$

です．

　ということで，畳込み用に係数の順番をひっくり返して，

幅 $s$ の２次Savitzky-Golay微分フィルタの係数は，
　 $\displaystyle a_k = -\frac{12 k}{s^3-s}$

です．

Savitzky-Golay平滑化フィルタ（上）と微分フィルタ（下）の係数（左）と周波数特性（右）．フィルタ幅 $s=9$ ．係数は，畳込みの場合と順番が逆になっています．

Savitzky-Golay平滑化フィルタ（上）と微分フィルタ（下）の係数（左）と周波数特性（右）．フィルタ幅 $s=9$ ．係数は，畳込みの場合と順番が逆になっています．

Rスクリプト

2次Savitzky-Golay平滑化フィルタの係数．

sg <- function(k,s){
 return((9*s^2-21-60*k^2)/(4*s^3-16*s))
}
# 幅
s <- 9
k <- (-(s-1)/2):((s-1)/2)
plot(k,sg(k,s),col=2,ylab="a_k",main="SG smoothing filter")

2次Savitzky-Golay微分フィルタの係数．

dsg <- function(k,s){
 return(-(12*k)/(s^3-s))
}
# 幅
s <- 9
k <- (-(s-1)/2):((s-1)/2)
plot(k,dsg(k,s),col=2,ylab="a_k",main="SG differentiation filter")

周波数特性もプロット

### 平滑化フィルタの係数
sg <- function(k,s){
 return((9*s^2-21-60*k^2)/(4*s^3-16*s))
}
###　微分フィルタの係数
dsg <- function(k,s){
 return((12*k)/(s^3-s))
}
########################
# フィルタ幅
s <- 9
########################
k <- (-(s-1)/2):((s-1)/2)
# 周波数
f.list <- exp(seq(log(0.001),log(0.5),length.out=1000))
# インパルス応答の計算
Freq.Resp <- c()
for(f in f.list){
 vcos <- cos(2*pi*f*k) 
 Freq.Resp <- c(Freq.Resp,sg(k,s) %*% vcos)
}
# plot
par(mfrow=c(2,2))
plot(k,sg(k,s),col=2,ylab="a_k",main="SG smoothing filter")
y.max <- max(Freq.Resp)
plot(f.list,abs(Freq.Resp),"l",col=2,log="xy",xlim=c(0.001,0.5),ylim=c(y.max/100,y.max),xaxs="i",lwd=2,xlab="f [Hz]",ylab="|A(f)|",main="SG smoothing filter")
################################
# インパルス応答の計算
Freq.Resp <- c()
for(f in f.list){
 vcos <- sin(2*pi*f*k) 
 Freq.Resp <- c(Freq.Resp,dsg(k,s) %*% vcos)
}
plot(k,dsg(k,s),col=2,ylab="a_k",main="SG differentiation filter")
y.max <- max(Freq.Resp)
plot(f.list,abs(Freq.Resp),"l",col=2,log="xy",xlim=c(0.001,0.5),ylim=c(y.max/100,y.max),xaxs="i",lwd=2,xlab="f [Hz]",ylab="|A(f)|",main="SG differentiation filter")