【Rで時系列解析】1/f 型時系列の生成：指数βを自由に設定

$1/f^{\beta}$ 型ゆらぎは、生体システムに広く見られる特徴的なフラクタルゆらぎです。心拍ゆらぎ、脳波パワー、呼吸リズム、歩行リズムなど、多くの生理信号では $\beta$ が 1 に近い 1/f ゆらぎが観測されます。 $1/f^{\beta}$ 型ゆらぎを評価する手法としては、パワースペクトル解析、Detrended Fluctuation Analysis (DFA)、Detrending Moving Average Analysis (DMA)、Higuchi のフラクタル次元 (HFD) など、複数の方法が提案されています。しかし、実際に時系列データを解析する際には、「どの手法を使うべきか？」と迷うことも少なくありません。

　本記事では、こうした解析法の特性を体験的に理解するために、FFT（高速フーリエ変換）を用いて、

$\displaystyle{ S(f) \propto \frac{1}{f^\beta} }$

に従う離散時系列を生成する方法を解説します。 $1/f^{\beta}$ 型ゆらぎのモデルとしては fractional Brownian motion、fractional Gaussian noise、Autoregressive fractionally integrated moving average (ARFIMA) などが知られていますが、これらのモデルにはパラメタの取りうる範囲に制約があります。それに対して、今回紹介する生成法では $\beta$ の値に制限がなく、幅広いスケーリング指数をもつゆらぎを自在に作れる点が大きな利点です。

　以下では、1/f 型時系列を FFT で生成する基本的な考え方、実際に利用できる R スクリプト、DMA と Higuchi のフラクタル解析を用いたデモンストレーション、について順に説明していきます。

1. 基本アイデア：ホワイトノイズのスペクトルを「1/f形に整形する」
2. Higuchi フラクタル次元と DMA スケーリング指数の性能評価
3.【Rスクリプト①】時系列の生成
4.【Rスクリプト②】DMAとHFDのスケーリング指数の推定可能領域

1. 基本アイデア：ホワイトノイズのスペクトルを「1/f形に整形する」

　Rスクリプトは解説記事の最後に掲載しました。ここでは、『3.【Rスクリプト①】時系列の生成』で示した時系列の生成方法について説明します。
　この R スクリプトでやっていることを、コードの流れに沿って整理すると：

長さ $N$ のホワイトノイズ時系列をつくる
そのフーリエ成分に、所望の PSD の平方根を掛けて「スペクトルの形」を $1/f^ \beta$ 型に変える
逆 FFT で時系列に戻す
得られた時系列の PSD を推定し、モデル PSD と重ねて確認する

というシンプルな構造になっています。

1.1 離散時系列と正規化周波数軸

　まず、長さ $N$ の離散時系列を扱います。

N <- 1001              # 元の時系列長
N <- round(N / 2) * 2 + 1
M <- (N - 1) / 2

N を必ず奇数 $N = 2M+1$ にそろえることで、周波数軸を「0, 正の周波数 M 個, 負の周波数 M 個」という 対称な形で扱えるようにしています。
　周波数軸は、離散フーリエ変換に対応する 正規化周波数（サンプリング周波数で割った単位）で

f.pos <- (0:M) / N      # 0, 1/N, ..., M/N
f.neg <- -(M:1) / N     # -M/N, ..., -1/N
f     <- c(f.pos, f.neg)

となっています。これが

$\displaystyle{ f = \left( 0, \frac{1}{N}, \dots, \frac{M}{N}, -\frac{M}{N}, \dots, -\frac{1}{N} \right) }$

という「両側スペクトル用」の周波数軸に対応します。

1.2 ターゲット PSD：正則化された $1/f^ \beta$ 型

　理想的な $\displaystyle{ 1/f^ \beta }$ スペクトル

$\displaystyle{ S(f) \propto \frac{1}{|f|^\beta} }$

は、 $\displaystyle{ f \to 0 }$ で発散します。有限長の数値シミュレーションではこれが扱いにくいので、このスクリプトでは次のように 低周波側をなだらかに丸めた形を使っています：

$\displaystyle{ S_{\text{model}}(f) \;\propto\; \frac{1}{\bigl(f_0^2 + f^2\bigr)^{\beta/2}}. }$

これを実装しているのが

PSD1.model <- function(f, f0, beta = 1) {
  1 / (f0^2 + f^2)^(beta / 2)
}

です。

