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

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

【時系列解析】フーリエ変換を使った時系列の分解

時系列解析で,弱定常とか,ガウス過程とかを前提とする場合は,パワースペクトルがとにかく基本.それを理解せずに,なんとなく数式をこねくり回したり,名前が格好いい解析法を中身を理解せずに使ったりしても,なんか学問として浅い感じがします.といっても道のりは意外と長いので,まずは,パワースペクトルを理解する前段階として,フーリエ変換を使って時系列を周期成分に分解することについて感じてもらいます.

Rで実感する時系列の分解

フーリエ変換が存在する条件は気にしなくても,我々が解析する時系列は,離散で,実数の有限の値をとり,長さも有限なので,絶対,離散フーリエ変換できます.ただし,欠損値がない場合です.

理屈は抜きにして,「周波数成分への分解」とか,「周波数成分の和で表せる」とかを体験するために,以下のRスクリプトを実行してみてください.

# 時系列の設定 (自由に決めてね)
# c(...)の中にコンマ区切りで書く.長くするとグラフの数が多くなるぞ.
x <- c(-1,1,-1,-2,3)
####################
# 時系列の長さ
N <- length(x)
####
i <- 0:(N-1)
i.c <- seq(0,(N-1),length.out=100) 
####################
x.w <- diff(range(x))
x.max <- max(x)+x.w/8
x.min <- min(x)-x.w/8
############
# 高速フーリエ変換
X <- fft(x)
if(N %% 2 == 1){ # 時系列の長さが奇数のとき
 par(mfrow=c((N+1)/2+2,1))
 par(mar=c(2,5,4,4))
 plot(i,x,"o",pch=16,col=2,main="元の時系列x",ylim=c(x.min,x.max))
 for(k in 1:((N+1)/2)){
  f <- (k-1)/N
  if(k == 1){
    x.cos <- rep(Re(X[1])/N,N)
    x.c <- rep(Re(X[1])/N,length(i.c))
    x.sum <- x.c 
  }else{
    x.cos <- 2*abs(X[k])/N*cos(2*pi*f*i+Arg(X[k]))
    x.c <- 2*abs(X[k])/N*cos(2*pi*f*i.c+Arg(X[k]))
    x.sum <- x.sum+x.c
  }
  plot(i,x.cos,pch=1,col=4,ylim=c(x.min,x.max),main=paste("分解した成分 ( f =",(k-1)/N,")"),ylab=paste("成分",k))
  lines(i.c,x.c,col=4,lty=2)
 }
}else{ # 時系列の長さが偶数のとき
 par(mfrow=c(N/2+3,1))
 par(mar=c(2,5,4,4))
 plot(i,x,"o",pch=16,col=2,main="元の時系列x",ylim=c(x.min,x.max))
 for(k in 1:(N/2+1)){
  f <- (k-1)/N
  if(k == 1){
    x.cos <- rep(Re(X[1])/N,N)
    x.c <- rep(Re(X[1])/N,length(i.c))
    x.sum <- x.c 
  }else if(k==N/2+1){
    x.cos <- Re(X[k])/N*cos(pi*i)
    x.c <- Re(X[k])/N*cos(pi*i.c)
    x.sum <- x.sum+x.c  
  }else{
    x.cos <- 2*abs(X[k])/N*cos(2*pi*f*i+Arg(X[k]))
    x.c <- 2*abs(X[k])/N*cos(2*pi*f*i.c+Arg(X[k]))
    x.sum <- x.sum+x.c
  }
  plot(i,x.cos,pch=1,col=4,ylim=c(x.min,x.max),main=paste("分解した成分 ( f =",(k-1)/N,")"),ylab=paste("成分",k))
  lines(i.c,x.c,col=4,lty=2)
 }
}
# 元信号の描画
plot(i,x,"o",pch=16,col=2,ylim=c(x.min,x.max),main="元の時系列(赤)と分解成分の和(青)の比較")
# 分解した成分の和
lines(i.c,x.sum,lty=2,col=4)

このRスクリプトを実行すると,結果は下の図のようになります.

周波数成分への分解

上のRスクリプトの3行目

x <- c(-1,1,-1,-2,3)

で時系列(数列)を設定しています(この値や長さを自分で変えて,いろいろ数値実験してみてください).上の図の最上段が,この時系列のプロットです.その下の3段が分解した成分(青破線).fは周波数なので,f=0の成分は,振動しない定数(分野によっては直流成分と呼ばれる)です.一番下の図が,元の時系列と上の青は線の成分を前部足したグラフの比較です.青い破線の上に赤い点が重なっているので,正しく分解できていることが確認できます.時系列の値や長さを変えても,必ず分解できることが体験できると思います(バグを見つけた場合は教えてください).

数式の説明はまたいつかします.