【Rでスペクトル解析３】パワースペクトル推定法のノンパラメトリックとパラメトリック

有限長の離散時系列データからパワースペクトル（Power Spectral Density: PSD）を推定する方法は、大きくノンパラメトリック法とパラメトリック法に分けることができます。
　ノンパラメトリック法では、時系列データに特定の確率モデルを仮定せず、時系列を有限区間でフーリエ変換してその振幅の二乗からパワースペクトルを推定する方法や、時系列から自己相関関数を推定し、それをフーリエ変換することでパワースペクトルを求める方法が用いられます。
　一方、パラメトリック法では、時系列が自己回帰（AR）モデルなどの確率モデルに従うと仮定し、そのモデルパラメタを推定した上で、得られたモデルの理論的なパワースペクトルを計算します。
　本稿では、これら二つのアプローチの違いを、Rを用いて実際にパワースペクトルを計算することを通して体験的に確認します。
　以下のRスクリプトを実行する前に、次の変数が既に与えられているものとします。

f.samp : サンプリング周波数 [Hz]
time : 時刻ベクトル（長さ N）
x : 時系列データ（長さ N）

手持ちのデータがない場合は付録の「3.【付録】心拍変動のLF、HF成分シミュレーション」に示したRスクリプトを実行してください。以下の解析では、DC成分（平均値）の除去を前処理として行っています。

1. ノンパラメトリック法
2. パラメトリック法
3. 付録：心拍変動のLF、HF成分シミュレーション
- シミュレーションのパラメタ設定
- 出力と利用方法

1. ノンパラメトリック法

　ここでは、パワースペクトル（PSD）の 3 つのノンパラメトリック推定法

ペリオドグラム（periodogram）
Welch 法（区間分割＋平均化）
自己相関のフーリエ変換（コレログラム: correlogram）

を、Rスクリプト例とともに示します。いずれも「両側スペクトル（two-sided）」として計算し、非負周波数側（0からNyquist周波数まで）だけを描画します（one-sided 化の 2 倍補正は行いません）。

1.1 ペリオドグラム（periodogram）

　ペリオドグラムは、有限長データをそのままフーリエ変換し、フーリエ振幅の二乗から PSD を推定する最も基本的な方法です。

長所：実装が簡単、周波数分解能が高い
短所：分散が大きく、スペクトルがギザギザになりやすい（単発の推定が不安定）

　以下が、Rスクリプト例です。両側PSDを計算し半分だけ描画する仕様ですので、厳密な片側PSDではないことに注意してください。

# ============================================================
# ペリオドグラムによるパワースペクトル推定
#
# 【入力として仮定する変数】
#   f.samp : サンプリング周波数 [Hz]
#   time   : 時刻ベクトル（長さ N）
#   x      : 時系列データ（長さ N）
#
# 【出力】
#   freq.h : 非負周波数側（0 ～ ナイキスト周波数）の周波数軸
#   Pperiod : パワースペクトル密度
#            （両側スペクトルのスケーリングのまま，半分のみを描画）
# ============================================================

# 時刻ベクトルと時系列データの長さが一致しているかを確認
stopifnot(length(time) == length(x))

# データ長 N とサンプリング間隔 dt
N  <- length(x)
dt <- 1 / f.samp

# ------------------------------------------------------------
# 直流成分（平均値）を除去
#   f = 0 におけるスペクトルの過度な突出を避けるため，
# ------------------------------------------------------------
x0 <- x - mean(x, na.rm = TRUE)

# ------------------------------------------------------------
# 高速フーリエ変換（FFT）
#   R の fft() は正規化されていない点に注意
# ------------------------------------------------------------
X <- fft(x0)

# ------------------------------------------------------------
# 両側ペリオドグラムの計算
#   周波数 f_k = k * f.samp / N における
#   パワースペクトル密度の推定値
#   S(f_k) = |X_k|^2 / (N * f.samp)
# ------------------------------------------------------------
P2 <- (Mod(X)^2) / (N * f.samp)