f0 <- 1 / (2 * N) は、「最も低い離散周波数ステップ（1/N）の約半分」程度の小さなカットオフ周波数で、 $f = 0$ 近傍の発散をやわらげる役割を持っています。
さらに

  PSD.model[1L] <- 0

として DC 成分（ $f=0$ ）のパワーは完全に 0 にしています。これにより、「巨大な DC 成分」や「極端なトレンド」が乗るのを避けています。

1.3 PSD の面積を 1 にそろえる

　同じ $\beta$ でも、f0 の値や周波数刻みによって、PSD の「総エネルギー」が変わってしまいます。スクリプトでは、比較のしやすさのために、PSD の面積を 1 に正規化しています：

if (normalize_psd) {
  df   <- 1 / N
  area <- sum(PSD.model) * df
  if (area > 0) {
    PSD.model <- PSD.model / area
  }
}

これは、離散近似として

$\displaystyle{ \sum_k S_{\text{model}}(f_k)\,\Delta f \;\approx\; 1\,\quad (\Delta f = 1/N) }$

を満たすように、PSD.model 全体を比例縮尺している、ということです。

こうしておくと、(\beta) を変えても 時系列の分散や振幅スケールが極端に変わりにくくなり、
DFA / DMA / HFD などの解析では、「形（スケーリング指数）の違い」に注目しやすくなります。

1.4 ホワイトノイズのスペクトル：平らな PSD

　次に、長さ $N$ の 白色ガウス雑音を生成します：

WN <- rnorm(N)
fft.WN <- fft(WN)

WN は平均 0・分散一定の独立なガウス乱数列。
その FFT である fft.WN は、周波数ごとの複素フーリエ成分 (W[k]) に対応します。

ホワイトノイズの PSD の特徴は：

$\displaystyle{ S_{\text{white}}(f) \approx \text{const.} }$

つまり、周波数によらずほぼ一定で、「フラットなスペクトル」を持っています。この「フラット」なスペクトルを、これから $\displaystyle{ 1/f^ \beta }$ 形に「整形」していきます。

1.5 スペクトル整形： $\displaystyle{\sqrt{S _ {\text{model}}(f _ k)}}$ を掛ける

　ここがこの方法の核心です。
　ホワイトノイズのフーリエ成分を $\displaystyle{ W[k] }$ 、ターゲット PSD を $\displaystyle{ S _ {\text{model}}(f _ k) }$ とすると、

$\displaystyle{ X_{\text{sim}}(k) = \sqrt{S_{\text{model}}(f_k)} \cdot W[k] }$

と定義します。すると、

$\displaystyle{ |X_{\text{sim}}(k)|^2 = S_{\text{model}}(f_k)\, |W[k]|^2. }$

ホワイトノイズでは、 $\displaystyle{ |W[k]|^ 2 }$ の期待値は「周波数に依存しない定数」なので、

$\displaystyle{ \mathbb{E}\bigl[|X_{\text{sim}}(k)|^2\bigr] \;\propto\; S_{\text{model}}(f_k) }$

となり、平均としてはターゲット PSD に従うスペクトルが得られます。これを R でそのまま書いたのが：

fft.sim <- sqrt(PSD.model) * fft.WN

です。

$\displaystyle{\sqrt{S _ {\text{model}}(f _ k)}}$ を掛けることで「振幅の大きさ」を $\displaystyle{1/f^ \beta}$ 型に整形し、
位相はホワイトノイズがもともと持っているランダムな位相をそのまま使うことで、現実的なランダム時系列になっています。

さらに PSD.model は 実数かつ偶関数なので、元の fft.WN が持っていた「共役対称性」を壊さず、結果として実数時系列が得られます。　　

1.6 逆 FFT で時系列へ戻す

スペクトル整形したフーリエ成分 fft.sim を、逆 FFT で時系列に戻します：

x.sim <- Re(fft(fft.sim, inverse = TRUE)) / N

fft(..., inverse = TRUE) は、数学的な逆 FFT の $N$ 倍になって返ってくるので、/ N で正規化しています。
小さな数値誤差によって虚部が出るので、Re() で実部だけを取り出しています。

この x.sim が、 $1/f^ \beta$ 型 PSD を持つ「合成ゆらぎ時系列」です。

1.7 合成時系列の PSD を推定し、平滑化する

　生成された時系列 x.sim から、実際にどんな PSD が得られているかを確認します。

X  <- fft(x.sim)
Sy <- Mod(X)^2 / N

これは $\displaystyle{|X[k]|^ 2/N}$ による単純なスペクトル推定です。
1 本の時系列からの推定なので、かなりギザギザします。

