1/fゆらぎの自己相関関数は対数減衰：自己相関はpower lawじゃない

1/fゆらぎ (1/f雑音、フリッカーノイズとも呼ばれます)は、自然界や複雑なシステムに広く観測されるゆらぎの特性です。この雑音は、電気回路、経済データ、生体信号、交通流、さらには音楽や芸術にも現れることが知られており、非常に多様な現象と関係しています。1/fゆらぎはの特徴は、パワースペクトルが周波数 $\displaystyle{ f }$ に対して反比例する、

$\displaystyle{ S(f) \sim \frac{1}{f} }$

という形をとる点にあり、この $1/f$ に比例する構造に由来して「1/f」ゆらぎと呼ばれます。

　通常、時系列データの性質を理解するためには、周波数領域でのパワースペクトルと時間領域での自己相関関数の関係を理解しておくことが重要です。パワースペクトルと自己相関関数は、ウィーナー＝ヒンチンの定理によって互いにフーリエ変換の関係にあり、片方の情報から他方を求めることができます。

　しかし、1/fスペクトルが、低周波数領域全体に続くと、周波数 $\displaystyle{ f \to 0 }$ で発散するため、全体としてのエネルギー（あるいは分散）が無限大になってしまい、数学的には非定常な過程となってしまいます (パーセバルの定理によりパワースペクトルの面積は時系列の分散に一致します)。これでは理論的に扱いにくく、現実的ではないため、通常は有限の周波数範囲（バンド）において1/f特性を持つ過程を考えます。

　今回は、この有限バンドでの1/fスペクトルに対応する自己相関関数 $\displaystyle{ C(\tau) }$ を導出します。導出の工夫として、1/fスペクトルが異なる時定数をもつローレンツ型スペクトルの重ね合わせとして得られるという仮定を用い、対応する時間領域の自己相関関数の形を求めます。その結果、自己相関関数が対数的に時間遅れ $\displaystyle{ \tau }$ に依存するという性質が導かれます。論文などでたまに、1/fゆらぎの自己相関関数はべき的に減衰すると書かれていることがありますが、間違いです。べき的に減衰しません。

　以下ではまず、1/fスペクトルの定義とその周波数領域での性質について確認し、その後、ローレンツスペクトルを重ね合わせる方法を用いて、自己相関関数の導出を進めていきます。

　参考にさせていただいた論文はこれです。

On the correlation function of 1/f noise - ScienceDirect

もくじ

有限バンドでの 1/f スペクトルの導入
ローレンツスペクトルの重ね合わせによる 1/f スペクトルの構成
- スペクトルの重ね合わせ
- 中間領域での arctan の近似
対応する自己相関関数の導出
- 積分の評価のために自己相関関数を微分する
- スケーリング領域での近似
Rで数値実験

有限バンドでの 1/f スペクトルの導入

　低周波数側の全領域で、厳密に

$\displaystyle{ S(f) = \frac{A}{f} }$

というスペクトルを仮定すると、積分

$\displaystyle{ \sigma^2 = \int_0^\infty S(f)\,df = A \int_0^\infty \frac{1}{f}\, df }$

が発散し、分散が無限大になってしまうため、定常過程ではなくなります。
そこで、実際には「低周波側 $\displaystyle{ f _ l }$ 以下では白色雑音(定数)、高周波側 $\displaystyle{ f _ h }$ 以上では $\displaystyle{ f^ {-2} }$ に漸近する」ようにカットオフを入れて

$\displaystyle{ S(f) = \begin{cases} \text{constant} & (0 \le f < f_l), \\ \dfrac{A}{f} & (f_l \le f \le f_h), \\ \propto f^{-2} & (f > f_h), \end{cases} }$

とし、有限のスケーリング領域と有限の分散をもつ、有限バンド内での $\displaystyle{ 1/f }$ スペクトルを考えます。以下では

$\displaystyle{ f_l = \frac{1}{\tau_2}, \quad f_h = \frac{1}{\tau_1} }$

