信号処理では，時系列に含まれる特定の周波数成分を取り出すために，ローパス (low pass)フィルタとか，ハイパス (high pass)フィルタとか，バンドパス (band pass)フィルタとかを使うことがあります．学生の皆さんには，デジタルフィルタの基礎の理論から実際の使い方まで理解しておいてほしいですが，それを学ぶための手っ取り早い参考書がないとも感じています．今回は，「まずは，デジタルフィルタを使ってみよー」ということで，基礎はすっ飛ばして，○○パスフィルタの使い方を簡単に説明しておきます．ここでは，Rのsignalパッケージを使います．
事前準備としてRにsignalパッケージをインストールしておいてください．インストールのコマンドは，
install.packages('signal')です．パッケージをダウンロードするサイトを選んで，インストールしてください．このパッケージを使うときは，最初に１回，
require(signal)
を実行してください．
ローパス (low pass)，ハイパス (high pass)，バンドパス (band pass)とは
時系列はフーリエ変換を使って周波数成分の和に分解できます．つまり，どんな時系列も，周波数成分の和で表せるということです．フィルタ処理では，時系列を構成する周波数成分に対して，必要な周波数領域の成分はそのまま通過させて，必要ない成分を小さくします．例えば，ローパスフィルタは，下の絵のような感じです．
上の図の真ん中の絵が，どの周波数を通過させて，どの周波数を通過させないのかを表しています．同様の図を，ローパス (low pass)，ハイパス (high pass)，バンドパス (band pass)で書いたものが下の図です．通過するか，しないかの分かれ目の周波数がカットオフ周波数 (コーナー周波数) です．
バターワースフィルタ特性をもつデジタルフィルタ
ローパス，ハイパス，バンドパスを実現するデジタルフィルタを使いたいときは，まず，「必要な周波数特性をもつデジタルフィルタを設計しろ」とか言われます．「設計」って言葉の響きが，難しくて，面倒くさそうと感じさせますが，実際は，Rとか，Pythonを使えば簡単です．ここでは，Ｒを使います．
ローパスフィルタの設計
バターワース (Butterworth)フィルタと呼ばれる周波数特性をもつローパスフィルタのゲインは，ローパスフィルタのカットオフ周波数 (コーナー周波数)を
，フィルタの次数を
として，
と表されます．
ここでは，バターワースフィルタの特性をまねたローパスフィルタを設計します．時系列のサンプリング周波数を，ローパスフィルタのカットオフ周波数 (コーナー周波数)を
，フィルタの次数を
とすれば，ローパスフィルタの設計は次のスクリプトで実現できます．
# パッケージの読み込み require(signal) # サンプリング周波数の設定 (20 Hzの例) f.s <- 20 # カットオフ周波数の設定 (2.5 Hzの例) f.c <- 2.5 # ２次のローパスフィルタの設計 bt <- butter(2, f.c/(f.s/2),type="low")
設計した結果がbtに入ってます．
butter(...)の使いこなしのポイントは，カットオフ周波の設定です．butter(... , f.c, ...)とするのではなく，カットオフ周波数をナイキスト周波数で割って，butter(... , f.c/(f.s/2), ...)とする必要があります．
btの中身を見てみると
$b [1] 0.09763107 0.19526215 0.09763107 $a [1] 1.0000000 -0.9428090 0.3333333 attr(,"class") [1] "Arma"
となってます．なんすかこれと疑問に思い，理解してみてーと，思ったあなたは勉強してみてください．
実際の時系列解析でのデジタルフィルタの設計とは，差分方程式の係数を決定することです．入力信号を (
)，フィルタを通過した信号を
とすれば，今回設定したデジタルフィルタって，要するに，
の係数を求めたってことです．つまり，
btの中身から，デジタルフィルタは
になります．この式では，の最初の２つの値
を決めることができません．しかたないので，適当な値を当てはめておいてください．
ハイパスフィルタの設計
時系列のサンプリング周波数を，ハイパスフィルタのカットオフ周波数 (コーナー周波数)を
，フィルタの次数を
とすれば，ハイパスフィルタの設計は次のスクリプトで実現できます．
# パッケージの読み込み require(signal) # サンプリング周波数の設定 (20 Hzの例) f.s <- 20 # カットオフ周波数の設定 (1.0 Hzの例) f.c <- 1.0 # ２次のローパスフィルタの設計 bt <- butter(2,f.c/(f.s/2),type="high")
設計した結果がbtに入ってます．
バンドパスフィルタの設計
時系列のサンプリング周波数を，バンドパスフィルタのカットオフ周波数 (コーナー周波数)を
，フィルタの次数を
とすれば，バンドパスフィルタの設計は次のスクリプトで実現できます．
# パッケージの読み込み require(signal) # サンプリング周波数の設定 (20 Hzの例) f.s <- 20 # カットオフ周波数の設定 ([1.0 Hz, 2.5 Hz]の例) f.c1 <- 1.0 f.c2 <- 2.5 # ２次のローパスフィルタの設計 bt <- butter(2, c(f.c1/(f.s/2), f.c2/(f.s/2)),type="pass")
設計した結果がbtに入ってます．
フィルタをかける
上では，デジタルフィルタを設計しました．ここでは，時系列にフィルタをかける方法を説明します．以下では，xを時系列とします．
Rのsignalパッケージを使ってローパスフィルタをかけたいときは，filter(...)を使います．例えば以下のようになります．
# サンプリング周波数の設定 (20 Hzの例) f.s <- 20 # カットオフ周波数の設定 (2.5 Hzの例) f.c <- 2.5 # ２次のローパスフィルタの設計 bt <- butter(2,f.c/(f.s/2),type="low") # 時系列データxにフィルタをかけてyに代入 y <- filter(bt,x)
このスクリプトでは，yにフィルタをかけた時系列が入っています．
ここまで，フィルタの特性としてゲインの話しかしませんでした．しかし，フィルタの特性には，周波数成分の振幅を変えるだけでなく，位相をずらす性質もあります．ここでは，詳しく説明しませんが，今回設計したデジタルフィルタをかけると，周波数に依存して，位相がずれてしまいます．
そのような位相のずれを小さくしたいときは，filter(...)の代わりにfiltfilt(...)を使います．つまり，
# 時系列データxにフィルタをかけてyに代入 y <- filtfilt(bt,x)
とします．
filtfilt(...)では，フィルタを時系列の前からかけて，その後に後ろからもかけます．こうすることで，フィルタを前からかけたときの位相のずれが，さらに，フィルタを後ろからかけることで元に戻されます．
数値実験
サンプル時系列として白色ノイズを生成してローパスフィルタをかけた結果が下の図です．
上の図では，出力信号の周波数成分の振幅を，入力信号の周波数成分の振幅で割って，フィルタの特性と合っているか検証しています．上図下段の赤破線が，デジタルフィルタの係数から計算したゲイン特性です．実際の振幅比 (ピンク実線)とよく一致しています．図中の青破線は，バターワースフィルタの理論的な特性です．ちょっとずれていることが分かります．
この数値実験に使ったスクリプトはこれです．
# 白色ノイズのサンプル時系列
N <- 10001
x <- rnorm(N)
################################
# デジタルフィルタの設計
# フィルタの次数
q <- 1
# サンプリング周波数 (1 Hzの例)
f.s <- 1
# カットオフ周波数 (0.05 Hzの例)
f.c <- 0.05
# ローパスフィルタの設計
bt <- butter(q,f.c/(f.s/2),type="low")
################################
# フィルタをかける
y <- filter(bt,x)
#################################
# 理論値の計算
###
# デジタルフィルタのゲイン
ord <- length(bt$a)
n.freq <- 500
freq <- seq(0,f.s/2,length.out=n.freq)
G.d <- c()
for(i in 1:n.freq){
w <- exp(-2*pi*1i*(0:(ord-1))*freq[i]/f.s)
H <- (bt$b %*% w)[1,1]/(bt$a %*% w)[1,1]
G.d[i] <- abs(H)
}
#################################
# グラフ描画
par(mfrow=c(2,2),mar=c(5,5,2,2),cex.lab=1.6,cex.axis=1.3)
# 時系列
plot(1:N,x,"l",xaxs="i",ylim=c(-4,4),col=4,xlab="i")
plot(1:N,y,"l",xaxs="i",ylim=c(-4,4),col="lightpink",xlab="i")
# xとyの振幅比
f <- 0:((N-1)/2)/N # 周波数の定義
gain <- sgolayfilt(abs(fft(y)/fft(x))[1:((N+1)/2)],0,3) # 振幅比を３点移動平均で平滑化
plot(f,gain,"l",xlab="f",ylab="Gain",las=1,xaxs="i",yaxs="i",col="lightpink",lwd=4,ylim=c(0,1.2),xlim=c(0,f.s/2))
# デジタルフィルタ特性
lines(freq,G.d,col="red",lwd=2,lty=2)
# Butterworthフィルタ特性
curve(1/sqrt(1+(x/f.c)^(2*q)),add=TRUE,col="blue",lty=2,lwd=2,xlim=c(0,f.s/2))
#
legend("topright",legend=c("fft(y)/fft(x)","digital filter","Butterworth"),col=c("lightpink","red","blue"),cex=1.2,lwd=c(4,2,2),lty=c(1,2,2))
###
# 対数目盛
plot(f[f > 0 & gain > 0],gain[f > 0 & gain > 0],"l",log="xy",xlab="f",ylab="Gain",las=1,xaxs="i",yaxs="i",col="lightpink",lwd=4,ylim=c(1e-2,1.2),xlim=c(4/N,f.s/2),xaxt="n",yaxt="n")
# デジタルフィルタ特性
lines(freq,G.d,col="red",lwd=2,lty=2)
# Butterworthフィルタ特性
curve(1/sqrt(1+(x/f.c)^(2*q)),add=TRUE,col="blue",lty=2,lwd=2,xlim=c(4/N,f.s/2))
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-4:1)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-4:1,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-4:1),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########
legend("bottomleft",legend=c("fft(y)/fft(x)","digital filter","Butterworth"),col=c("lightpink","red","blue"),cex=1.2,lwd=c(4,2,2),lty=c(1,2,2))