そこで、正の周波数だけ抜き出し、

idx.pos       <- which(f > 0)
f.pos.only    <- f[idx.pos]
Sy.pos        <- Sy[idx.pos]
Sy.pos.smooth <- running_mean(Sy.pos, smooth_window)

とし、running_mean() で ローカルな移動平均平滑化をかけて、 $1/f^ \beta$ 型の傾きが見やすくなるようにしています。

running_mean() は自前実装で、累積和を使って $O(N)$ で計算する素朴な移動平均です。

1.8 両対数プロット

　最後に、モデル PSD と推定 PSD を 両対数プロットで比較します。

2. Higuchi フラクタル次元と DMA スケーリング指数の性能評価

　ここからは、生成した時系列データを使って、Higuchi フラクタル次元 (HFD)と DMA スケーリング指数が、真のべき指数 $\beta$ に対してどのように振る舞うかを系統的に調べる R スクリプトの概要を説明します。実際のスクリプトは、最後に『4.【Rスクリプト②】DMAとHFDのスケーリング指数の推定可能領域』として掲載してあります。

Higuchi フラクタル次元 (HFD)と DMA スケーリング指数の性能評価の結果。HFDは $1 \lt \beta \lt 3$ 、DMAは $-1 \lt \beta \lt 3$ の範囲であれば適用可能。

Higuchi フラクタル次元 (HFD)と DMA スケーリング指数の性能評価の結果。HFDは $1 \lt \beta \lt 3$ 、DMAは $-1 \lt \beta \lt 3$ の範囲であれば適用可能。

2.1 $1/f ^\beta$ 型時系列生成関数の再利用

　まず、前節と同じ形のモデル PSD

$\displaystyle{ S(f) \propto \frac{1}{(f_0^2 + f^2)^{\beta/2}} }$

を与える関数 PSD1.model() を定義し、それを用いた時系列生成関数 simulate_1overf_ts() を用意します。

simulate_1overf_ts(
  N             = 10001,
  beta          = 1,
  f0            = NULL,
  normalize_psd = TRUE
)

N：系列長（内部で必ず奇数 (2M+1) に補正）
beta：目標とする $1/f ^\beta$ の指数
f0：低周波の正則化用カットオフ（指定しなければ $1/(2N)$ ）
normalize_psd：PSD の面積を 1 にそろえるかどうか

内部では、

周波数軸 (f) を作る
PSD1.model() でターゲット PSD を計算し、DC を 0 に
面積が 1 になるように正規化（任意）
ホワイトノイズの FFT 結果を $1/f ^\beta$ 型に変形
逆 FFT でサンプル時系列 x.sim を得る

という処理を行います。

2.2 Higuchi フラクタル次元 $D _ H$ の推定

　higuchi_fd() 関数は、与えられた時系列 x に対して Higuchi の方法でフラクタル次元 $D _ H$ を推定します。

higuchi_fd <- function(x, k.max = 100) { ... }

時系列を、間引き間隔 $k = 1, 2, \dots$ ごとにサブ系列に分解
各 (k) について「キザキザの長さ」 $L(k)$ を求める
$L(k) \propto k^{-D}$ を仮定し、

$\displaystyle{ \log L(k) = \text{const} + D \log(1/k) }$

の直線当てはめの傾きから $D$ を推定
　スクリプトでは、 $k$ を 1〜k.max すべてではなく、 exp(seq(0, log(k.max), length.out=20)) のように 対数間隔で間引いて評価し、計算量を抑えています。

2.3 DMA スケーリング指数 $\alpha$ の推定

　dma_alpha() 関数は、Detrending Moving Average (DMA) によるスケーリング解析を行います。

dma_alpha <- function(
  x,
  scale.min = 9,
  scale.max = NULL,
  n.scales  = 20
) { ... }

まず時系列 x に対し、窓長 $s$ ごとの中心移動平均（stats::filter）を計算
元の信号から移動平均を引いた残差の RMS を

$\displaystyle{ F(s) = \sqrt{\langle r^2 \rangle} }$

として求める

$\displaystyle{F(s) \propto s^ \alpha}$ を仮定し、

$\displaystyle{ \log F(s) = \text{const} + \alpha \log s }$

の傾きからスケーリング指数 $\alpha$ を推定

ここでは scale.min〜scale.max の範囲で対数間隔にスケールをとり、奇数長の窓（.../2)*2+1）になるように調整しています。

なお、メインループでは dma_alpha(cumsum(x)) として、プロファイル（累積和）に DMA をかける形になっていることに注意してください。

