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

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

【Rで時系列解析】ヒルベルト変換で振動のエンベロープを抽出

振動する時系列の中には,振動の振幅の変化に役立つ情報がある場合があります.例えば,呼吸信号(換気量や気流)は振動します.そのとき,心拍数に影響を与えるのは,振動の周期ではなく,振幅です.ということで,今回はその「振幅」の情報を取り出すために,ヒルベルト (Hilbert)変換を使ってみます.下の図上の時系列に対して,そのエンベロープ (赤実線)を計算するということです.

ヒルベルト変換で計算したエンベロープ.(上) 元信号 (青実線).(下) 元信号 (青実線)と解析信号の絶対値 (赤実線).

ヒルベルト (Hilbert)変換

 詳しい説明は他の方のホームページなどを見てください.要は,ある周波数で振動している時系列をフーリエ変換して,振動の振幅と位相を計算してやり,元の信号の位相を90°変えた虚部を加えてやるだけです.下の図からイメージしてください.

ヒルベルト変換の考え方

Rで計算

 使ったRスクリプトは下にまとめておきました.

 時系列 \{x(t)\}ヒルベルト変換の手順は以下です.
1. 時系列xフーリエ変換X(f)を計算.Rでは,

X <- fft(x)

2.  i \, {\rm sign}(f)\, X(f)の逆フーリエ変換を計算し,その実部をx_{\rm H}(t)とする.ここで,

 \begin{eqnarray}
{\rm sign}(f) = \left\{\begin{array}{ll} 
 1 & f > 0 \\
 0 & f = 0 \\
 -1 & f < 0 \\
\end{array} \right.
\end{eqnarray}

 Rを使えば,以下のようになります.

# Xの各周波数の位相を90°ずらす
X.H <- 1i*sign(f)*X
# フーリエ逆変換して振動成分の虚部をえる
x.H <- Re(fft(X.H,inverse=TRUE))/N

3. ヒルベルト変換した解析信号z(t)

 z(t) = x(t) + i \, x_{\rm H}(t)

 になります.

 振幅は,|z(t)|です.

気流信号のエンベロープの計算例.

Rスクリプト

ヒルベルト変換の部分

# 事前に時系列を入れておくこと
# x <- c(...)
###
# xの長さ
N <- length(x)
###############################
# Hilbert変換
# 周波数を定義
f <- (0:(N-1))/N
f[f > 1/2] <- f[f > 1/2]-1
###
# xをフーリエ変換
X <- fft(x)
# Xの各周波数の位相を90°ずらす
X.H <- 1i*sign(f)*X
# フーリエ逆変換して振動成分の虚部をえる
x.h <- Re(fft(X.H,inverse=TRUE))/N
###
# ヒルベルト変換の結果
# 虚部の値x.hをもとの信号xに加える
z <- x+1i*x.h
################################
# zの絶対値を計算すると振幅がえられる
amp <- abs(z)


サンプル時系列を使って分析した例.

##############################
# サンプル時系列の生成
# 時系列の長さ
N <- 1000
# サンプリング間隔1で時間を定義
time <- 0:(N-1)
# 振動の総回数
n.osci <- 20
x <- (1-cos(2*pi*time/(N-1)))*sin(2*pi*time/(N/n.osci))
###############################
# Hilbert変換
# 周波数を定義
f <- (0:(N-1))/N
f[f > 1/2] <- f[f > 1/2]-1
###
# xをフーリエ変換
X <- fft(x)
# Xの各周波数の位相を90°ずらす
X.H <- 1i*sign(f)*X
# フーリエ逆変換して振動成分の虚部をえる
x.h <- Re(fft(X.H,inverse=TRUE))/N
###
# ヒルベルト変換の結果
# 虚部の値x.hをもとの信号xに加える
z <- x+1i*x.h
# 振幅
amp <- abs(z)
################################
# プロット
par(mfrow=c(1,1),pty="m")
plot(1:N,x,"l",col=4,lwd=2,las=1,xaxs="i",xlab="time",ylab="x")
# 上のエンベロープ
lines(time,amp,"l",col=2,lwd=2)
# 下のエンベロープ
lines(time,-amp,"l",col=3,lwd=2)