# ------------------------------------------------------------
# 周波数軸（両側スペクトル用）
#   k = 0, 1, ..., N-1 に対応
# ------------------------------------------------------------
freq <- (0:(N - 1)) * (f.samp / N)

# ------------------------------------------------------------
# 非負周波数側（0 ～ ナイキスト）のみを取り出す
# ※ one-sided スペクトルのような 2 倍補正は行わない
# ------------------------------------------------------------
k_max <- floor(N / 2)
idx   <- 1:(k_max + 1)   # R は 1 始まりのインデックス
freq.h <- freq[idx]
Pperiod <- P2[idx]

# ------------------------------------------------------------
# 結果の描画
# ------------------------------------------------------------
ymax <- max(P2)
plot(freq.h, Pperiod,
     type = "l", col = 2, lwd = 2,
     xlab = "Frequency [Hz]",
     ylab = "PSD",
     main = "Periodogram (two-sided, half plotted)",
     ylim = c(0, ymax * 1.2),
     xaxs = "i", yaxs = "i")

ペリオドグラムの計算例。時系列は「3. 付録」に示したRスクリプトを使用して事前に生成。

1.2 Welch 法（区間分割＋平均化）

　Welch 法は、ペリオドグラムの弱点である激しいギザギザを改善するために、

データを複数区間に分割（通常は重なりあり）
各区間で（窓掛けした）ペリオドグラムを計算
それらを平均して PSD を得る

という手順を使います。ここでは 区間数 n.sub を指定し、50%オーバーラップで分割します。n.sub を増やせば、ギザギザは減りますが、周波数の解像度が下がります。つまり、

長所：分散が低下し、滑らかで再現性の高いスペクトルになる
短所：区間長が短くなるため周波数分解能は低下する（分散と分解能のトレードオフ）

また、切り出し部分時系列には、Hann窓をかけています。私の個人的意見としては、確率過程のサンプル時系列解析では、矩形窓の切り出しでも全く問題ないと思います。
　以下が、Rスクリプト例です。両側PSDを計算し半分だけ描画する仕様ですので、厳密な片側PSDではないことに注意してください。

n.sub <- 5   # Welch 法で平均化に用いる区間数
# ============================================================
# Welch 法によるパワースペクトル推定（50% オーバーラップ）
#
# 【入力として仮定する変数】
#   f.samp : サンプリング周波数 [Hz]
#   time   : 時刻ベクトル（長さ N）※ 本スクリプトでは直接は使用しない
#   x      : 時系列データ（長さ N）
#   n.sub  : 分割する区間数（Welch 平均の回数）
#
# 【出力】
#   freq.h : 非負周波数側（0 ～ ナイキスト周波数）の周波数軸
#   Pwelch : Welch 法で平均化したパワースペクトル密度
#            （両側スペクトルのスケーリングのまま，半分のみを描画）
# ============================================================

# ------------------------------------------------------------
# 入力の簡単なチェック
# ------------------------------------------------------------
stopifnot(is.numeric(f.samp), length(f.samp) == 1, f.samp > 0)
stopifnot(is.numeric(n.sub), length(n.sub) == 1, n.sub >= 1)
stopifnot(length(x) >= 4)

# ------------------------------------------------------------
# 直流成分（平均値）の除去
#   f = 0 成分の過大な影響を避けるため，事前に平均との差分を取る
# ------------------------------------------------------------
x0 <- x - mean(x, na.rm = TRUE)
N  <- length(x0)

# ------------------------------------------------------------
# 区間長 L の決定
#
# 50% オーバーラップの場合，
#   区間の開始点の間隔は step = L/2
#
# n.sub 個の区間を配置するための条件は
#   N >= L + (n.sub - 1) * (L / 2)
#      = L * (n.sub + 1) / 2
#
# よって，
#   L <= 2N / (n.sub + 1)
# ------------------------------------------------------------
L <- floor(2 * N / (n.sub + 1))

