【Rで時系列解析】1/fノイズ (ピンクノイズ)時系列の生成：中点変位法

今回は，中点変位法 (midpoint displacement algorithm)で，1/fノイズ (ピンクノイズ)のサンプル時系列を生成します．元々，中点変位法は，非整数ブラウン運動のサンプルパス (サンプル時系列)を作る方法ですが，非整数ブラウン運動をわずかにはみ出た領域を考えれば，1/fノイズ (ピンクノイズ)のサンプル時系列を作れます．1/fノイズ (ピンクノイズ)は，ハースト指数が $H=0$ の非整数ブラウン運動じゃないの？と考える人がいるかもしれませんが，非整数ブラウン運動の定義では $0 < H < 1$ なので違います．とはいえ，非整数ブラウン運動用の中点変位法で $H=0$ とすれば，1/fノイズ (ピンクノイズ)を生成できます．前回，中点変位法を説明した記事は，これです．
chaos-kiyono.hatenablog.com

中点変位法を使った1/fノイズ (ピンクノイズ)の生成

　以下の手順で生成できます．1/fノイズ (ピンクノイズ)の時系列を， $\{x(t)\}$ ( $0 \le t \le 1$ )とします．
１．２点 $(0,0)$ ， $(1,0)$ を用意します．つまり， $x(0)=0$ ， $x(1)=0$ とします．

２．２点 $(0,0)$ ， $(1,0)$ の中点に，平均0，標準偏差 $\sigma^2$ の正規乱数を足します．つまり， $x(1/2)$ に，平均0，標準偏差 $\sigma^2$ 正規乱数を代入します． $\sigma^2$ は，1に固定していいです．

３．隣り合う２点の中点に新たな点をとり，その点 $x$ に平均0，標準偏差 $\sigma^2$ の正規乱数を足します．この作業を繰り返すと，1/fノイズ (ピンクノイズ)のサンプル時系列が作れます．

文章の説明だと言葉たらずでポイントが伝わらないかもしれないので，サンプル時系列を作るRスクリプトを，下に掲載しておきます．めちゃくちゃ簡単です．

　もう一つ，パワースペクトルを推定して，ちゃんと，1/fになっているかを確認するスクリプトも載せておきました．結果はこれです．赤い破線が，1/fの傾きですので．ちゃんと，1/fノイズ (ピンクノイズ)になっていることが確認できます．

　この動画も参考にしてください．

www.youtube.com

Rスクリプト

　時系列の生成用のスクリプトです．

################
# 分割する回数
# 15~18くらいを超えると計算時間がかかります
n.step <- 14
################
# 最初の２点
x <- c(0,1)
y <- c(0,0) 
################
i <- 2
for(j in 1:n.step){
 h <- (1/2)^j
 x.tmp <- x[1:i]
 y.tmp <- y[1:i]
 x[(1:(i-1))*2] <- (x.tmp[1:(i-1)]+x.tmp[2:i])/2
 y[(1:(i-1))*2] <- (y.tmp[1:(i-1)]+y.tmp[2:i])/2+rnorm(i-1)
 x[(1:i)*2-1] <- x.tmp
 y[(1:i)*2-1] <- y.tmp
 i <- i + (i-1)
}
# グラフの描画
plot(x,y,"l",col=4,las=1,lwd=2,main=paste("1/f noise (pink noise)"))

　パワースペクトルを推定して検証するスクリプトです．

###########################
# 平均をとるサンプル時系列数
N.samps <- 100
################
# 分割する回数
# 15~18くらいを超えると計算時間がかかります
n.step <- 14
################
# 数値実験
for(i in 1:N.samps){
　#########################
　# 時系列の生成
　################
　# 最初の２点
　x <- c(0,1)
　y <- c(0,0) 
　################
　k <- 2
　for(j in 1:n.step){
　 h <- (1/2)^j
　 x.tmp <- x[1:k]
　 y.tmp <- y[1:k]
　 x[(1:(k-1))*2] <- (x.tmp[1:(k-1)]+x.tmp[2:k])/2
　 y[(1:(k-1))*2] <- (y.tmp[1:(k-1)]+y.tmp[2:k])/2+rnorm(k-1)
　 x[(1:k)*2-1] <- x.tmp
　 y[(1:k)*2-1] <- y.tmp
　 k <- k + (k-1)
　}
 x.sim <- y
 N <- length(x.sim)
　#########################
　# PSDの推定
　tmp <- fft(x.sim)
　PSD <- abs(tmp)^2/N
 if(i == 1){
  PSD.samps <- data.frame(PSD)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
 }
}
# サンプルの平均
PSD <- rowMeans(PSD.samps)
f <- c((0:N)/N)
par(mfrow=c(1,2),cex.lab=1.4)
plot(x,x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main=paste("Synthetic Time Series: m =",m),lwd=2,las=1)
#
plot(f[f > 0 & f <= 1/2],PSD[f > 0 & f <= 1/2],"l",log="xy",col=4,pch=16,xlim=c(f[2],0.5),xaxs="i",xlab="f",ylab="",main="Power Spectrum",lwd=4,las=1)
curve(0.5/x,add=TRUE,lty=2,col=2,lwd=2)
axis(1,at=10^(-15:15)%x%(1:9),label=FALSE,tck=-0.015)
axis(2,at=10^(-15:15)%x%(1:9),label=FALSE,tck=-0.015)