2.4 $\beta$ を振って平均を取る　

　メイン部分では、

betas <- seq(-2, 5, by = 0.2)  # 真の beta のリスト
n.ens <- 20                    # 各 beta での試行回数（ensemble サイズ）

とし、 $\beta = -2, -1.8, \dots, 5$ の範囲で系統的に調べます。

二重ループの中身は：

simulate_1overf_ts() で $1/f^ \beta$ 型時系列 x を生成
平均を引いて中心化：x <- x - mean(x)
higuchi_fd(x) で Higuchi 次元 $D _ H$ を推定
dma_alpha(cumsum(x)) で DMA スケーリング指数 $\alpha$ を推定
それぞれの結果を行列 D.higuchi.mat, alpha.dma.mat に格納

最後に、各 $\beta$ について平均をとります：

D.higuchi.mean <- rowMeans(D.higuchi.mat, na.rm = TRUE)
alpha.dma.mean <- rowMeans(alpha.dma.mat, na.rm = TRUE)

2.5 $\beta$ vs 推定スケーリング指数の関係をプロット

　最後に、横軸に真の $\beta$ 、縦軸に

DMA で得た $\alpha$
Higuchi フラクタル次元 $D _ H$

の平均をプロットします。

plot(
  betas, alpha.dma.mean,
  type = "b",
  pch  = 16,
  col  = 4,
  ylim = ylim.all,
  xlab = expression(beta),
  ylab = "Estimated scaling exponent",
  main = "Higuchi FD and DMA exponent vs 1/f^beta (ensemble mean)"
)

lines(betas, D.higuchi.mean, type = "b", pch = 17, col = 2)

この図を見ることで、

真の $\beta$ に対して DMA の $\alpha$ や Higuchi の $D _ H$ が どの範囲で線形に反応するか
どのレンジで バイアスや飽和が出るか

といった、各解析法の「測定レンジ」や「クセ」を、 $1/f^ \beta$ 型の合成データを使って直感的に把握できます。

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。

3.【Rスクリプト①】時系列の生成

############################################################
# 1/f^beta 型ゆらぎ時系列シミュレーション（FFT 法）
#  - モデル PSD: S(f) ∝ 1 / (f0^2 + f^2)^(beta/2)
#  - 白色雑音のフーリエ成分に sqrt(PSD) を掛けて整形
#  - PSD の数値積分（面積）を 1 に正規化（任意）
#  - 推定 PSD を周波数近傍でローカル平滑化（移動平均）
############################################################

## -----------------------------
## 基本パラメータ設定
## -----------------------------

N <- 1001              # 元の時系列長
# N を必ず奇数 (2M+1) にする
N <- round(N / 2) * 2 + 1
M <- (N - 1) / 2        # 片側の点数（周波数の正側の個数）

beta <- 1               # 1/f^beta のべき指数
f0   <- 1 / (2 * N)     # 正則化用カットオフ周波数

# PSD の面積を 1 に正規化するかどうか
normalize_psd <- TRUE

# PSD を平滑化する窓幅（正の周波数側の Sy に対して）
# 例: 5 → 自分を含めて前後 2 点ずつの移動平均
smooth_window <- 5

## -----------------------------
## 1/f ゆらぎの近似パワースペクトル
##   S(f) = 1 / (f0^2 + f^2)^(beta/2)
## -----------------------------

PSD1.model <- function(f, f0, beta = 1) {
  # f   : 周波数（-0.5〜0.5 の正規化周波数を想定）
  # f0  : カットオフ周波数（0 での発散を避けるための定数）
  # beta: 1/f^beta の指数
  1 / (f0^2 + f^2)^(beta / 2)
}

## 周波数軸の生成（両側スペクトル：0, 1/N, ..., M/N, ..., -1/N）
f.pos <- (0:M) / N         # 0 〜 正の周波数
f.neg <- -(M:1) / N        # 負の周波数
f     <- c(f.pos, f.neg)   # 長さ N (= 2M+1)

## モデル PSD（両側）
PSD.model <- PSD1.model(f, f0, beta)

# DC 成分（f = 0）は無相関成分として 0 にしておく
PSD.model[1L] <- 0

## -----------------------------
## PSD の面積を 1 に正規化（数値積分）
##   面積 ≈ sum_k S(f_k) * Δf
##   ここでは Δf = 1/N
## -----------------------------

if (normalize_psd) {
  df   <- 1 / N
  area <- sum(PSD.model) * df
  if (area > 0) {
    PSD.model <- PSD.model / area
  }
}