# 50% オーバーラップを整数サンプルで実現するため，
# 区間長 L は偶数にしておく
if (L %% 2 == 1) L <- L - 1

stopifnot(L >= 4)  # 区間長が極端に短くならないことを確認

# 区間の開始点と終了点
step   <- L / 2
starts <- 1 + (0:(n.sub - 1)) * step
ends   <- starts + L - 1

# データ長が不足していないかを確認
if (max(ends) > N) {
  stop("データ長が不足しています。n.sub を減らすか，より長いデータを使用してください。")
}

# ------------------------------------------------------------
# 窓関数の設定（Welch 法の標準：Hann 窓）
#
# 窓を掛けることで区間端での不連続を緩和し，
# スペクトルリークを抑える
# ------------------------------------------------------------
n <- 0:(L - 1)
w <- 0.5 - 0.5 * cos(2 * pi * n / (L - 1))  # Hann 窓
U <- sum(w^2)                               # 窓のエネルギー（正規化用）

# ------------------------------------------------------------
# この区間長 L に対応する周波数軸（両側）
# ------------------------------------------------------------
freq <- (0:(L - 1)) * (f.samp / L)

# ------------------------------------------------------------
# 各区間のペリオドグラムを加算していくための初期化
# ------------------------------------------------------------
Pacc <- rep(0, L)

# ------------------------------------------------------------
# Welch 法の本体：
#   1) 各区間を切り出す
#   2) 窓を掛ける
#   3) FFT を計算
#   4) 修正ペリオドグラムを計算
#   5) 全区間で平均
# ------------------------------------------------------------
for (i in seq_len(n.sub)) {

  # 区間の切り出し
  seg <- x0[starts[i]:ends[i]]

  # 窓掛け
  segw <- seg * w

  # FFT（R の fft() は正規化されていない）
  X <- fft(segw)

  # 修正ペリオドグラム
  #   S(f) = |X|^2 / (f.samp * sum(w^2))
  # とすることで，窓によるエネルギー変化を補正
  Pseg <- (Mod(X)^2) / (f.samp * U)

  # 区間ごとのスペクトルを加算
  Pacc <- Pacc + Pseg
}

# ------------------------------------------------------------
# Welch 平均（両側スペクトル）
# ------------------------------------------------------------
Pwelch2 <- Pacc / n.sub

# ------------------------------------------------------------
# 非負周波数側（0 ～ ナイキスト）のみを取り出す
# ※ one-sided スペクトルのような 2 倍補正は行わない
# ------------------------------------------------------------
k_max  <- floor(L / 2)
idx.h  <- 1:(k_max + 1)
freq.h <- freq[idx.h]
Pwelch <- Pwelch2[idx.h]

# ------------------------------------------------------------
# 結果の描画
# ------------------------------------------------------------
ymax <- max(Pwelch, na.rm = TRUE)
plot(freq.h, Pwelch,
     type = "l", col = 2, lwd = 2,
     xlab = "Frequency [Hz]",
     ylab = "PSD",
     main = sprintf("Welch PSD (n.sub = %d, 50%% overlap, Hann window)", n.sub),
     ylim = c(0, ymax * 1.2),
     xaxs = "i", yaxs = "i")

Welch 法で推定したパワースペクトル。時系列は「3. 付録」に示したRスクリプトを使用して事前に生成。

1.3 自己相関のフーリエ変換（コレログラム法）

　Wiener–Khinchin の定理により、定常過程では PSD は自己相関関数のフーリエ変換として定義されます。この定義に従い、

自己相関（自己共分散）を有限ラグまで推定
推定された自己相関にラグ窓を掛けて安定化
それをフーリエ変換して PSD を得る

という手順でPSDを推定します。

長所：定義に忠実で、理論（相関構造）と対応が明確
短所：ラグが大きい自己相関推定は不安定なので、有限ラグ打ち切り・ラグ窓が実質的に必須
備考：非循環自己相関を用いる限り、ペリオドグラムと厳密一致はしません（推定法が異なるため）