と定義し、この範囲で $\displaystyle{ \frac{A}{f} }$ になるように設定しています。

ローレンツスペクトルの重ね合わせによる 1/f スペクトルの構成

1次自己回帰過程のように自己相関関数が、指数関数的に減衰する確率過程を考えます。

$\displaystyle{ C_L(\tau) = \sigma^2 e^{-\tau/\tau_0} }$

ここで、 $\displaystyle{ \tau _ 0 }$ は時定数で、相関がほぼ0になるまでの、時間を調整するパラメタです。

　この過程のパワースペクトルは、ローレンツ型スペクトルとして知られる、

$\displaystyle{ S_L(f) = \frac{4\tau_0 \sigma^2}{1 + (2\pi f)^2 \tau_0^2} }$

で与えられます。これは、ウィーナーヒンチンの定理で計算すれば求められます。

　ここでは、 $\displaystyle{ \tau _ 0 }$ が異なるたくさんのローレンツスペクトルを足し合わせることで、1/fスペクトルを表現します。

　具体的には， $\displaystyle{ \tau _ 0 }$ を $\displaystyle{ \tau _ 1 }$ （最大周波数側）から $\displaystyle{ \tau _ 2 }$ （最小周波数側）まで変化させるとき，分布関数 $\displaystyle{ g(\tau _ 0) }$ を

$\displaystyle{ g(\tau_0) = \frac{1}{\ln(\tau_2 / \tau_1)} \cdot \frac{1}{\tau_0} }$

の形で与えると、想定した周波数領域で $\displaystyle{ \propto 1/f }$ となるスペクトルを構成できます。

スペクトルの重ね合わせ

有限バンドでのスペクトル $\displaystyle{ S(f) }$ は

$\displaystyle{ \begin{aligned} S(f) &= \int_{\tau_1}^{\tau_2} S_L(f)\, g(\tau_0)\,d\tau_0 \\ &= \int_{\tau_1}^{\tau_2} \frac{4\tau_0 \sigma^2}{1 + (2\pi f)^2\tau_0^2} \cdot \frac{1}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{\tau_0}\,d\tau_0 \\ &= \frac{4\sigma^2}{\ln(\tau_2/\tau_1)} \int_{\tau_1}^{\tau_2} \frac{1}{1 + (2\pi f)^2\tau_0^2} \,d\tau_0 \end{aligned} }$

この積分は、公式

$\displaystyle{ \int \frac{1}{1+(ax)^2}\,dx = \frac{1}{a}\,\arctan(ax) }$

を用いると

$\displaystyle{ \int_{\tau_1}^{\tau_2} \frac{d\tau_0}{1 + (2\pi f)^2\tau_0^2} = \left[ \frac{1}{2\pi f}\arctan(2\pi f \,\tau_0) \right]_{\tau_0 = \tau_1}^{\tau_2}. }$

となることがわかりますので、

$\displaystyle{ S(f) = \frac{4 \sigma^2}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{2\pi f} \bigl[\arctan(2\pi f\,\tau_0)\bigr]_{\tau_1}^{\tau_2} }$

となります。

中間領域での arctan の近似

ここで、

$\displaystyle{ f_l = \frac{1}{\tau_2} \approx 0 \ll f \ll \frac{1}{\tau_1} = f_h \approx \frac{1}{2} }$

のように、可能な限り広いスケーリング領域をとるとすれば、

$\displaystyle{ \arctan(2\pi f\,\tau_1) \approx 0 }$

$\displaystyle{ \arctan(2\pi f\,\tau_2) \approx \frac{\pi}{2} }$

と近似できます。したがって、

$\displaystyle{ \bigl[\arctan(2\pi f\,\tau_0)\bigr]_{\tau_1}^{\tau_2} \approx \frac{\pi}{2} - 0 = \frac{\pi}{2} }$

です。

　その結果、1/fスペクトル