## -----------------------------
## 白色雑音の生成とフーリエ変換
## -----------------------------

# 白色ガウス雑音
WN <- rnorm(N)

# 白色雑音の FFT（両側スペクトル）
fft.WN <- fft(WN)

## -----------------------------
## 1/f スペクトルを持つフーリエ成分を生成
##   X_sim(k) = sqrt(PSD(k)) * X_WN(k)
##   → PSD に従うスペクトル形状を付与
## -----------------------------

# PSD が実数・偶関数なので、実数時系列の共役対称性は維持される
fft.sim <- sqrt(PSD.model) * fft.WN

## -----------------------------
## 逆 FFT で時系列を生成
##   R の fft(inverse=TRUE) は 1/N を掛けていないので、
##   ここで 1/N を掛けて正規化
## -----------------------------

x.sim <- Re(fft(fft.sim, inverse = TRUE)) / N

## -----------------------------
## 得られた時系列のパワースペクトル（推定）
##   S_y(f_k) = |X_k|^2 / N
## -----------------------------

X  <- fft(x.sim)
Sy <- Mod(X)^2 / N   # |X|^2 / N

## -----------------------------
## ローカル平滑化用の移動平均関数
##   - パッケージを使わず O(N) で計算
##   - 端点では利用可能な範囲だけで平均
## -----------------------------

running_mean <- function(y, window) {
  n <- length(y)
  if (window <= 1L || n == 0L) return(y)

  # 奇数窓にしておくと「中心」が定義しやすい
  if (window %% 2L == 0L) {
    window <- window + 1L
  }
  half <- (window - 1L) / 2L

  cs  <- c(0, cumsum(y))  # cs の長さは n+1
  out <- numeric(n)

  for (i in seq_len(n)) {
    left  <- max(1L, i - half)
    right <- min(n, i + half)
    # 区間 [left, right] の和 = cs[right+1] - cs[left]
    out[i] <- (cs[right + 1L] - cs[left]) / (right - left + 1L)
  }
  out
}

## 正の周波数だけを抜き出す（DC を除外）
idx.pos       <- which(f > 0)
f.pos.only    <- f[idx.pos]
Sy.pos        <- Sy[idx.pos]
Sy.pos.smooth <- running_mean(Sy.pos, smooth_window)

## -----------------------------
## 描画設定
## -----------------------------

par(mfrow = c(1, 2), cex.lab=1.3)

## -----------------------------
## (1) 推定 PSD とモデル PSD（両対数プロット）
## -----------------------------

# 0 より大きい範囲だけを見る（log 軸のため）
x.range <- range(f.pos.only[f.pos.only > 0])
y.range <- range(Sy.pos.smooth[Sy.pos.smooth > 0])

# 10 の何乗か（指数）に変換
x.e.min <- floor(log10(x.range[1]))
x.e.max <- ceiling(log10(x.range[2]))
y.e.min <- floor(log10(y.range[1]))
y.e.max <- ceiling(log10(y.range[2]))

# 主目盛（10^n）と補助目盛（10^n × 1〜9）
x.major <- 10^(x.e.min:x.e.max)
y.major <- 10^(y.e.min:y.e.max)

x.minor <- as.vector(outer(10^(x.e.min:x.e.max), 1:9))
x.minor <- x.minor[x.minor >= x.range[1] & x.minor <= x.range[2]]

y.minor <- as.vector(outer(10^(y.e.min:y.e.max), 1:9))
y.minor <- y.minor[y.minor >= 10^y.e.min & y.minor <= 10^y.e.max]

# 両対数プロット本体（軸目盛りは自分で描くので xaxt, yaxt = "n"）
plot(
  f.pos.only,
  Sy.pos,
  type = "l",
  log  = "xy",
  ylim = 10^c(y.e.min, y.e.max),
  xlab = "frequency (normalized)",
  ylab = "Power spectral density",
  main = sprintf("Simulated 1/f^%.2f time series", beta),
  xaxs = "i",
  yaxs = "i",
  col  = "grey",
  lwd  = 1,
  xaxt = "n",
  yaxt = "n"
)

lines(f.pos.only, Sy.pos.smooth, col = 4, lwd = 2)
lines(f.pos.only, PSD.model[idx.pos], col = 2, lwd = 2, lty = 2)

legend(
  "bottomleft",
  legend = c("Raw estimated PSD", "Smoothed PSD", "Model PSD"),
  col    = c("grey", 4, 2),
  lwd    = c(1, 2, 2),
  lty    = c(1, 1, 2),
  bty    = "n"
)