注意すべき点として、実際の有限長データに対しては、コレログラム法とペリオドグラムは本質的に大きく異なる推定法ではありません。本稿では自己相関推定に対して平滑化を行っていませんが、自己相関に三角窓（Bartlett 窓）などのラグ窓を掛けて高ラグ成分の減衰を強めると、周波数領域では平滑化が強く働き、結果としてペリオドグラムとの差は大きくなります。
　以下が、Rスクリプト例です。両側PSDを計算し半分だけ描画する仕様ですので、厳密な片側PSDではないことに注意してください。

# ============================================================
# 自己相関のフーリエ変換による PSD 推定
# （コレログラム法：non-circular correlogram）
#
# 【入力として仮定する変数】
#   f.samp : サンプリング周波数 [Hz]
#   time   : 時刻ベクトル（長さ N）※ 本スクリプトでは直接は使用しない
#   x      : 時系列データ（長さ N）
#
# 【出力】
#   freq.h : 非負周波数側（0 ～ ナイキスト）の周波数軸
#   Pcor   : 非循環自己相関 → FFT による PSD
# ============================================================

stopifnot(is.numeric(f.samp), length(f.samp) == 1, f.samp > 0)
stopifnot(length(x) >= 8)

# ------------------------------------------------------------
# 前処理：直流成分（平均値）を除去
# ------------------------------------------------------------
x0 <- x - mean(x, na.rm = TRUE)
N  <- length(x0)

# ------------------------------------------------------------
# 非循環自己相関（自己共分散）の推定
# 定義：
#   r[k] = Σ_{n=0}^{N-k-1} x[n] x[n+k] / N
# ------------------------------------------------------------
lag.max <- floor(N / 4)   # 有限ラグで打ち切る
r.pos <- numeric(lag.max + 1)

for (k in 0:lag.max) {
  r.pos[k + 1] <- sum(x0[1:(N - k)] * x0[(1 + k):N]) / N
}

# ------------------------------------------------------------
# 偶関数として負ラグ側を復元
#   [r[0], r[1], ..., r[K], r[K-1], ..., r[1]]
# ------------------------------------------------------------
r.two <- c(r.pos, rev(r.pos[-1]))
M <- length(r.two)

# ------------------------------------------------------------
# 自己相関のフーリエ変換（コレログラム）
# ------------------------------------------------------------
S2 <- Re(fft(r.two))

# ------------------------------------------------------------
# PSD としてのスケーリング
# r[k] を 1/N で定義しているため，
# 周波数密度に変換するために dt = 1 / f.samp を掛ける
# ------------------------------------------------------------
dt <- 1 / f.samp
Pcor2 <- S2 * dt

# ------------------------------------------------------------
# 周波数軸
#   周波数刻みは f.samp / M
# ------------------------------------------------------------
freq <- (0:(M - 1)) * (f.samp / M)

# ------------------------------------------------------------
# 非負周波数側（0 ～ ナイキスト）のみを描画
#   ※ one-sided の 2 倍補正は行わない
# ------------------------------------------------------------
k_max <- floor(M / 2)
idx.h <- 1:(k_max + 1)
freq.h <- freq[idx.h]
Pcor   <- Pcor2[idx.h]

# ------------------------------------------------------------
# 描画
# ------------------------------------------------------------
ymax <- max(Pcor, na.rm = TRUE)
plot(freq.h, Pcor, type = "l", col = 2, lwd = 2,
     xlab = "Frequency [Hz]", ylab = "PSD",
     main = "PSD from ACF (non-circular correlogram)",
     ylim = c(0, ymax * 1.2), xaxs = "i", yaxs = "i")

コレログラム法で推定したパワースペクトル。時系列は「3. 付録」に示したRスクリプトを使用して事前に生成。

2. パラメトリック法

