【Rを使った時系列解析】有限の周波数領域で1/fゆらぎ（ピンクノイズ）になる時系列を生成

今回は，周波数 $f$ が $f_1 < f < f_2$ の領域で1/fゆらぎ（ピンクノイズ）の特性を持つサンプル時系列の生成をやります．

詳しい説明はしませんが，1/fゆらぎの生成メカニズムとして，自己共分散関数が指数関数的な減衰を示すサブシステムの和を考えるものが提案されています（下の図参照）．連続時間近似で説明すると，自己共分散関数が
　 $C(\tau) = \sigma^2 e^{-\tau/\tau0}$ 　（ $\sigma^2$ は分散， $\tau_0$ は緩和時間）
のサブシステムを考え，緩和時間 $\tau_0$ が $2 \pi/f_1 > \tau_0 > 2 \pi/f_2$ の領域で， $1/\tau_0$ に比例するように分布していれば，この領域でパワースペクトルは1/fに比例します．

この場合，サブシステムの和のパワースペクトルは，以下の特徴を持ちます（上の図参照）．

$f < f_1$ の領域で，定数
$f_1 < f < f_2$ の領域で， $1/f$ に比例
$f_2 < f < 1/2$ の領域で， $1/f^2$ に比例

この特徴を一つの関数で表すことを考えると，例えば以下のものが考えられます（私が適当に考えたので，既に論文で発表されているかどうかは分かりません）．
　 $S(f) = \displaystyle \frac{1}{(f_1^q+f^q)^{1/q}} \frac{f_2}{(f_2^q+f^q)^{1/q}}$
$q \ge 2$ は，上の3つの領域の接続を制御するパラメタです（値が大きいほど早く変化します）．ここでは，定数の係数は省略してありますので，時系列を生成した後に分散を調整する必要があります．

ここで導入した $S(f)$ のグラフを確かめるＲスクリプトが以下です．

# 区分的に分けた表現
PSD.piecewise <- function(f,f1,f2){
 Sf <- numeric(length(f))
 Sf[f < f1] <- 1/f1 # 低周波数側でホワイト
 Sf[f >= f1 & f <= f2] <- 1/(f[f >= f1 & f <= f2]) #　中間で1/f
 Sf[f > f2] <- f2/f[f > f2]^2 #　高周波数側で1/f^2
 return(Sf)
}
# 一つの関数での表現
PSD.continuous <- function(f,f1,f2,q=2){
 Sf <- numeric(length(f))
 Sf <- 1/((f1^q+f^q)^(1/q))*f2/(f2^q+f^q)^(1/q)
 return(Sf)
}
############################################
# 想定する時系列の長さN
N <- 10001
f <- c((0:((N-1)/2))/N)
# パラメタの設定
f1 <- 5/N
f2 <- 1/2/5
q <- 4
###
psd1 <- PSD.piecewise(f,f1,f2)
psd2 <- PSD.continuous(f,f1,f2,q) 
############################################
# プロット
par(mfrow=c(1,2))
# linear scale
plot(f[f>0],psd1[f>0],"l",col=4,lwd=2,xlab="f",ylab="S(f)")
lines(f[f>0],psd2[f>0],"l",col=2,lty=2,lwd=2)
abline(v=c(f1,f2),lty=2,col=8,lwd=2)
# log-log scale
plot(f[f>0],psd1[f>0],"l",log="xy",col=4,lwd=2,xlab="f",ylab="S(f)")
lines(f[f>0],psd2[f>0],"l",col=2,lty=2,lwd=2)
abline(v=c(f1,f2),lty=2,col=8,lwd=2)

実行結果は下の図のようになり，青実線が区分的に領域を分けたパワースペクトル，赤破線がここで導入した $S(f)$ です．qの値を2，6，10などにしてグラフの変化を調べてください．qの値を大きくすると，青と赤の線が折れ曲がり部分で重なるようになることが分かると思います．