###########################################
# 対数軸の描画（主目盛＝10^n、補助目盛＝10^n×1〜9）
###########################################

## x 軸：補助目盛（ラベルなし）
axis(1, at = x.minor, labels = FALSE, tck = -0.01)

## x 軸：主目盛（10^n 表記）
exponents.x <- x.e.min:x.e.max               # 指数ベクトル
x.labels    <- parse(text = paste0("10^", exponents.x))

axis(
  1,
  at       = x.major,
  labels   = x.labels,
  las      = 1,
  cex.axis = 1.2,
  tck      = -0.02
)

## y 軸：補助目盛（ラベルなし）
axis(2, at = y.minor, labels = FALSE, tck = -0.01)

## y 軸：主目盛（10^n 表記）
exponents.y <- y.e.min:y.e.max
y.labels    <- parse(text = paste0("10^", exponents.y))

axis(
  2,
  at       = y.major,
  labels   = y.labels,
  las      = 1,
  cex.axis = 1.2,
  tck      = -0.02
)
### 対数目盛（ここまで）

## -----------------------------
## (2) 時系列波形
## -----------------------------
plot(
  seq_len(N),
  x.sim,
  type = "l",
  xlim = c(0,N),
  xlab = "n",
  ylab = "x[n]",
  main = sprintf("Simulated 1/f^%.2f time series", beta),
  xaxs = "i",
  col  = 4,
  lwd  = 2
)

4.【Rスクリプト②】DMAとHFDのスケーリング指数の推定可能領域

############################################################
# 1/f^beta 型ゆらぎ時系列シミュレーション（FFT 法）
# + Higuchi フラクタル次元 & DMA スケーリング解析
#
#  - beta を -2 〜 3 まで 0.2 刻みで変化
#  - 各 beta について n.ens 本の時系列を生成
#  - Higuchi 次元 D_H と DMA exponent alpha を推定
#  - ensemble 平均を beta の関数としてプロット
#
#  パッケージは使用しない (stats は base)
############################################################

## --------------------------------------------------------
## モデル PSD: S(f) ∝ 1 / (f0^2 + f^2)^(beta/2)
## --------------------------------------------------------
PSD1.model <- function(f, f0, beta = 1) {
  # f   : 周波数（-0.5〜0.5 の正規化周波数）
  # f0  : カットオフ周波数（低周波での正則化）
  # beta: 1/f^beta のべき指数
  1 / (f0^2 + f^2)^(beta / 2)
}

## --------------------------------------------------------
## ローカル平滑化用の移動平均（O(N)）
##   - 端点では利用可能な範囲だけで平均
##   - window が偶数なら 1 つ増やして奇数にする
## --------------------------------------------------------
running_mean <- function(y, window) {
  n <- length(y)
  if (window <= 1L || n == 0L) return(y)

  if (window %% 2L == 0L) {
    window <- window + 1L
  }
  half <- (window - 1L) / 2L

  cs <- c(0, cumsum(y))  # 長さ n+1 の累積和
  out <- numeric(n)

  for (i in seq_len(n)) {
    left  <- max(1L, i - half)
    right <- min(n, i + half)
    out[i] <- (cs[right + 1L] - cs[left]) / (right - left + 1L)
  }
  out
}

## --------------------------------------------------------
## 1/f^beta 型時系列を生成する関数（プロットなし）
## --------------------------------------------------------
simulate_1overf_ts <- function(
  N             = 10001,   # 元の時系列長（必ず奇数に補正される）
  beta          = 1,       # 1/f^beta の指数
  f0            = NULL,    # カットオフ周波数（NULL なら 1/(2N)）
  normalize_psd = TRUE     # PSD の面積を 1 に正規化するか
) {

  ## ---- N を必ず奇数 (2M+1) にする
  N <- round(N / 2) * 2 + 1
  M <- (N - 1) / 2  # 片側の点数（周波数の正側の個数）

  ## ---- f0 が指定されていなければデフォルト値
  if (is.null(f0)) {
    f0 <- 1 / (2 * N)
  }

  ## ---- 周波数軸の生成（両側スペクトル）
  f.pos <- (0:M) / N        # 0 〜 正の周波数
  f.neg <- -(M:1) / N       # 負の周波数
  f     <- c(f.pos, f.neg)  # 長さ N (= 2M+1)

  ## ---- モデル PSD（両側）
  PSD.model <- PSD1.model(f, f0, beta = beta)

  # DC 成分（f = 0）は 0 にしておく（無相関成分として除去）
  PSD.model[1L] <- 0

  ## ---- PSD の面積を 1 に正規化（任意）
  if (normalize_psd) {
    df   <- 1 / N
    area <- sum(PSD.model) * df
    if (area > 0) {
      PSD.model <- PSD.model / area
    }
  }

  ## ---- 白色ガウス雑音とその FFT
  WN     <- rnorm(N)
  fft.WN <- fft(WN)

  ## ---- スペクトル整形：sqrt(PSD) を掛ける
  fft.sim <- sqrt(PSD.model) * fft.WN

  ## ---- 逆 FFT で時系列を生成
  x.sim <- Re(fft(fft.sim, inverse = TRUE)) / N

  x.sim
}