　次に、パラメトリック推定法として、自己回帰（AR）モデルに基づく方法を、Rスクリプト例とともに示します。パラメトリック法では、時系列がある確率モデルに従うと仮定し、その有限個のパラメタ（例：AR 係数と雑音分散）を推定した上で、モデルが与える理論スペクトルを計算します。
　ノンパラメトリック法が「データから直接」PSD を作るのに対し、パラメトリック法は「モデルを介して」PSD を求める点が本質的に異なります。
　ここでは、AR スペクトル推定の代表的手法として

Yule–Walker 法（自己相関方程式）
Burg 法（反射係数を用いた推定）
AR 次数の選択（AIC など）を併用した AR スペクトル推定

の 3 つを扱います。いずれも「両側スペクトル（two-sided）」として計算し、非負周波数側（0からNyquist周波数まで）だけを描画します（one-sided 化の 2 倍補正は行いません）。
　パラメトリック法は、 * 短いデータでもピーク周波数を鋭く推定したい場合に有効、 * 一方で、推定結果は AR次数 に大きく依存し、次数を大きくすると偽ピークが現れやすくなる、という特徴があります。
　真のスペクトル構造が事前に分からない実データ解析では、ノンパラメトリック法・パラメトリック法のいずれか一つに依存するのではなく、複数の推定結果を比較しながら、整合的で妥当な解釈を行うことが重要です。

2.1 ARモデルによるPSD（Yule–Walker 法）

　AR( $\displaystyle{p}$ ) モデルは、時系列 $\displaystyle{ x[t] }$ が過去 $\displaystyle{ p }$ 点の線形結合と白色雑音で記述できると仮定するモデルです。

$\displaystyle{ x[t] = \sum_{k=1}^{p} a_k x[t-k] + \varepsilon[t], \qquad \varepsilon[t] \sim (0,\sigma^2) }$

このとき、AR モデルの理論的 PSD は

$\displaystyle{ S(f) = \frac{\sigma^2}{\left|1-\sum_{k=1}^{p} a_k e^{-i2\pi f k/f_s}\right|^2} }$

で与えられます。Yule–Walker 法は、推定した自己相関（あるいは共分散）から AR 係数を求める方法です。

長所：計算が安定で実装が簡単、古典的で理解しやすい
短所：次数 p の選び方に結果が強く依存、ピークが過剰に立つ場合がある

　以下が R スクリプト例です。両側 PSD を計算し半分だけ描画する仕様ですので、厳密な片側 PSD ではないことに注意してください。ARモデルの次数 p の値を大きくしたり小さくしたりして、推定結果の変化を観察してみてください。

# ============================================================
# ARモデル（Yule–Walker 法）によるPSD推定
#
# 【入力として仮定する変数】
#   f.samp : サンプリング周波数 [Hz]
#   time   : 時刻ベクトル（長さ N）※ 本スクリプトでは直接は使用しない
#   x      : 時系列データ（長さ N）
#
# 【設定】
#   p      : AR次数（ユーザが指定）
#
# 【出力】
#   freq.h : 非負周波数側（0 ～ ナイキスト）の周波数軸
#   Par    : ARモデル由来のPSD（両側定義のまま半分を描画）
# ============================================================

p <- 15  # AR次数（例）

stopifnot(is.numeric(f.samp), length(f.samp) == 1, f.samp > 0)
stopifnot(length(x) >= (p + 5))

# 前処理：平均除去（DC成分除去）
x0 <- x - mean(x, na.rm = TRUE)
N  <- length(x0)

# Yule–Walker 法でAR係数を推定
fit <- ar.yw(x0, order.max = p, aic = FALSE)

a     <- as.numeric(fit$ar)          # a1..ap
sigma <- as.numeric(fit$var.pred)    # 雑音分散 σ^2（予測誤差分散）

# 周波数軸（N点：ペリオドグラムと同じ刻みで評価）
freq <- (0:(N - 1)) * (f.samp / N)

# ARスペクトルの計算
#   A(e^{-iω}) = 1 - Σ a_k e^{-iωk},  ω = 2π f / f.samp
omega <- 2 * pi * freq / f.samp
den <- rep(0+0i, N)
for (k in 1:p) den <- den + a[k] * exp(-1i * omega * k)
H <- 1 - den
Par2 <- sigma / (Mod(H)^2)           # 両側PSD