$\displaystyle{ S(f) \approx \frac{4\sigma^2}{\ln(\tau_2/\tau_1)}\cdot \frac{1}{2\pi f} \cdot \frac{\pi}{2} = \frac{\sigma^2}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{f} }$

がえられます。

対応する自己相関関数の導出

　パワースペクトルと同様に、自己相関関数 $\displaystyle{ C(\tau) }$ もローレンツ型自己相関 $\displaystyle{ C _ L(\tau) }$ を同じ重み付け $\displaystyle{ g(\tau _ 0) }$ で重ね合わせることで求められます。

　計算してみると、

$\displaystyle{ \begin{aligned} C(\tau) &= \int_{\tau_1}^{\tau_2} C_L(\tau)\,g(\tau_0)\, d\tau_0 \\ &= \int_{\tau_1}^{\tau_2} \sigma^2 e^{-\tau/\tau_0} \cdot \frac{1}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{\tau_0} \, d\tau_0 \\ &= \frac{\sigma^2}{\ln(\tau_2/\tau_1)} \int_{\tau_1}^{\tau_2} \frac{1}{\tau_0} e^{-\tau/\tau_0} \, d\tau_0 \end{aligned} }$

となります。最後の積分はぱっと見、めんどくさそうです。そこで、直接評価するのをやめて、 $\displaystyle{ \tau }$ で微分するテクニックを使います。

積分の評価のために自己相関関数を微分する

　上で登場した積分を直接評価するのをやめて、一度、 $\displaystyle{ \tau }$ で微分してから、 $\displaystyle{ \tau _ 0 }$ で積分してみます。どの文字で微分して、どの文字で積分するのかに注意して計算してください。以下では、積分と微分の順序を入れ替えるテクニックを使っています。

$\displaystyle{ \begin{aligned} \frac{dC(\tau)}{d\tau} &= \frac{d}{d\tau} \Biggl[ \frac{\sigma^2}{\ln(\tau_2/\tau_1)} \int_{\tau_1}^{\tau_2} \frac{1}{\tau_0} e^{-\tau/\tau_0}\, d\tau_0 \Biggr]\\ &= \frac{\sigma^2}{\ln(\tau_2/\tau_1)} \int_{\tau_1}^{\tau_2} \frac{1}{\tau_0}\, \frac{d}{d\tau}\bigl(e^{-\tau/\tau_0}\bigr) \, d\tau_0\\ &= -\frac{\sigma^2}{\ln(\tau_2/\tau_1)} \int_{\tau_1}^{\tau_2} \frac{1}{\tau_0^2} e^{-\tau/\tau_0} \, d\tau_0\\ &= -\frac{\sigma^2}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{\tau} \Bigl[ e^{-\tau/\tau_0} \Bigr]_{\tau_0=\tau_1}^{\tau_2}\\ &= -\frac{\sigma^2}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{\tau} \Bigl( e^{-\tau/\tau_2} - e^{-\tau/\tau_1} \Bigr)\\ \end{aligned} }$

のように、微分したものは解析的に積分できます。

スケーリング領域での近似

興味があるのは， $\displaystyle{ \tau _ 1 }$ より十分大きく，かつ、 $\displaystyle{ \tau _ 2 }$ より十分小さい $\displaystyle{ \tau }$ の範囲での振る舞いです．大雑把に思い切った近似をして，

$\displaystyle{ \tau \gg \tau_1 \;\Rightarrow\; e^{-\tau/\tau_1} \approx e^{-\infty} \approx 0 }$

$\displaystyle{ \tau \ll \tau_2 \;\Rightarrow\; e^{-\tau/\tau_2} \approx e^{- 0} \approx 1 }$

とします。

　こうすると、スケーリング領域の主要な部分では、 $\tau$ の依存性は無視出来て、

$\displaystyle{ e^{-\tau/\tau_2} - e^{-\tau/\tau_1} \approx 1 - 0 = 1 }$

となり、自己相関関数の微分が、

