【Rで時系列解析】fftを使ってパワースペクトルを推定（ピリオドグラムを計算）

Rの高速フーリエ変換のコマンドfftを使った結果の見方と，その結果を使ったパワースペクトルの計算について説明します．今回は，以前の
chaos-kiyono.hatenablog.com
のつづき（補足）です．パワースペクトルの推定は奥が深い（いろいろとテクニックがある）のですが，ここでは，そんな奥深さは無視して，一番単純なアプローチである「ピリオドグラム」と呼ばれるパワースペクトルの推定法を考えます．ピリオドグラムについては，インターネットで検索して勉強してください．

ここでは実践から学ぶことにします．Rで短い時系列を乱数で作って，fftを計算し，パワースペクトル密度を推定してみます．以下のRスクリプトを実行してみてください．

# 時系列の長さN
N <- 5
# 時系列 x[1], x[2], ..., x[N]に正規乱数を代入
x <- rnorm(N)
####################################
# 離散フーリエ変換 (fftの結果を表示するためにprintを使った)
print(X <- fft(x))
# 周波数 (f.sampはサンプリング周波数，以下では1に設定)
f.samp <- 1
fk <- (0:(N-1))/N * f.samp
# ナイキスト周波数のレンジに移動
fk[fk > 0.5*f.samp] <- fk[fk > 0.5*f.samp] - 1*f.samp
####################################
# ピリオドグラムの計算 （両側パワースペクトル密度の推定）
PSD <- abs(X)^2/N/f.samp
####################################
# 結果のまとめ
data.frame(k=0:(N-1),fk,Element=paste("X[",1:N,"]",sep=""),X=X[1:N],PSD)

上のスクリプトでは，3行目の

N <- 5

で時系列の長さNを設定しています．まずは，奇数の長さにしてみてください．
実行すると，以下のような結果が出力されます（数字XとPSDの数値は各自異なります）．

  k   fk Element                   X       PSD
1 0  0.0    X[1] -2.475980+0.000000i 1.2260949
2 1  0.2    X[2]  1.229301+1.175788i 0.5787317
3 2  0.4    X[3]  0.297968-2.509305i 1.2770792
4 3 -0.4    X[4]  0.297968+2.509305i 1.2770792
5 4 -0.2    X[5]  1.229301-1.175788i 0.5787317

ここで，Xが，xのfft結果で，PSDがパワースペクトル密度の値です．

上の結果で注目してほしい部分は，周波数fkです．マイナスの周波数があります．そうです．パワースペクトルは，正と負の周波数上で定義されたy軸対称なグラフです．軸対称なグラフを描くと，論文や教科書のスペースがもったいないので，右側半分だけ描かれることが多いです．ただし，単に周波数がマイナスの領域を無視するのではなく，その領域のパワースペクトル密度を，プラス側に足してください．つまり，上の例の場合，