時系列解析では，ホワイトノイズとか，ピンクノイズ，レッドノイズとか登場します．私は，ホワイトノイズ以外の呼び方は，ほとんど使いません（色じゃないし）．これらの色の名前がついたノイズは，それらの時系列のパワースペクトルの形状と色のスペクトルの類似性に由来しています．今回は，理屈はおいておいて，Rを使って，ホワイトノイズとか，ピンクノイズを音として体験してみます．
ここでは，tuneRパッケージを使います．Rでパッケージをインストールしていない場合は，
install.packages("tuneR")を実行してください．
データを保存するフォルダは，事前に各自作成して下さい．
例えば，
D:\Document\時系列解析
の場合は，バックスラッシュ("\"，日本語環境では￥マーク)を，"/"に書き直す必要がありますので，下記のRスクリプト3行目にある
DIR <- "...... ここに各自，データを保存したいフォルダへのパスを指定 ......"
の部分を
DIR <- "D:/Document/時系列解析"
に変更してください．
以上の点（パッケージのインストールとDIRの指定）に注意して，このスクリプトを実行してみてください．
############################
# ファイルの保存場所の指定
DIR <- "...... ここに各自，データを保存したいフォルダへのパスを指定 ......"
############################
### 最初に"tuneR"をインストール
#install.packages("tuneR")
require(tuneR)
############################
# 長さ (sec)
L <- 10
# サンプリングレート
f.samp <- 44100
############################
# ノイズ生成
white.nz <- noise(kind="white",stereo=FALSE,duration = L,xunit="time",samp.rate = f.samp)
pink.nz <- noise(kind="pink",stereo=FALSE,duration = L,xunit="time",samp.rate = f.samp)
#############################
# ファイルの書き込み
setwd(DIR)
writeWave(white.nz,file.path("white_noise.wav"))
writeWave(pink.nz,file.path("pink_noise.wav"))
############################
# グラフ 【注意】最初の1000点のみプロットする
par(mfrow=c(2,2)) # 2x2のアレンジでプロット
par(cex.lab=1.4,cex.axis=1.2) # ラベルなどのサイズ指定
par(pty="m")# プロットが長方形
plot((0:1000)/f.samp,white.nz@left[1:1001], type='l',col="gray", xlab='Time [s]',ylab="x",xaxs="i",main="White noise")
plot((0:1000)/f.samp,pink.nz@left[1:1001], type='l',col="pink", xlab='Time [s]',ylab="x",xaxs="i",main="Pink noise")
var.white <- var(white.nz@left)
var.pink <- var(pink.nz@left)
####################################
# パワースペクトルの推定
###
# 解析する時系列 (white noiseの左チャンネル)
x <- white.nz@left
# 時系列の長さ
N <- length(x)
# スペクトル推定のために時系列を分割する部分区間の数
# (精度は高いのは１０以上，最低でも４推奨)
n.psd <- 20
####
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
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)
}
}
# fごとにPSDの平均を計算
psd.mean <- rowMeans(psd)
# psdの片側面積を分散と等しく設定
df <- f[2]-f[1]
psd.a <- sum(psd.mean)*df
psd.est <- psd.mean/psd.a*var.white
# log-logプロットに直線をフィット
fit <- lm(log10(psd.est)~log10(f))
# グラフを描画
plot(log10(f),log10(psd.est),xlab="log10 f",ylab="log10 S(f)","l",las=1,col="gray",main=paste("Estimated PSD (slope=",format(fit$coefficients[[2]],digits=2,nsmall=2),")",sep=""))
# あてはめた直線を破線でプロット
abline(fit,lwd=2,lty=2,col=2)
###
# 解析する時系列 (pink noiseの左チャンネル)
x <- pink.nz@left
# 時系列の長さ
N <- length(x)
####
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
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)
}
}
# fごとにPSDの平均を計算
psd.mean <- rowMeans(psd)
# psdの片側面積を分散と等しく設定
df <- f[2]-f[1]
psd.a <- sum(psd.mean)*df
psd.est <- psd.mean/psd.a*var.white
# log-logプロットに直線をフィット
fit <- lm(log10(psd.est)~log10(f))
# グラフを描画
plot(log10(f),log10(psd.est),xlab="log10 f",ylab="log10 S(f)","l",las=1,col="pink",main=paste("Estimated PSD (slope=",format(fit$coefficients[[2]],digits=2,nsmall=2),")",sep=""))
# あてはめた直線を破線でプロット
abline(fit,lwd=2,lty=2,col=2)
白色ノイズについては，以下のYoutube動画を参考にしてください．
www.youtube.com