############################################################
# Higuchi のフラクタル次元推定
############################################################
higuchi_fd <- function(x, k.max = 100) {
  x <- as.numeric(x)
  N <- length(x)
  if (N < 2) stop("length(x) must be >= 2")

  k.max <- min(k.max, floor((N - 1) / 10))  # 極端に大きすぎる k を避ける
  k.scale <- unique(round(exp(seq(0,log(k.max),length.out=20))))
  Lk <- numeric(length(k.scale))

  for (k in k.scale) {
    Lm <- numeric(k)

    # m = 1,...,k ごとに L_m(k) を計算
    for (m in 1:k) {
      idx <- seq(m, N, by = k)
      n.m <- length(idx)
      if (n.m > 1) {
        diff.sum <- sum(abs(diff(x[idx])))
        # Higuchi の正規化
        Lm[m] <- diff.sum * (N - 1) / ((n.m - 1) * k)
      } else {
        Lm[m] <- NA_real_
      }
    }

    # 有効な m の平均
    Lk[k] <- mean(Lm, na.rm = TRUE)
  }

  # 0 や NA を除外
  valid <- which(Lk > 0 & is.finite(Lk))
  if (length(valid) < 2) return(NA_real_)

  Lk <- Lk[valid]
  kk <- valid

  # L(k) ∝ k^{-D} → log L(k) = const + D * log(1/k)
  y <- log(Lk)
  x.reg <- log(1 / kk)

  D <- coef(lm(y ~ x.reg))[2]+1
  as.numeric(D)
}

############################################################
# DMA スケーリング指数 (alpha) 推定
############################################################
dma_alpha <- function(
  x,
  scale.min = 9,
  scale.max = NULL,
  n.scales  = 20
) {
  x <- as.numeric(x)
  N <- length(x)
  if (is.null(scale.max)) {
    scale.max <- floor(N / 10)   # 長さの 1/10 くらいまで
  }

  # 対数間隔のスケールを生成（重複は除去）
  scales <- unique(round(exp(seq(
    log(scale.min),
    log(scale.max),
    length.out = n.scales
  ))/2))*2+1
  scales <- scales[scales >= 2]
  nS <- length(scales)
  if (nS < 2) stop("Too few valid scales for DMA")

  Fv <- numeric(nS)  # F(s)

  for (i in seq_along(scales)) {
    s <- scales[i]

    # 中心移動平均（sides = 2）
    ma <- stats::filter(x, rep(1 / s, s), sides = 2)

    # 残差
    r <- x - ma

    # NA を除外して RMS を取る
    Fv[i] <- sqrt(mean(r^2, na.rm = TRUE))
  }

  # 対数回帰: log F(s) = const + alpha * log s
  y <- log(Fv)
  x.reg <- log(scales)
  fit <- lm(y ~ x.reg)
  alpha <- coef(fit)[2]

  list(alpha = as.numeric(alpha), scales = scales, F = Fv)
}

############################################################
# メイン：beta を振って ensemble 平均のスケーリング指数を計算
############################################################

# --- 解析条件の設定 --------------------------------------
N     <- 1001    # 時系列長
k.max <- 100     # Higuchi の最大 k
n.ens <- 20      # ensemble 本数

# beta = -2, -1.8, ..., 5  (0.2 刻み)
betas <- seq(-2, 5, by = 0.2)
n.beta <- length(betas)

# 結果を格納する行列
D.higuchi.mat <- matrix(NA_real_, nrow = n.beta, ncol = n.ens)
alpha.dma.mat <- matrix(NA_real_, nrow = n.beta, ncol = n.ens)

