時系列データ解析 上級編 其の四 - Power Spectrum 再考 (2008年11月03日)

時系列データ解析 上級編 其の一 - Hurst Exponent (2008年02月11日) の回で、R-R Interval のパワースペクトラム(power spectrum)の傾き β とトレース(trace)のパワースペクトラム(power spectrum)の傾き α が理論的な値に一致しないという問題にぶち当たった。 ここでは、その原因を解明してみたい。

まず、R-R Interval もあるデータのトレース(trace)であると考えて R-R Interval の微分を取ってみる。 微分を取ると言っても難しいことではなく、単純に隣同士のデータの差を取っただけ。 下の図-43が微分成分をプロットしたもの。

図-43

これだけでは何が何だかわからないので、統計解析(statistical analysis) - 其の壱 (2007年11月05日) で行ったのと同じように 自己相関係数(auto-correlation coefficient) を取ってみる(下の図-44)。 R-R Interval1/fβ のような周波数成分からのみ成っているのであれば、際立ったピークは現れないはずなのだけど、ぱっと観ただけでも たくさんの周期成分が見てとれる。

図-44

微分成分(図-43)のパワースペクトラム(power spectrum)をプロットしてみると、下の図-45のようになる。 オリジナルの R-R Interval データと同じく 周期約3.5秒と15秒辺りに大きなピークがあるのがはっきり分かる。

図-45

そこで、パワースペクトラム(power spectrum)が下記(51)のような特性をしていると仮定して 各パラメータ a, β, ai, bi, ci を最小二乗法でフィットしてみる。 ただし、変数 t は周期(周波数 f の逆数)を表す変数、最初の項が 1/fβ の周波数成分、第二項が複数の周期成分を表す。

(51)

図-44に現れている周期成分を参考にしつつ (51)を図-45の主要成分に対してフィットしてみたのが下の図-46だ。 緑の直線が(51)第一項の 1/fβ 周波数成分、その他の山形の曲線が第二項の複数の周期成分に当たる。

図-46

そして、全ての成分を加えてプロットしたのが下の図-47の緑の曲線になる。 まあまあ、いい感じにフィットできている。 周期約7秒辺りにもう一つ周期成分がありそうだけど、まぁそこは愛嬌って事で。

図-47

さてと、知りたいのはオリジナルの R-R Interval データの 1/fβ 周波数成分だ。 R-R Interval のパワースペクトラム(power spectrum)を改めてプロットしてみると 下の図-48のようになる。

図-48

この図からだとなかなか周期成分が見てとりにくいのだけど、我慢強く丁寧にフィットしていく。 手順は単純で、

  1. 1/fβ 周波数成分のパラメータ a, β を 図-46,47の結果からおおまかに推測する
  2. ピーク値の大きい順に周期成分 ai, bi, ci をフィットしていく
  3. 主要な周期成分をフィットできたら 改めて a, β, ai をフィットする
という感じになる。 a, β の初期値は、それぞれ 2.0e-2 及び 0.7344 とした。

下の図-49,50がその結果になる。 図-49の緑の直線が 1/fβ 周波数成分を表し、欲しかったその傾きは β=0.829742(+/-0.1567) となって、それらしい値が得られた。

β=0.83(44)に適用すると α=2.83 となって図-32右の結果 α=2.7 となんとなく一致する。(笑) めでたしめでたし・・・なのか(?)(汗)

図-49

図-50