# 非負周波数側（0..Nyquist）のみ
k_max <- floor(N / 2)
idx.h <- 1:(k_max + 1)
freq.h <- freq[idx.h]
Par <- Par2[idx.h]

# 描画
ymax <- max(Par, na.rm = TRUE)
plot(freq.h, Par, type = "l", col = 2, lwd = 2,
     xlab = "Frequency [Hz]", ylab = "PSD",
     main = sprintf("AR PSD (Yule–Walker, p=%d; two-sided half)", p),
     ylim = c(0, ymax * 1.2), xaxs = "i", yaxs = "i")

ARモデルの当てはめを使って推定したパワースペクトル。時系列は「3. 付録」に示したRスクリプトを使用して事前に生成。

2.2 ARモデルによるPSD（Burg 法）

　Burg 法は、反射係数（PARCOR）を用いて AR 係数を推定する方法で、短いデータでもピーク位置が鋭く推定されやすいことからスペクトル推定でよく用いられます。Rでは ar.burg() が利用できます。

長所：短いデータでも分解能が高くなりやすい、推定が比較的安定
短所：次数 (p) が大きいと偽ピーク（スパイク）が出やすい

　以下が R スクリプト例です。

# ============================================================
# ARモデル（Burg 法）によるPSD推定
# ============================================================

p <- 15  # AR次数（例）

stopifnot(is.numeric(f.samp), length(f.samp) == 1, f.samp > 0)
stopifnot(length(x) >= (p + 5))

x0 <- x - mean(x, na.rm = TRUE)
N  <- length(x0)

# Burg 法でAR係数を推定
fit <- ar.burg(x0, order.max = p)

a     <- as.numeric(fit$ar)
sigma <- as.numeric(fit$var.pred)

# 周波数軸（N点）
freq <- (0:(N - 1)) * (f.samp / N)
omega <- 2 * pi * freq / f.samp

# ARスペクトルの計算
den <- rep(0+0i, N)
for (k in 1:p) den <- den + a[k] * exp(-1i * omega * k)
H <- 1 - den
Par2 <- sigma / (Mod(H)^2)

# 非負周波数側のみ
k_max <- floor(N / 2)
idx.h <- 1:(k_max + 1)
freq.h <- freq[idx.h]
Par <- Par2[idx.h]

# 描画
ymax <- max(Par, na.rm = TRUE)
plot(freq.h, Par, type = "l", col = 2, lwd = 2,
     xlab = "Frequency [Hz]", ylab = "PSD",
     main = sprintf("AR PSD (Burg, p=%d; two-sided half)", p),
     ylim = c(0, ymax * 1.2), xaxs = "i", yaxs = "i")

2.3 AR次数の選択（AIC）を併用したARスペクトル推定

　パラメトリック法では 次数 $\displaystyle{p}$ の選択が結果を左右します。小さ過ぎるとピークが潰れ、大き過ぎると偽ピークが現れます。そこで、情報量規準（AIC など）で次数を自動選択し、その次数で AR PSD を計算する方法が実務的です。
　Rの ar()（method="yw" や "burg"）は、AIC に基づく次数選択を行えます。

長所：次数選択の恣意性を減らせる、再現性が高い
短所：AIC の最適次数が必ずしも目的（ピーク推定など）に最適とは限らない

　以下が R スクリプト例です（Burg 法＋AIC を例示）。

# ============================================================
# AIC による次数選択を併用した AR PSD 推定（例：Burg）
# ============================================================

stopifnot(is.numeric(f.samp), length(f.samp) == 1, f.samp > 0)
stopifnot(length(x) >= 16)

x0 <- x - mean(x, na.rm = TRUE)
N  <- length(x0)

# ar() による次数選択（AIC最小の次数を採用）
# method は "yw"（Yule–Walker）や "burg" を選択可能
fit <- ar(x0, method = "burg", aic = TRUE)

