振動する時系列の中には，振動の振幅の変化に役立つ情報がある場合があります．例えば，呼吸信号（換気量や気流）は振動します．そのとき，心拍数に影響を与えるのは，振動の周期ではなく，振幅です．ということで，今回はその「振幅」の情報を取り出すために，ヒルベルト (Hilbert)変換を使ってみます．下の図上の時系列に対して，そのエンベロープ (赤実線)を計算するということです．
ヒルベルト (Hilbert)変換
詳しい説明は他の方のホームページなどを見てください．要は，ある周波数で振動している時系列をフーリエ変換して，振動の振幅と位相を計算してやり，元の信号の位相を90°変えた虚部を加えてやるだけです．下の図からイメージしてください．
Rで計算
使ったRスクリプトは下にまとめておきました．
時系列 のヒルベルト変換の手順は以下です．
1. 時系列xのフーリエ変換を計算．Rでは，
X <- fft(x)
2. の逆フーリエ変換を計算し，その実部を
とする．ここで，
Rを使えば，以下のようになります．
# Xの各周波数の位相を90°ずらす X.H <- 1i*sign(f)*X # フーリエ逆変換して振動成分の虚部をえる x.H <- Re(fft(X.H,inverse=TRUE))/N
3. ヒルベルト変換した解析信号は
になります．
振幅は，です．
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)