以前説明した「自己回帰 (AR)パワースペクトルの分解：波素分析」（詳細は以下のリンク参照）
自己回帰 (AR)パワースペクトルの分解：波素分析 - ケィオスの時系列解析メモランダム
をRで実際にやってみようという記事です。ARスペクトル分解（波素分析）は、ARモデルの数理とパワースペクトル推定、複素解析、代数定理などを理解していれば、すっきりと理解できるものなのですが、最近は数学嫌いな研究者が増えてなかなか使いこなせないようです。理論的なことはよくわからないけどとにかく試したいという方のためにRスクリプトを掲載しておきます。
専門家の方への注意ですが、私が過去に説明したARスペクトル分解（波素分析）と以下のRスクリプトでは、パワースペクトルの形状を連続時間・連続周波数近似しています。ですので、離散時系列のパワースペクトルとしては厳密には正しくありません。離散版のARスペクトル分解（波素分析）も構築可能ですので、興味がある方は自分で理論を導出してください。代数の基本定理や、部分分数分解がローラン展開と留数定理で計算できることがわかっていれば、結果は見通せると思います。離散過程のAR(1)とAR(2)の厳密な自己相関関数とパワースペクトルは私の過去の記事に書いてあると思います。
ARスペクトル分解（波素分析）への最短ルート
Rを使ったことがある人は、下のRスクリプトをコピーして実行してみてください。
############################################################
# AR スペクトル分解：波素分析
#
# 目的
# 1) Burg 法で AR(p) を推定
# 2) 推定 AR から 1-sided PSD P(f) を計算
# 3) AR の特性方程式の根（実根・複素根）に基づき、PSD を成分に分解
# 4) 成分和（psd_sum）が全体 PSD（psd_total）と一致することを確認
#
# 前提（ここでの AR の書き方）
# x(t) = Σ_{m=1..M} a_m x(t - mΔt) + ε(t), Var(ε)=σ^2
#
# 全体 PSD（片側ではなく連続の理論式の形）
# P(f) = (σ^2 Δt) / | 1 - Σ a_m exp(-i 2π f m Δt) |^2
#
# 分解の要点
# 特性方程式（多項式）の根 z を用いる（|z|<1 の安定根）
# - 実根 z>0: 減衰率 α = -log(z)/Δt のローレンツ型
# - 複素根: 減衰率 β = -log(|z|)/Δt と中心周波数 f0 = Arg(z)/(2πΔt)
# の共鳴型（±f0 を両側に持つ）
#
# 注意
# ・実根の z<=0 は log(z) が定義できないため、必要なら呼び出し側で
# 監視してください（例：負の実根が出たら警告）。
############################################################
# ==========================================================
# 1. AR スペクトル分解（根に基づく成分分解）
# ==========================================================
ar_spectral_decomposition <- function(x, p, fs = 1.0, n_freq = 1000) {
# ---- 基本設定 ----
dt <- 1 / fs # サンプリング間隔 Δt
freq_axis <- seq(0, fs/2, length.out = n_freq) # 0..ナイキストの片側周波数
# ---- AR 係数推定（Burg 法）----
# R の ar() が返す ar は (a1,...,aM) を想定（上のモデル式に対応）
ar_model <- ar(x, aic = FALSE, order.max = p, method = "burg")
a <- as.numeric(ar_model$ar) # a_1, ..., a_M
sigma2 <- ar_model$var.pred # 予測誤差分散 σ^2
# ---- 特性方程式の根を求める ----
# ここでは 1 - a1 z - a2 z^2 - ... - aM z^M = 0 の根を polyroot() で計算する。
# polyroot は「定数項→高次項」の順で係数を渡すため、
# 1 - a1 z - ... - aM z^M
# の係数ベクトルは c(1, -a1, -a2, ..., -aM)。
poly_coeffs <- c(-rev(a), 1)
z_roots <- polyroot(poly_coeffs)
# ---- D(z) を定義 ----
# D(z) = 1 / [ (1 - Σ a_m z^m) * (M - Σ a_m (M-m) z^{-m}) ]
calc_D <- function(z, a, M) {
m_seq <- 1:M
term1 <- 1 - sum(a * (z^m_seq))
term2 <- M - sum(a * (M - m_seq) * (z^(-m_seq)))
1 / (term1 * term2)
}
# ---- 分解スペクトル成分を格納（列を後から追加）----
psd_components <- matrix(0, nrow = n_freq, ncol = 0)
# ---- 根を実根・複素根に分類 ----
# 実根判定（数値誤差に対して許容幅を設定）
is_real <- abs(Im(z_roots)) < 1e-8
# ======================================================
# 1-A. 実根による成分（j = 1..m1）
# ======================================================
real_roots <- z_roots[is_real]
if (length(real_roots) > 0) {
for (zm in Re(real_roots)) {
# 実根 z から減衰率 α を定義： z = exp(-α Δt) → α = -log(z)/Δt
alpha_j <- -log(zm) / dt
# Fj = σ^2 Re[D(zm)]（実根の場合は本来 D が実数になる想定）
Fj <- sigma2 * Re(calc_D(zm, a, p))
# 成分 PSD: Pj(f) = (2 α Fj) / (α^2 + (2πf)^2)
Pj <- (2 * alpha_j * Fj) / (alpha_j^2 + 4 * pi^2 * freq_axis^2)
# 列として追加
psd_components <- cbind(psd_components, Pj)
}
}
# ======================================================
# 1-B. 複素根（Im>0 のみ採用）による成分（l = 1..m2）
# ======================================================
# 複素根は共役対で現れるので Im>0 側だけ取り、二重カウントを避ける。
comp_roots <- z_roots[!is_real & Im(z_roots) > 0]
if (length(comp_roots) > 0) {
for (zl in comp_roots) {
# 複素根 z = |z| exp(iθ) から
# |z| = exp(-β Δt) → β = -log(|z|)/Δt
# θ = 2π f0 Δt → f0 = Arg(z)/(2πΔt)
beta_l <- -log(Mod(zl)) / dt
fl <- Arg(zl) / (2 * pi * dt)
# D(z) から Gl, Hl を定義
Dz <- calc_D(zl, a, p)
Gl <- 2 * sigma2 * Re(Dz)
Hl <- 2 * sigma2 * Im(Dz)
# 成分 PSD（±fl の両側に共鳴を持つ形）
dp <- beta_l^2 + 4 * pi^2 * (freq_axis + fl)^2
dm <- beta_l^2 + 4 * pi^2 * (freq_axis - fl)^2
Pl <- Gl * (beta_l/dp + beta_l/dm) -
Hl * (2*pi*(freq_axis + fl)/dp - 2*pi*(freq_axis - fl)/dm)
# 列として追加
psd_components <- cbind(psd_components, Pl)
}
}
# ======================================================
# 1-C. 全体 PSD（直接計算）
# ======================================================
# P(f) = (σ^2 Δt) / | 1 - Σ a_m exp(-i 2π f m Δt) |^2
psd_total <- sapply(freq_axis, function(f) {
m_seq <- 1:p
(sigma2 * dt) /
Mod(1 - sum(a * exp(-complex(imaginary = 2*pi*f*m_seq*dt))))^2
})
# 分解成分の和（検算用）
psd_sum <- rowSums(psd_components)
return(list(
freq = freq_axis,
psd_total = psd_total,
psd_sum = psd_sum,
psd_components = psd_components,
opt_p = p
))
}
# ==========================================================
# 2. モデル生成と解析実行（例）
# ==========================================================
# ---- 例：5 次 AR 係数（ユーザー指定）----
a_5th <- c(1.3067, -1.1903, 1.1467, -0.9882, 0.3058)
# ---- 疑似データ生成 ----
# filter(..., method="recursive") は「再帰フィルタ」で AR 的系列を作る簡便な方法。
n_samples <- 2000
x <- filter(rnorm(n_samples), filter = a_5th, method = "recursive")
# ---- FPE による次数選択 ----
# FPE(M) = σ^2 * (1 + M/N) / (1 - M/N)
max_p <- 15
fpe_values <- numeric(max_p)
n <- length(x)
for (p in 1:max_p) {
model <- ar(x, aic = FALSE, order.max = p, method = "burg")
fpe_values[p] <- model$var.pred * (1 + p/n) / (1 - p/n)
}
opt_p <- which.min(fpe_values)
# ---- 分解実行 ----
res <- ar_spectral_decomposition(x, p = opt_p, fs = 1.0, n_freq = 1000)
# ==========================================================
# 3. プロット
# ==========================================================
# パネルのレイアウト設定
par(mfrow = c(1, 3), mar = c(5, 5, 3, 1),cex.lab=1.5,cex.axis=1.3, cex.main=1.7)
# -------------------------
# パネル 1: FPE プロット
# -------------------------
plot(1:max_p, fpe_values, type = "b", pch = 19, col = 4,
xlab = "Order M", ylab = "FPE(M)",
main = "1. Order Estimation (FPE)")
abline(v = opt_p, col = "red", lty = 2)
text(opt_p, max(fpe_values),
labels = paste("Best M =", opt_p),
pos = 4, col = "red")
# -------------------------
# パネル 2: 全体 PSD
# -------------------------
plot(res$freq, res$psd_total, type = "l", lwd = 2, col = 2,
xlab = "Frequency f", ylab = "P(f)",
main = "2. Power Spectral Density")
# -------------------------
# パネル 3: 分解成分の重ね描き
# -------------------------
# 全体スペクトルを背景として描画
plot(res$freq, res$psd_total, type = "l", lwd = 2, col = "gray85",
xlab = "Frequency f", ylab = "P(f)",
main = "3. Spectral Decomposition")
# 視認性の良い色（Set1 系）と透明度の設定
n_comp <- ncol(res$psd_components)
base_cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33")
if (n_comp > length(base_cols)) {
base_cols <- rainbow(n_comp)
}
# 半透明色（塗りつぶし用）と線色（境界線用）
poly_cols <- adjustcolor(base_cols[1:n_comp], alpha.f = 0.4)
line_cols <- base_cols[1:n_comp]
# 各成分を「半透明ポリゴン + 輪郭線」で重ねる
for (i in 1:n_comp) {
polygon(c(res$freq, rev(res$freq)),
c(res$psd_components[, i], rep(0, length(res$freq))),
col = poly_cols[i], border = NA)
lines(res$freq, res$psd_components[, i], col = line_cols[i], lwd = 1)
}
# 成分和（検算）：赤の点線
lines(res$freq, res$psd_sum, col = 2, lty = 2, lwd = 2)
# 凡例
legend("topright",
legend = c("Total PSD", "Sum of AR components"),
col = c("gray85", 2),
lwd = c(2, 2), lty = c(1, 2),
cex = 1.2, bty = "n")
このスクリプトを実行すれば、下図が描かれます。
■ ５次ARモデルで生成された時系列の解析
上のRスクリプトでは、既知の AR モデルから人工的に生成した時系列データを用いて、AR スペクトル分解を試しています。
a_5th <- c(1.3067, -1.1903, 1.1467, -0.9882, 0.3058) n_samples <- 2000 x <- filter(rnorm(n_samples), filter = a_5th, method = "recursive")
ここでは、5 次の AR モデル
を定義しています。ここで
：時刻
の時系列値
：平均 0、分散
の白色雑音
を表します。
このモデルから生成された時系列に対し、最終予測誤差 (FPE) を計算し、AR 次数を未知として次数 を推定します。
max_p <- 15 fpe_values <- numeric(max_p) n <- length(x) for (p in 1:max_p) { model <- ar(x, aic = FALSE, order.max = p, method = "burg") fpe_values[p] <- model$var.pred * (1 + p/n) / (1 - p/n) } opt_p <- which.min(fpe_values)
この部分では、FPE が最小となる次数を opt_p ()としています。Rスクリプトを実行して出力される図の左側が、FPEの次数
依存性です。
AR スペクトル分解は、選ばれた次数 opt_p を用いて、以下の関数を呼び出すことで実行できます。
res <- ar_spectral_decomposition(x, p = opt_p, fs = 1.0, n_freq = 1000)
この処理によって得られる主な結果は次の 3 つです。
res$psd_total推定された AR モデルから直接計算した 全体のパワースペクトルres$psd_components特性方程式の各根に対応する スペクトル成分res$psd_sum分解成分をすべて足し合わせたもの
出力される図の真ん中が、時系列にARモデルをフィットして推定された一般的なパワースペクトルです。
右側の図には、
- 灰色の太線：全体の PSD
- 半透明の色付き領域：分解された各ARスペクトル成分
- 赤の点線：成分の和（
psd_sum）
が描かれています。AR スペクトルの分解結果から、2 つの振動成分と 1 つの指数的な減衰成分の存在が読み取れます。
自分のデータを分析してみる
自分が用意したデータを分析してみるには、上のRスクリプトの
# ---- 疑似データ生成 ---- # filter(..., method="recursive") は「再帰フィルタ」で AR 的系列を作る簡便な方法。 n_samples <- 2000 x <- filter(rnorm(n_samples), filter = a_5th, method = "recursive")
の部分を削除し、xに自分のデータを代入してください。そして、実際のサンプリング周波数を以下のコマンド部分に設定してください。
res <- ar_spectral_decomposition(x, p = opt_p, fs = サンプリング周波数, n_freq = 1000)
たとえば、サンプリング周波数 fs が、32 Hzのときは、
res <- ar_spectral_decomposition(x, p = opt_p, fs = 32, n_freq = 1000)
とします。
■ 心拍変動への適用例
心拍変動時系列にARスペクトル分解（波素分析）を適用した例をのせておきます。ここでは、以下のようなRスクリプトを使いました。
############################################################
# 4-panel figure:
# Top (1 panel): HRV time series during Head-Up Tilt Test
# Bottom (3 panels): AR spectral decomposition (FPE / PSD / Components)
############################################################
DIR <- r"(ファイルのあるフォルダへのパスを指定)"
setwd(DIR)
options(digits.secs = 3)
FN <- "読み込むファイル名を指定"
TMP <- read.csv(FN, header = TRUE)
N <- nrow(TMP)
TMP$time <- as.POSIXct(TMP$time, tz = "Asia/Tokyo")
# Resampling at 2 Hz
fs <- 2
time.resamp <- seq(TMP$time[1], TMP$time[N], by = 1/fs)
RRI.resamp <- approx(TMP$time, TMP$RRI, xout = time.resamp)$y
# ==========================================================
# 1. AR spectral decomposition (root-based)
# ==========================================================
ar_spectral_decomposition <- function(x, p, fs = 1.0, n_freq = 1000) {
dt <- 1 / fs
freq_axis <- seq(0, fs/2, length.out = n_freq)
# AR coeffs by Burg
ar_model <- ar(x, aic = FALSE, order.max = p, method = "burg")
a <- as.numeric(ar_model$ar)
sigma2 <- ar_model$var.pred
# Roots of: 1 - a1 z - a2 z^2 - ... - aM z^M = 0
poly_coeffs <- c(-rev(a), 1)
z_roots <- polyroot(poly_coeffs)
calc_D <- function(z, a, M) {
m_seq <- 1:M
term1 <- 1 - sum(a * (z^m_seq))
term2 <- M - sum(a * (M - m_seq) * (z^(-m_seq)))
1 / (term1 * term2)
}
psd_components <- matrix(0, nrow = n_freq, ncol = 0)
is_real <- abs(Im(z_roots)) < 1e-8
# ---- real roots ----
real_roots <- z_roots[is_real]
if (length(real_roots) > 0) {
for (zm in Re(real_roots)) {
# NOTE: If zm <= 0, log(zm) is undefined; skip or warn.
if (zm <= 0) next
alpha_j <- -log(zm) / dt
Fj <- sigma2 * Re(calc_D(zm, a, p))
Pj <- (2 * alpha_j * Fj) / (alpha_j^2 + 4 * pi^2 * freq_axis^2)
psd_components <- cbind(psd_components, Pj)
}
}
# ---- complex roots (Im>0 only) ----
comp_roots <- z_roots[!is_real & Im(z_roots) > 0]
if (length(comp_roots) > 0) {
for (zl in comp_roots) {
beta_l <- -log(Mod(zl)) / dt
fl <- Arg(zl) / (2 * pi * dt)
Dz <- calc_D(zl, a, p)
Gl <- 2 * sigma2 * Re(Dz)
Hl <- 2 * sigma2 * Im(Dz)
dp <- beta_l^2 + 4 * pi^2 * (freq_axis + fl)^2
dm <- beta_l^2 + 4 * pi^2 * (freq_axis - fl)^2
Pl <- Gl * (beta_l/dp + beta_l/dm) -
Hl * (2*pi*(freq_axis + fl)/dp - 2*pi*(freq_axis - fl)/dm)
psd_components <- cbind(psd_components, Pl)
}
}
# ---- total PSD ----
psd_total <- sapply(freq_axis, function(f) {
m_seq <- 1:p
(sigma2 * dt) /
Mod(1 - sum(a * exp(-complex(imaginary = 2*pi*f*m_seq*dt))))^2
})
psd_sum <- rowSums(psd_components)
list(
freq = freq_axis,
psd_total = psd_total,
psd_sum = psd_sum,
psd_components = psd_components,
opt_p = p
)
}
# ==========================================================
# 2. Prepare data and select AR order by FPE
# ==========================================================
f.samp <- 2
x <- RRI.resamp
max_p <- 20
fpe_values <- numeric(max_p)
n <- length(x)
for (p in 1:max_p) {
model <- ar(x, aic = FALSE, order.max = p, method = "burg")
fpe_values[p] <- model$var.pred * (1 + p/n) / (1 - p/n)
}
opt_p <- which.min(fpe_values)
res <- ar_spectral_decomposition(x, p = opt_p, fs = f.samp, n_freq = 1000)
# ==========================================================
# 3. 4-panel plot layout (top=1, bottom=3)
# ==========================================================
op <- par(no.readonly = TRUE)
on.exit(par(op))
# Layout: 2 rows x 3 cols
# - Top row spans all 3 columns
# - Bottom row has 3 separate panels
layout(matrix(c(1, 1, 1,
2, 3, 4),
nrow = 2, byrow = TRUE),
heights = c(1.1, 1.0))
# Keep your bottom-panel design settings
par(mar = c(5, 5, 3, 1), cex.lab = 1.5, cex.axis = 1.3, cex.main = 1.7)
# -------------------------
# Panel 1 (Top): HRV time series
# -------------------------
plot(
time.resamp, RRI.resamp,
type = "l",
lwd = 2,
col = 2,
xaxs = "i",
xlab = "Time",
ylab = "R–R interval (ms)",
main = "(a) HRV during Head-Up Tilt Test"
)
# -------------------------
# Panel 2 (Bottom-left): FPE
# -------------------------
plot(1:max_p, fpe_values, type = "b", pch = 19, col = 4,
xaxs = "i",
xlab = "Order M", ylab = "FPE(M)",
main = "(b) Order Estimation (FPE)")
abline(v = opt_p, col = 2, lty = 2, lwd = 2)
text(opt_p, max(fpe_values),
labels = paste("Best M =", opt_p),
pos = 4, col = 2)
# -------------------------
# Panel 3 (Bottom-middle): Total PSD
# -------------------------
plot(res$freq, res$psd_total, type = "l", lwd = 2, col = 2,
xaxs = "i",
xlab = "Frequency f", ylab = "P(f)",
main = "(c) Power Spectral Density")
# -------------------------
# Panel 4 (Bottom-right): Components overlay
# -------------------------
plot(res$freq, res$psd_total, type = "l", lwd = 2, col = "gray85",
xaxs = "i",
xlab = "Frequency f", ylab = "P(f)",
main = "(d) Spectral Decomposition")
n_comp <- ncol(res$psd_components)
base_cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33")
if (n_comp > length(base_cols)) base_cols <- rainbow(n_comp)
poly_cols <- adjustcolor(base_cols[1:n_comp], alpha.f = 0.4)
line_cols <- base_cols[1:n_comp]
for (i in 1:n_comp) {
polygon(c(res$freq, rev(res$freq)),
c(res$psd_components[, i], rep(0, length(res$freq))),
col = poly_cols[i], border = NA)
lines(res$freq, res$psd_components[, i], col = line_cols[i], lwd = 1)
}
# Sum of components (check) in the requested style: lwd=2, col=2
lines(res$freq, res$psd_sum, col = 2, lty = 2, lwd = 2)
legend("topright",
legend = c("Total PSD", "Sum of AR components"),
col = c("gray85", 2),
lwd = c(2, 2), lty = c(1, 2),
cex = 1.2, bty = "n")
解析結果の例は以下です。
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