p     <- fit$order
a     <- as.numeric(fit$ar)
sigma <- as.numeric(fit$var.pred)

cat(sprintf("Selected AR order (AIC): p = %d\n", p))

# 周波数軸（N点）
freq <- (0:(N - 1)) * (f.samp / N)
omega <- 2 * pi * freq / f.samp

# ARスペクトルの計算
den <- rep(0+0i, N)
for (k in 1:p) den <- den + a[k] * exp(-1i * omega * k)
H <- 1 - den
Par2 <- sigma / (Mod(H)^2)

# 非負周波数側のみ
k_max <- floor(N / 2)
idx.h <- 1:(k_max + 1)
freq.h <- freq[idx.h]
Par <- Par2[idx.h]

# 描画
ymax <- max(Par, na.rm = TRUE)
plot(freq.h, Par, type = "l", col = 2, lwd = 2,
     xlab = "Frequency [Hz]", ylab = "PSD",
     main = sprintf("AR PSD (Burg + AIC, p=%d; two-sided half)", p),
     ylim = c(0, ymax * 1.2), xaxs = "i", yaxs = "i")

AIC を用いて AR モデルの次数を決定し、推定された AR モデルのパラメタに基づいて算出したパワースペクトル。時系列は「3. 付録」に示したRスクリプトを使用して事前に生成。

3. 付録：心拍変動のLF、HF成分シミュレーション

　本付録では、心拍変動（HRV）に典型的な LF（低周波）成分および HF（高周波）成分を含む時系列データを、簡易的な数理モデルにより生成するための R スクリプトを示します。本シミュレーションは、スペクトル解析手法の挙動を確認するための擬似データ生成を目的としており、生理学的メカニズムの厳密な再現を意図したものではありません。
　HRV における LF 成分および HF 成分は、それぞれ

LF：おおよそ10秒周期 (0.1 HZ) の振動をカバーする 0.04–0.15 Hz（血圧反射などに関連）
HF：おおよそ4秒周期 (0.25 HZ) の振動をカバーする 0.15–0.40 Hz（呼吸性変動）

に対応する振動成分として知られています。ここでは、AR(2) モデルを用いて、指定した中心周波数にピークを持つ振動成分を個別に生成し、それらを線形に合成することで HRV 時系列を構成します。

シミュレーションのパラメタ設定

　関数 hrv_lfhf() では、LF 成分と HF 成分を独立に生成し、それらを合成して HRV の擬似時系列を作ります。主なパラメタは以下のとおりです。

fs：サンプリング周波数 [Hz] （例：4 Hz は補間後の HRV 解析で一般的）
minutes：信号長 [分]
f_lf, f_hf：LF / HF 成分の中心周波数 [Hz]
r_lf, r_hf：LF / HF 成分の極半径
- 大きいほどスペクトルピークが鋭くなります
LFHF：指定する LF/HF 比（振幅ベース）
- LF 成分の振幅を $\displaystyle{\sqrt{\mathrm{LFHF}}}$ 倍することで制御しています
sd_ms：最終的な HRV 時系列の標準偏差 [ms]

ここでの LF/HF 比は、解析結果として得られる厳密なパワー比を保証するものではなく、あくまで「LF 成分と HF 成分の相対的な強さを制御するための簡易パラメタ」です。

出力と利用方法

関数の出力はリスト形式で、

t：時間軸 [s]
x_ms：HRV タコグラム [ms]
LFHF：設定した LF/HF 比

などを含みます。本章で示した各種スペクトル推定法（ペリオドグラム、Welch 法、コレログラム法、AR スペクトル推定）の動作確認や比較に、そのまま利用できます。

############################################################
# 心拍変動（HRV）の LF / HF 成分シミュレーション
#
# 【目的】
#  ・LF（低周波）/ HF（高周波）帯域に対応する振動成分を持つ
#    HRV 時系列データを生成する
#  ・LF/HF 比をパラメタとして直接指定できる簡易モデルを構成する
#
# 【基本方針】
#  ・LF 成分、HF 成分をそれぞれ AR(2) モデルで生成
#  ・両者を線形結合し、HRV タコグラムを構成
############################################################