$\displaystyle{ \frac{dC(\tau)}{d\tau} \simeq -\frac{\sigma^2}{\ln(\tau_2/\tau_1)} \cdot \frac{1}{\tau} }$

となります。

　この結果を使って、両辺を $\tau$ で積分すると

$\displaystyle{ \begin{aligned} C(\tau) &\approx -\frac{\sigma^2}{\ln(\tau_2/\tau_1)} \int \frac{d\tau}{\tau} + C_0 \\ &= -\frac{\sigma^2}{\ln(\tau_2/\tau_1)} \ln|\tau| + C_0 \\ &\sim -\ln|\tau| \end{aligned} }$

となります。ここで、 $C _ 0$ は定数です。

　十分に広い周波数範囲で1/f スペクトルになる場合には、対数的に減衰する自己相関関数が導かれます。

Rで数値実験

　数値実験でここで導いたことを確認するRスクリプトをおまけでのせておきます。結果は、下の図のようになると思います。

　上の図では、左と中央がパワースペクトルの両対数プロットです。中央の周波数領域でばっちり1/fになってます。右の図が自己相関関数です。横軸のみ対数をとってありますので、直線に見える領域が対数的減衰になっています。

###############################
# シミュレーション設定と関数定義
###############################

# --- 時系列データの長さ設定 ---
N <- 10001            # 基本となる時系列長
if (N %% 2 == 0) {    # 奇数でない場合は奇数に調整
  N <- N + 1
}
M <- (N - 1) / 2      # 正の周波数成分の数（0～0.5程度の分解能）

# --- 1/fゆらぎ近似パワースペクトル関数 ---
# 入力:
#   f  : 周波数ベクトル（正負両方）
#   f1 : 下限周波数（この付近より低いと1/f特性が始まる）
#   f2 : 上限周波数（この付近より高いと減衰が始まる）
#   q  : 形状パラメータ（大きいほど急峻な切れ方）
PSD_model_func <- function(f, f1, f2, q = 2) {
  # fが負の場合も同じ値となるように絶対値を使用
  1 / ((f1^q + abs(f)^q)^(1/q)) * f2 / ((f2^q + abs(f)^q)^(1/q))
}

# --- パラメータの設定 ---
f1 <- 5 / N         # 下限周波数
f2 <- 1 / 4         # 上限周波数
q  <- 4             # パワースペクトルの形状パラメータ

# 周波数ベクトルの作成（正の周波数と対称な負の周波数）
f <- c((0:M) / N, -(M:1) / N)
# 理論上のパワースペクトルを計算
PSD_theo <- PSD_model_func(f, f1, f2, q)


###############################
# シミュレーションによる時系列生成
###############################

N_ens <- 100        # アンサンブル数（生成するサンプル数）

# --- 積算用変数の事前初期化 ---
PSD_sum  <- numeric(N)              # 各サンプルの周波数領域PSDを積算
lag_max <- round(N / 4)              # 自己相関を計算する最大ラグ
ACOR_sum <- numeric(lag_max + 1)      # 自己相関関数（lag 0～lag_max）の積算

# --- アンサンブルループ ---
for(i in 1:N_ens){
  # ① 白色ノイズ（正規乱数）による時系列生成
  WN <- rnorm(N)
  
  # ② 白色ノイズのフーリエ変換
  fft_WN <- fft(WN)
  
  # ③ 理論PSD形状を持たせるため、各周波数成分に sqrt(PSD_theo) を乗算
  fft_sim <- sqrt(PSD_theo) * fft_WN
  
  # ④ 逆フーリエ変換により時系列を再構成（実数部分のみ抽出、Rのifftは1/Nで割る）
  x_sim <- Re(fft(fft_sim, inverse = TRUE)) / N
  
  # ⑤ 自己相関関数の計算（ラグ0～lag_max; plot=FALSEでプロットは作成しない）
  acf_obj <- acf(x_sim, plot = FALSE, type = "covariance", lag.max = lag_max)
  acor <- as.vector(acf_obj$acf[, 1, 1])
  
  # ⑥ x_simからのPSD推定（フーリエ変換後のパワーを算出）
  PSD_est <- abs(fft(x_sim))^2 / N
  
  # ⑦ 各アンサンブルの結果を積算
  PSD_sum  <- PSD_sum + PSD_est
  ACOR_sum <- ACOR_sum + acor
}