ここで導入した $S(f)$ を使ってサンプル時系列を生成し，意図通りにできているかを確かめるＲスクリプトが以下です．

# 時系列の長さ
N <- 1001
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# 1/fゆらぎの近似パワースペクトル
# [f1,f2]の内側で1/f
PSD0.model <- function(f,f1,f2,q=2){
 Sf <- numeric(length(f))
 Sf <- 1/((f1^q+f^q)^(1/q))*f2/(f2^q+f^q)^(1/q)
 return(Sf)
}
#########################
# パラメタの設定
f1 <- 3/N
f2 <- 1/2/3
q <- 4
f <- c((0:M)/N,-(M:1)/(N))
PSD.model <- PSD0.model(f,f1,f2,q)
#PSD.model[1] <- 0
##############################################
# サンプル時系列の総数
N.ens <- 500
PSD <- c()
ACOR <- c() 
MEAN <- c()
for(i in 1:N.ens){
#########################
# 白色ノイズの生成
WN <- rnorm(N)
#　ホワイトノイズのフーリエ変換
fft.WN <- fft(WN)
#########################
fft.sim <-  sqrt(PSD.model)*fft.WN
#########################
x.sim <- Re(fft(fft.sim,inverse=TRUE))/N
###############################################
# 以下は検証用プロット
################
TMP <- acf(x.sim,plot=FALSE,type="covariance",lag.max=round(N/10))
Ac.y <- TMP$acf[,,1]
X <- fft(x.sim)
S.y <- abs(X)^2/N
 if(i == 1){
  PSD <- S.y
  ACOR <- Ac.y
  MEAN <- mean(x.sim)
 }else{
  PSD <- PSD+S.y
  ACOR <- ACOR + Ac.y
  MEAN <- MEAN + mean(x.sim)
 }
}
PSD <- PSD/N.ens
ACOR <- ACOR/N.ens
MEAN <- MEAN/N.ens
################################################
# グラフ用パラメタ設定
par(mfrow=c(2,2),pty="m",cex=1.1)
plot(f[order(f)],PSD[order(f)],log="y",col=4,"l",lwd=2,xlim=c(-0.5,0.5),xaxs="i",xlab="f",ylab="PSD",main="Estimated Power Spectrum")
lines(f[order(f)],PSD.model[order(f)],col=2,lwd=2,lty=1)
legend("topright", legend=c("Model Power Spectrum"), col=c(2), lty=c(1),lwd=c(2))
plot(f[f>0],PSD[f>0],log="xy",col=4,"l",lwd=2,xlab="f",ylab="PSD",main="Estimated Power Spectrum")
lines(f[f>0],PSD.model[f>0],col=2,lwd=2,lty=1)
abline(v=c(f1,f2),lty=2,col=8,lwd=2)
legend("topright", legend=c("Model Power Spectrum"), col=c(2), lty=c(1),lwd=c(2))
#
fit <- lm(log10(PSD[f > f1 & f < f2])~log10(f[f > f1 & f < f2]))
plot(log10(f[f>0]),log10(PSD[f>0]),col=4,pch=16,"l",lwd=2,xlab="f",ylab="PSD",main="Power Spectrum")
abline(fit,lty=2,col=2,lwd=2)
abline(v=c(log10(f1),log10(f2)),lty=2,col=8,lwd=2)
legend("topright", legend=c(paste("Fitted line with slope =",format(fit$coefficients[[2]],digit=2))), col=c(2), lty=c(2),lwd=c(2))
#################################################
k.lag <- 1:length(ACOR)-1
plot(k.lag[k.lag>0],ACOR[k.lag>0],pch=16,col=4,"o",log="x",lwd=2,xlim=c(1,N/10),xaxs="i",xlab="k",ylab="Autocovariance",main="Estimated Autocovariance")

実行結果の例が以下です．図中の2本の灰色破線が $f_1$ と $f_2$ で指定した領域で，この領域で1/fゆらぎになっていることが分かります．右下は自己共分散関数の推定結果です．