【講義資料7/23】フーリエ変換を体験 - ケィオスの時系列解析メモランダム

　理屈は抜きにして，時系列の「周波数成分への分解」を体験するために，以下の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スクリプトの３行目

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

時系列（数列）を設定しています（この値や長さを自分で変えて，いろいろ数値実験してみてください）．

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

心拍変動のパワースペクトル

　次に，フーリエ変換を使って，心拍変動のパワースペクトルを推定してみます．

# 各自，RRIデータが含まれるフォルダのパスを指定すること
DIR <- r"(D:\Document\大阪大学\講義解説\学問への扉_マチカネゼミ\心拍変動解析やってみる\HUT)"
FN <- c("RRI_HUT0deg.csv","RRI_HUT20deg.csv","RRI_HUT40deg.csv","RRI_HUT60deg.csv")
###################
if (!requireNamespace("signal", quietly = TRUE)) {
  install.packages("signal")
}
# パッケージを読み込む
library(signal)
###################
# 作業フォルダの指定
setwd(DIR)
# ミリ秒まで表示
options(digits.secs=3)
# plotのオプション指定
par(mfrow=c(2,4),cex=1,mar=c(5,5,2,2),pty="m")
###################
# ベッドの角度が0°のデータ
id <- 1
# データの読み込み
TMP <- read.csv(FN[id],header=TRUE)
RRI.0deg <- TMP$RRI
time.0deg <- cumsum(TMP$RRI/1000)-TMP$RRI[1]/1000
# plot
par(pty="m")
plot(time.0deg,RRI.0deg,"l",col=4,ylim=c(600,1100),lwd=2,las=1,xlab="Time [s]",ylab="RRI [ms]",xaxs="i",xlim=c(0,300),main="0 deg")
# パワースペクトルの推定
t.max <- max(time.0deg)
f.samp <- 2
time.resamp <- seq(0,t.max,1/f.samp)
RRI.resamp <- approx(time.0deg,RRI.0deg,xout=time.resamp)$y/1000
psd <- spectrum(RRI.resamp,plot=FALSE)
psd$freq <- psd$freq*f.samp
psd$spec <- sgolayfilt(psd$spec*var(RRI.resamp)/(sum(psd$spec)*(psd$freq[2]-psd$freq[1])),p=2,n=17)
par(pty="s")
plot(psd$freq,psd$spec,"l",xlim=c(0,0.5),xlab="f [Hz]",ylab="S(f)",ylim=c(0,0.035),xaxs="i",yaxs="i",col=4,lwd=2,main="0 deg")
#######################################
# ベッドの角度が20°のデータ
id <- 2
# データの読み込み
TMP <- read.csv(FN[id],header=TRUE)
RRI.20deg <- TMP$RRI
time.20deg <- cumsum(TMP$RRI/1000)-TMP$RRI[1]/1000
# plot
par(pty="m")
plot(time.20deg,RRI.20deg,"l",col=4,ylim=c(600,1100),lwd=2,las=1,xlab="Time [s]",ylab="RRI [ms]",xaxs="i",xlim=c(0,300),main="20 deg")
# パワースペクトルの推定
t.max <- max(time.20deg)
time.resamp <- seq(0,t.max,1/f.samp)
RRI.resamp <- approx(time.20deg,RRI.20deg,xout=time.resamp)$y/1000
psd <- spectrum(RRI.resamp,plot=FALSE)
psd$freq <- psd$freq*f.samp
psd$spec <- sgolayfilt(psd$spec*var(RRI.resamp)/(sum(psd$spec)*(psd$freq[2]-psd$freq[1])),p=2,n=17)
par(pty="s")
plot(psd$freq,psd$spec,"l",xlim=c(0,0.5),xlab="f [Hz]",ylab="S(f)",ylim=c(0,0.035),xaxs="i",yaxs="i",col=4,lwd=2,main="20 deg")
#######################################
# ベッドの角度が40°のデータ
id <- 3
# データの読み込み
TMP <- read.csv(FN[id],header=TRUE)
RRI.40deg <- TMP$RRI
time.40deg <- cumsum(TMP$RRI/1000)-TMP$RRI[1]/1000
# plot
par(pty="m")
plot(time.40deg,RRI.40deg,"l",col=4,ylim=c(600,1100),lwd=2,las=1,xlab="Time [s]",ylab="RRI [ms]",xaxs="i",xlim=c(0,300),main="40 deg")
# パワースペクトルの推定
t.max <- max(time.40deg)
time.resamp <- seq(0,t.max,1/f.samp)
RRI.resamp <- approx(time.40deg,RRI.40deg,xout=time.resamp)$y/1000
psd <- spectrum(RRI.resamp,plot=FALSE)
psd$freq <- psd$freq*f.samp
psd$spec <- sgolayfilt(psd$spec*var(RRI.resamp)/(sum(psd$spec)*(psd$freq[2]-psd$freq[1])),p=2,n=17)
par(pty="s")
plot(psd$freq,psd$spec,"l",xlim=c(0,0.5),xlab="f [Hz]",ylab="S(f)",ylim=c(0,0.035),xaxs="i",yaxs="i",col=4,lwd=2,main="40 deg")
#######################################
# ベッドの角度が60°のデータ
id <- 4
# データの読み込み
TMP <- read.csv(FN[id],header=TRUE)
RRI.60deg <- TMP$RRI
time.60deg <- cumsum(TMP$RRI/1000)-TMP$RRI[1]/1000
# plot
par(pty="m")
plot(time.60deg,RRI.60deg,"l",col=4,ylim=c(600,1100),lwd=2,las=1,xlab="Time [s]",ylab="RRI [ms]",xaxs="i",xlim=c(0,300),main="60 deg")
# パワースペクトルの推定
t.max <- max(time.60deg)
time.resamp <- seq(0,t.max,1/f.samp)
RRI.resamp <- approx(time.60deg,RRI.60deg,xout=time.resamp)$y/1000
psd <- spectrum(RRI.resamp,plot=FALSE)
psd$freq <- psd$freq*f.samp
psd$spec <- sgolayfilt(psd$spec*var(RRI.resamp)/(sum(psd$spec)*(psd$freq[2]-psd$freq[1])),p=2,n=17)
par(pty="s")
plot(psd$freq,psd$spec,"l",xlim=c(0,0.5),xlab="f [Hz]",ylab="S(f)",ylim=c(0,0.035),xaxs="i",yaxs="i",col=4,lwd=2,main="60 deg")
#######################################