# ==========================================================
# AR(2) モデルによる単一ピーク振動成分の生成
#
# モデル：
#   x[n] = 2 r cos(2π f0 / fs) x[n-1]
#        - r^2 x[n-2] + e[n]
#
# 極が複素共役対となるため、周波数 f0 に鋭いピークを持つ
# ==========================================================
sim_ar2 <- function(N, fs, f0, r = 0.98, burn = 2000) {

  # 角周波数
  w <- 2 * pi * f0 / fs

  # AR(2) 係数
  a1 <- 2 * r * cos(w)
  a2 <- -r^2

  # 初期過渡を捨てるために余分に生成
  M <- N + burn

  # 配列の初期化
  x <- numeric(M)
  e <- rnorm(M)     # 白色雑音

  # 初期値
  x[1] <- e[1]
  x[2] <- a1 * x[1] + e[2]

  # AR(2) 再帰計算
  for (n in 3:M) {
    x[n] <- a1 * x[n - 1] + a2 * x[n - 2] + e[n]
  }

  # 初期過渡を除去し、平均値を 0 に揃える
  x <- x[(burn + 1):M]
  x - mean(x)
}


# ==========================================================
# HRV（LF + HF）成分を含む時系列の生成
#
# 【引数】
#  fs      : サンプリング周波数 [Hz]
#  minutes : 信号長 [分]
#  f_lf    : LF 成分の中心周波数 [Hz]
#  f_hf    : HF 成分の中心周波数 [Hz]
#  r_lf    : LF 成分の極半径（ピークの鋭さ）
#  r_hf    : HF 成分の極半径（ピークの鋭さ）
#  LFHF    : 指定する LF/HF 比（振幅ベース）
#  sd_ms   : 出力 HRV の標準偏差 [ms]
# ==========================================================
hrv_lfhf <- function(fs = 4, minutes = 5,
                     f_lf = 0.10, f_hf = 0.25,
                     r_lf = 0.992, r_hf = 0.985,
                     LFHF = 1.0, sd_ms = 50) {

  # データ点数
  N <- fs * 60 * minutes

  # 時間軸 [s]
  t <- (0:(N - 1)) / fs

  # --------------------------------------------------------
  # LF / HF 成分をそれぞれ独立に生成
  # --------------------------------------------------------
  x_lf <- sim_ar2(N, fs, f_lf, r_lf)
  x_hf <- sim_ar2(N, fs, f_hf, r_hf)

  # --------------------------------------------------------
  # LF/HF 比の制御
  #   振幅比として LF を sqrt(LFHF) 倍
  #   （パワー比がおおよそ LFHF になる）
  # --------------------------------------------------------
  g_lf <- sqrt(LFHF)
  g_hf <- 1

  # LF + HF の合成
  x <- g_lf * x_lf + g_hf * x_hf

  # --------------------------------------------------------
  # 振幅スケーリング
  #   平均 0、標準偏差 sd_ms [ms] に正規化
  # --------------------------------------------------------
  x <- as.numeric(scale(x)) * sd_ms

  # 結果をリストで返す
  list(
    fs   = fs,    # サンプリング周波数 [Hz]
    t    = t,     # 時間軸 [s]
    x_ms = x,     # HRV タコグラム [ms]
    LFHF = LFHF   # 設定した LF/HF 比
  )
}


# ==========================================================
# 使用例
# ==========================================================

# サンプリング周波数
f.samp <- 4

# HRV シミュレーションの実行
HRV.sim <- hrv_lfhf(
  fs      = f.samp,
  minutes = 5,
  LFHF    = 0.3,
  sd_ms   = 50
)

# 解析用変数として取り出し
time <- HRV.sim$t
x    <- HRV.sim$x_ms

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