# --- ループ本体 ------------------------------------------
for (i in seq_along(betas)) {
  beta <- betas[i]
  cat(sprintf("beta = %5.2f ...\n", beta))

  for (j in 1:n.ens) {

    # 1/f^beta 時系列を生成（simulate_1overf_ts）
    x <- simulate_1overf_ts(
      N             = N,
      beta          = beta,
      normalize_psd = TRUE
    )

    # 念のため平均を引いておく
    x <- x - mean(x)

    # Higuchi フラクタル次元
    D.higuchi.mat[i, j] <- higuchi_fd(x, k.max = k.max)

    # DMA スケーリング指数
    dma.res <- dma_alpha(cumsum(x))
    alpha.dma.mat[i, j] <- dma.res$alpha
  }
}

# --- ensemble 平均 ---------------------------------------
D.higuchi.mean <- rowMeans(D.higuchi.mat, na.rm = TRUE)
alpha.dma.mean <- rowMeans(alpha.dma.mat, na.rm = TRUE)

############################################################
# プロット：横軸 beta，縦軸 スケーリング指数（Higuchi & DMA）
############################################################
par(mfrow=c(1,1),las=1)
ylim.all <- range(c(D.higuchi.mean, alpha.dma.mean), na.rm = TRUE)

plot(
  betas, alpha.dma.mean,
  type = "b",
  pch  = 16,
  col  = 4,
  ylim = ylim.all,
  xlab = expression(beta),
  ylab = "Estimated scaling exponent",
  main = "Higuchi FD and DMA exponent vs 1/f^beta (ensemble mean)"
)

lines(betas, D.higuchi.mean, type = "b", pch = 17, col = 2)

legend(
  "bottomright",
  legend = c("DMA exponent (alpha)", "Higuchi fractal dimension D"),
  col    = c(4, 2),
  pch    = c(16, 17),
  lty    = 1,
  bty    = "n"
)

ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで時系列解析】1/f 型時系列の生成：指数βを自由に設定

1. 基本アイデア：ホワイトノイズのスペクトルを「1/f形に整形する」

1.1 離散時系列と正規化周波数軸

1.2 ターゲット PSD：正則化された $1/f^ \beta$ 型

1.3 PSD の面積を 1 にそろえる

1.4 ホワイトノイズのスペクトル：平らな PSD

1.5 スペクトル整形： $\displaystyle{\sqrt{S _ {\text{model}}(f _ k)}}$ を掛ける

1.6 逆 FFT で時系列へ戻す

1.7 合成時系列の PSD を推定し、平滑化する

1.8 両対数プロット

2. Higuchi フラクタル次元と DMA スケーリング指数の性能評価

2.1 $1/f ^\beta$ 型時系列生成関数の再利用

2.2 Higuchi フラクタル次元 $D _ H$ の推定

2.3 DMA スケーリング指数 $\alpha$ の推定

2.4 $\beta$ を振って平均を取る

2.5 $\beta$ vs 推定スケーリング指数の関係をプロット

3.【Rスクリプト①】時系列の生成

4.【Rスクリプト②】DMAとHFDのスケーリング指数の推定可能領域

1. 基本アイデア：ホワイトノイズのスペクトルを「1/f形に整形する」

1.1 離散時系列と正規化周波数軸

1.2 ターゲット PSD：正則化された 型

1.3 PSD の面積を 1 にそろえる

1.4 ホワイトノイズのスペクトル：平らな PSD

1.5 スペクトル整形： を掛ける

1.6 逆 FFT で時系列へ戻す

1.7 合成時系列の PSD を推定し、平滑化する

1.8 両対数プロット

2. Higuchi フラクタル次元と DMA スケーリング指数の性能評価

2.1 型時系列生成関数の再利用

2.2 Higuchi フラクタル次元 の推定

2.3 DMA スケーリング指数 の推定

2.4 を振って平均を取る

2.5 vs 推定スケーリング指数の関係をプロット

3.【Rスクリプト①】時系列の生成

4.【Rスクリプト②】DMAとHFDのスケーリング指数の推定可能領域

1.2 ターゲット PSD：正則化された $1/f^ \beta$ 型

1.5 スペクトル整形： $\displaystyle{\sqrt{S _ {\text{model}}(f _ k)}}$ を掛ける

2.1 $1/f ^\beta$ 型時系列生成関数の再利用

2.2 Higuchi フラクタル次元 $D _ H$ の推定

2.3 DMA スケーリング指数 $\alpha$ の推定

2.4 $\beta$ を振って平均を取る　

2.5 $\beta$ vs 推定スケーリング指数の関係をプロット