教科書的には，時系列解析の基本はパワースペクトルです．とはいってもこれは，ガウス過程（正規分布に従う過程）中心主義，そして，自己相関関数中心主義の立場では，パワースペクトルが世界のリーダーだよね，ということです（中心主義とかは私が勝手に権威批判的に言っているだけです）．パワースペクトルのグラフは，時系列を周波数の異なる成分（sinやcosの波）に分解して，横軸に周波数，縦軸をその周波数成分の振幅の2乗に比例する量を描いたものです．時系列解析では，パワースペクトルを波への分解だけのイメージでとらえるだけでは，真の解釈にはたどり着けないというのが私の持論です．しかし，そんなのどーでもいいという人も多いと思いますので，ここでは，波の成分とパワースペクトルの関係，音の高さと周波数の関係を感じられるRスクリプトを載せておきます．
ここでは，tuneRパッケージを使います．Rでパッケージをインストールしていない場合は，
install.packages("tuneR")を実行してください．
今回のRスクリプトを実行すると音を鳴らすためのファイルが生成されますので，そのデータを保存するフォルダを，事前に各自作成して下さい．
例えば，
D:\Document\時系列解析
の場合は，バックスラッシュ("\"，日本語環境では￥マーク)を，"/"に書き直す必要がありますので，下記のRスクリプト3行目にある
DIR <- "...... ここに各自，データを保存したいフォルダへのパスを指定 ......"
の部分を
DIR <- "D:/Document/時系列解析"
に変更してください．
以上に注意して，下記のスクリプトを実行すると，"A2_tone.wav"，"A4_tone.wav"，"A7_tone.wav"という3つのファイルが，DIRで指定したフォルダに作られます．これらのファイルを音楽再生アプリ（メディアプレーヤーなど）で再生すると純音を聴くことができます．「ピー」と音がなるので，音量を小さめにしておいてください．
############################
# ファイルが保存されるフォルダの指定
DIR <- "...... ここに各自，データを保存したいフォルダへのパスを指定 ......"
############################
### 最初に"tuneR"をインストール
#install.packages("tuneR")
require(tuneR)
############################
# 長さ (sec)
L <- 10
# サンプリングレート
f.samp <- 44100
############################
# 純音生成
sin.A2 <- sine(freq=110,stereo=TRUE,duration = L,xunit="time",samp.rate = f.samp)
sin.A4 <- sine(freq=440,stereo=TRUE,duration = L,xunit="time",samp.rate = f.samp)
sin.A7 <- sine(freq=3520,stereo=TRUE,duration = L,xunit="time",samp.rate = f.samp)
#############################
# ファイルの書き込み
setwd(DIR)
writeWave(sin.A2,file.path("A2_tone.wav"))
writeWave(sin.A4,file.path("A4_tone.wav"))
writeWave(sin.A7,file.path("A7_tone.wav"))
############################
# グラフ 【注意】最初の1000点のみプロットする
par(mfrow=c(2,3)) # 2x2のアレンジでプロット
par(cex.lab=1.4,cex.axis=1.2) # ラベルなどのサイズ指定
par(pty="m")# プロットが長方形
rbw <- rainbow(6)
plot((0:1000)/f.samp,sin.A2@left[1:1001], type='l',col=rbw[1], xlab='Time [s]',ylab="x",xaxs="i",main="A2 tone (110 Hz)",las=1)
plot((0:1000)/f.samp,sin.A4@left[1:1001], type='l',col=rbw[3], xlab='Time [s]',ylab="x",xaxs="i",main="A4 tone (440 Hz)",las=1)
plot((0:1000)/f.samp,sin.A7@left[1:1001], type='l',col=rbw[6], xlab='Time [s]',ylab="x",xaxs="i",main="A7 tone (3520 Hz)",las=1)
####################################
# パワースペクトルの推定
###
# 解析する時系列
x <- sin.A2@left
# 時系列の長さ
N <- length(x)
# スペクトル推定のために時系列を分割する部分区間の数
# (周波数の解像度を上げるため今回は分割なし)
n.psd <- 1
####
n.sub <- floor(N/n.psd)
for(i in 1:n.psd){
if(i == 1){
f <- spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$freq*f.samp
psd <- data.frame(spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
}else{
psd <- cbind(psd,spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
}
}
psd.mean <- rowMeans(psd)/f.samp
# グラフを描画
plot(f[-1],psd.mean[-1],log="x",xlab="f",ylab="S(f)","l",las=1,col=rbw[1],main="Estimated PSD",xlim=c(50,22000))
###
# 解析する時系列
x <- sin.A4@left
# 時系列の長さ
N <- length(x)
# スペクトル推定のために時系列を分割する部分区間の数
# (周波数の解像度を上げるため今回は分割なし)
n.psd <- 1
####
n.sub <- floor(N/n.psd)
for(i in 1:n.psd){
if(i == 1){
f <- spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$freq*f.samp
psd <- data.frame(spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
}else{
psd <- cbind(psd,spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
}
}
psd.mean <- rowMeans(psd)/f.samp
# log-logプロットに直線をフィット
fit <- lm(log10(psd.mean)~log10(f))
# グラフを描画
plot(f[-1],psd.mean[-1],log="x",xlab="f",ylab="S(f)","l",las=1,col=rbw[3],main="Estimated PSD",xlim=c(50,22000))
###
# 解析する時系列
x <- sin.A7@left
# 時系列の長さ
N <- length(x)
# スペクトル推定のために時系列を分割する部分区間の数
# (周波数の解像度を上げるため今回は分割なし)
n.psd <- 1
####
n.sub <- floor(N/n.psd)
for(i in 1:n.psd){
if(i == 1){
f <- spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$freq*f.samp
psd <- data.frame(spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
}else{
psd <- cbind(psd,spectrum(x[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
}
}
psd.mean <- rowMeans(psd)/f.samp
# log-logプロットに直線をフィット
fit <- lm(log10(psd.mean)~log10(f))
# グラフを描画
plot(f[-1],psd.mean[-1],log="x",xlab="f",ylab="S(f)","l",las=1,col=rbw[6],main="Estimated PSD",xlim=c(50,22000))
