Rで時系列のパワースペクトルを推定したいとき，もっともお気楽なのは，"spectrum"というコマンドを使うことだと思います．今回は，心拍変動 (RR間隔時系列)のスペクトル解析で"spectrum"コマンドを使った実例を紹介します．ちなみに，例として使用したRR間隔時系列は，１０分の長さのもので，4 Hz (0.25 秒間隔)でリサンプリングされています (リサンプリングについては関連記事参照)．ヘッドアップティルト試験で，０°，20°，40°の３条件で計測されたものです．今回は，HFとか，LFの説明は省略します．真のパワースペクトルは，どんな形か想像して推定結果を観察してください．
Rの"spectrum"コマンドでは，ペリオドグラムの計算の前処理として，トレンドを除く，窓関数をかける処理が含まれています．後処理として，全パワーも調整されているようです (私の印象です)．
1. 自己回帰過程の当てはめを使ったパラメトリック推定
まずは，時系列に自己回帰過程を当てはめ，推定された自己回帰過程の係数を使ってパワースペクトルを推定した結果です．
サンプリング間隔dtの時系列をxとすれば，
tmp <- spectrum(x, plot="false",method="ar") freq <- tmp$freq/dt psd <- tmp$spec*dt
2. ペリオドグラムを使ったノンパラメトリック推定
ペリオドグラムの平滑化でいろいろと設定を変えた結果を以下になれべておきます．ここでは，
tmp <- spectrum(x, plot="false",spans=c( ... )) freq <- tmp$freq/dt psd <- tmp$spec*dt
として，"spans"オプションの設定を変えています．