# --- アンサンブル平均の算出 ---
PSD_avg  <- PSD_sum / N_ens
ACOR_avg <- ACOR_sum / N_ens


###############################
# プロット（3つのパネル）
###############################

# プロット領域を横に3分割して設定
par(mfrow = c(1, 3),mar=c(5,5,5,2), pty = "m", cex = 1.1, las = 1,cex.lab=1.5,cex.axis=1.3)

## (1) 推定PSDと理論PSDの比較
# 正の周波数成分のみプロット（対数-対数プロット）
plot(f[f > 0], PSD_avg[f > 0], log = "xy", type = "l", col = 4, lwd = 2,
     xlab = "f", ylab = "PSD", main = "Estimated Power Spectrum", xaxt = "n", yaxt = "n")
# 理論パワースペクトルを重ね描き
lines(f[f > 0], PSD_theo[f > 0], col = 2, lwd = 2, lty = 1)
# f1, f2 の位置に垂直線を描画
abline(v = c(f1, f2), lty = 2, col = 8, lwd = 2)
# 軸はシンプルに描画（必要に応じてカスタムラベル可）
axis(1)
axis(2)
legend("topright", legend = "Model Power Spectrum", col = 2, lty = 1, lwd = 2)


## (2) 対数変換後のPSDと線形フィッティング
# 正の周波数成分をlog10スケールに変換してプロット
log_f <- log10(f[f > 0])
log_PSD <- log10(PSD_avg[f > 0])
plot(log_f, log_PSD, type = "l", col = 4, lwd = 2,
     xlab = "log10(f)", ylab = "log10(PSD)", main = "Power Spectrum (log-log)")
# f1～f2の範囲内（かつ正の周波数）のデータを抽出して線形回帰
sel <- f > f1 & f < f2 & f > 0
fit <- lm(log10(PSD_avg[sel]) ~ log10(f[sel]))
abline(fit, lty = 2, col = 2, lwd = 2)
# f1, f2 に対応する位置（log10スケール）に垂直線を追加
abline(v = log10(c(f1, f2)), lty = 2, col = 8, lwd = 2)
legend("topright", legend = paste("Fitted slope =", format(coef(fit)[2], digits = 2)),
       col = 2, lty = 2, lwd = 2)


## (3) 自己相関関数のプロットとフィッティング
# ラグ（0～lag_max）の定義
lags <- 0:lag_max
# ラグ0はlog10(0)が無限大になるため、1以上のラグのみをプロット
plot(log10(lags[-1]), ACOR_avg[-1], type = "o", pch = 16, col = 4, lwd = 2,
     xlim = log10(c(1, N / 3)), xlab = expression(log[10](tau)), ylab = expression(C(tau)),
     main = "Autocorrelation Function")
# 線形フィッティングのため、ラグ5以上の部分を補間して抽出
log10tau <- seq(log10(1/f2), log10(1/f1/5), length.out = 20)
interp_acor <- approx(log10(lags[lags > 5]), ACOR_avg[lags > 5], xout = log10tau)$y
fit_acor <- lm(interp_acor ~ log10tau)
points(log10tau,interp_acor,pch=16,col="lightpink")
abline(fit_acor, col = 2, lwd = 2, lty = 2)
# f1, f2 の逆数に対応するラグ位置（対数スケール）に垂直線を描画
abline(v = log10(c(1/f1, 1/f2)), lty = 2, col = 8, lwd = 2)
# 自己相関がゼロになる位置に水平線を追加
abline(h = 0, lty = 2, col = 8, lwd = 2)