時系列のパワースペクトルの推定は、時系列に含まれる周期性を検出したり、背後にある確率過程の基本特性を調べたりするための、最も代表的な解析法の一つで、「スペクトル解析」と呼ばれます。今回はスペクトル解析の基礎を解説します。
1. 真のパワースペクトル
スペクトル解析を理解するうえで、最初に押さえておくべき重要な考え方があります。それは「真のパワースペクトル」というものが現象の背後に存在し、それを観測された時系列データから推定するのがスペクトル解析であるという立場です。
ここでいう「真のパワースペクトル」とは、ある確率過程が本来もっている周波数特性を表すパワースペクトルのことです。もし無限に長い時系列を観測でき、確率的な平均を正確に取ることができるなら、その過程のパワースペクトルは一意に定まります。弱定常性と呼ばれる、数学的に理想化され単純化された過程においては、パワースペクトルは自己相関関数のフーリエ変換として定義される量であり、これは確率過程の「DNA」とも言える生成ルールを表しています。
一方、私たちが実際に扱うのは、有限の長さをもつ単一の時系列データにすぎません。そこには観測ノイズや偶然の揺らぎが含まれ、確率平均を直接評価することはできません。したがって、実データからそのまま「真のパワースペクトル」を見ることは不可能です。それでもなお、限られた時系列データから背後にある真のパワースペクトルをできる限り妥当に推定すること——これこそがスペクトル解析の目的です。
下の図は、２次の自己回帰過程 AR(2) により生成された時系列データに対して、（左）spectrum() 関数をデフォルト設定のまま用いた推定結果と、（右）Welch 法による区間平均に平滑化を加えた推定結果を、それぞれ理論的なパワースペクトル（赤破線）と重ねて描いたものです。
spectrum() コマンドを単純に使用した結果、右はWelch法と平滑化というテクニックを使った結果。spectrum() による推定では、理論的には一つしか存在しないはずのピーク全体がギザギザしており、ピークの高さや幅が不明瞭です。これは推定誤差というより、確率過程の実現値がもつ確率的な特性が、そのまま可視化された結果です。真のパワースペクトルは、期待値をとらないと現れないのです。
一方、上の右の図では、同じデータを Welch 法で区間に分割し、各区間のスペクトルを平均したうえで平滑化を施しています。その結果、スペクトル全体の揺らぎは大きく抑えられ、理論的なパワースペクトルに対応する単一のピーク構造が明瞭に現れています。
この例が示しているのは、同じ時系列データであっても、推定手法によって「見えるスペクトル」は大きく異なるという事実です。現実のデータに対して、どの推定結果が「正しい」と単純に判断することはできません。重要なのは、どのような仮定のもとで推定が行われ、どの情報が強調され、どの情報が犠牲になっているのかを理解したうえで、結果を解釈することです。スペクトル解析とは、単にグラフを描く作業ではありません。有限なデータから本来は見えない構造を推測しようとする、試行錯誤の過程なのです。
【R】ここで使用したRスクリプト
############################################################
# AR(2)過程のパワースペクトル推定と真のパワースペクトル
# (1) spectrum() デフォルト
# (2) Welch（区間平均） + 平滑化（移動平均）
############################################################
set.seed(1) # 乱数を固定
# ---- sampling & AR(2) parameters ----
fs <- 10
N <- 4096
f0 <- 2.3
r <- 0.85
sigma <- 1.0
sigma2 <- sigma^2
# AR(2) coefficients (complex-conjugate poles)
w0 <- 2*pi*f0/fs
a1 <- 2*r*cos(w0)
a2 <- -r^2
# ---- simulate AR(2) ----
x <- as.numeric(arima.sim(model = list(ar = c(a1, a2)), n = N, sd = sigma))
# ==========================================================
# (1) spectrum() default
# ==========================================================
sp1 <- spectrum(x, plot = FALSE) # freq: cycles/sample
f1 <- sp1$freq * fs # Hz
S1 <- sp1$spec / fs # -> x^2/Hz
# theoretical PSD on the same grid (x^2/Hz)
z1 <- exp(-1i*2*pi*f1/fs)
z2 <- exp(-1i*4*pi*f1/fs)
S1_theory <- (sigma2 / Mod(1 - a1*z1 - a2*z2)^2) / fs
# ==========================================================
# (2) Welch PSD (segment-averaged) + smoothing
# ==========================================================
L <- 512
ov <- 0.5
st <- as.integer(L*(1-ov))
idx0 <- seq(1, N - L + 1, by = st)
Pacc <- NULL
for(i in idx0){
seg <- x[i:(i+L-1)]
pg <- spec.pgram(seg, taper = 0.1, detrend = TRUE, plot = FALSE, fast = TRUE)
if(is.null(Pacc)){
Pacc <- pg$spec
f2 <- pg$freq * fs
} else {
Pacc <- Pacc + pg$spec
}
}
Pwelch <- (Pacc / length(idx0)) / fs # -> x^2/Hz（Hz表示に合わせる）
# smoothing: moving average
m <- 15
w <- rep(1/m, m)
Pwelch_s <- as.numeric(stats::filter(Pwelch, w, sides = 2))
na <- is.na(Pwelch_s)
Pwelch_s[na] <- Pwelch[na]
# theoretical PSD on Welch grid (x^2/Hz)
z1b <- exp(-1i*2*pi*f2/fs)
z2b <- exp(-1i*4*pi*f2/fs)
S2_theory <- (sigma2 / Mod(1 - a1*z1b - a2*z2b)^2) / fs
# ==========================================================
# plot
# ==========================================================
op <- par(mfrow = c(1,2), mar = c(4,4,2,1), pty="m")
plot(f1, S1, type = "l", col = 4, xaxs = "i",
xlab = "Frequency [Hz]", ylab = "PSD [x^2/Hz]",
main = "Default spectrum(x)")
lines(f1, S1_theory, lty = 2, lwd = 2, col = 2)
legend("topright",
legend = c("estimated PSD", "theoretical PSD"),
lty = c(1,2), lwd = c(1,2), bty = "n", col = c(4,2))
plot(f2, Pwelch_s, type = "l", col = 4, xaxs = "i",
xlab = "Frequency [Hz]", ylab = "PSD [x^2/Hz]",
main = "Welch + smoothing")
lines(f2, S2_theory, lty = 2, lwd = 2, col = 2)
legend("topright",
legend = c("estimated PSD", "theoretical PSD"),
lty = c(1,2), lwd = c(1,2), bty = "n", col = c(4,2))
par(op)
2. spectrum() を使いこなす
R でスペクトル解析をするとき、お世話になるのが fft() と spectrum() です。どちらを使っても、パワースペクトルを推定できますが、結果は少し違います。ここでは、その理由を説明することで、spectrum() では何をしているのかを理解していきます。
2.1 まずは fft()：最も素朴なパワースペクトル推定
時系列 に対して FFT を計算し、その絶対値の二乗を取れば、いわゆるピリオドグラム（periodogram）と呼ばれるパワースペクトルの最も素朴な推定結果が得られます。
R の fft() を使ってパワースペクトルを計算するには以下のようにします。
# sampling frequency [Hz] fs <- 10 # FFT-based (naive) power spectrum N <- length(x) x0 <- x - mean(x) X <- fft(x0) df <- fs / N P2 <- (Mod(X)^2) / (N * fs) k <- 0:floor(N/2) # indices for 0..fs/2 f <- k * df P <- P2[k + 1]
ここで、まじめな実務家に向けて注意をしておきますが、この結果は「正式な」片側スペクトルになっていません。「正式な」片側スペクトルは、 とNyquist周波数以外の値を２倍します。今回は、
spectrum() が正式な片側スペクトルを与えないので、あえてまがい物の片側のみ描画するスペクトルを考えます。
fft() を使った単純なスペクトル推定と spectrum() を使った推定の比較。微妙に違うが、spectrum() が良いとも感じない。
2.2 spectrum() のデフォルトは前処理を含む
次に spectrum(x) を使うと、同じデータでも fft() で求めたピリオドグラムとは少し違う形になります。これはバグではなく、spectrum() が内部でいくつかの処理を自動的に行っているためです。fftを使った上の結果と一致させたければ、
spectrum(x, plot = FALSE, detrend=FALSE, taper=FALSE)
とします。detrend=FALSE, taper=FALSEを指定しない場合、spectrum() ではデフォルトで以下の前処理がされます。
(1) 直線トレンドの除去（detrend）
平均値だけでなく、緩やかな直線トレンドがあると、低周波成分が過大に評価されます。spectrum() はデフォルトで 直線トレンドを除去し、低周波の歪みを抑えます。fft() でパワースペクトルを計算するときは、Rで以下のようにすれば直線トレンド除けます。
# xに含まれる直線トレンドを除く t <- seq_along(x) x_detrended <- as.numeric(resid(lm(x ~ t)))
(2) 軽いテーパ（窓）の付与（taper）
有限長データをそのまま FFT すると、端点の不連続によりリーク（偽の周期成分）が生じます。spectrum() はデフォルトで、データの端に 弱いテーパ（窓）をかけ、リークを軽減します。fft() でパワースペクトルを計算するときは、Rで以下のようにすれば似た処理ができます。
# xの端を幅ゼロに絞る taper <- 0.1 n <- length(x); m <- floor(n*taper/2) w <- rep(1,n) if(m>0){ k <- seq_len(m) w[c(k,n-k+1)] <- c(0.5*(1-cos(pi*k/(m+1))), rev(0.5*(1-cos(pi*k/(m+1))))) } x_tapered <- as.numeric(x)*w
【R】ここで使用したRスクリプト
############################################################
# AR(2)過程のスペクトル解析
############################################################
set.seed(1)
# ---- sampling & AR(2) parameters ----
fs <- 10
N <- 4096
f0 <- 2.3
r <- 0.85
sigma <- 1.0
sigma2 <- sigma^2
w0 <- 2*pi*f0/fs
a1 <- 2*r*cos(w0)
a2 <- -r^2
# ---- simulate AR(2) ----
x <- as.numeric(arima.sim(model = list(ar = c(a1, a2)), n = N, sd = sigma))
# ==========================================================
# (1) FFT-based PSD: two-sided (x^2/Hz), plot 0..fs/2 only
# ==========================================================
x0 <- x - mean(x)
X <- fft(x0) # unnormalized DFT
df <- fs / N # Hz
P2 <- (Mod(X)^2) / (N * fs) # two-sided PSD [x^2/Hz]
k <- 0:floor(N/2) # indices for 0..fs/2
f <- k * df
P <- P2[k + 1] # NO doubling
# ==========================================================
# (2) spectrum() default (no options)
# ==========================================================
sp <- spectrum(x, plot = FALSE) # freq: cycles/sample (0..0.5)
f_sp <- sp$freq * fs # Hz
S_sp <- sp$spec / fs # convert density to per-Hz
# ==========================================================
# theoretical PSD in per-Hz (two-sided definition)
# S_f(f) = (1/fs) * sigma^2 / |1 - a1 e^{-j2π f/fs} - a2 e^{-j4π f/fs}|^2
# ==========================================================
S_theory_f <- function(f_hz){
z1 <- exp(-1i*2*pi*f_hz/fs)
z2 <- exp(-1i*4*pi*f_hz/fs)
(sigma2 / Mod(1 - a1*z1 - a2*z2)^2) / fs
}
S_th_fft <- S_theory_f(f)
S_th_sp <- S_theory_f(f_sp)
# ---- plot (two panels) ----
op <- par(mfrow = c(1,2), mar = c(4,4,2,1), pty = "s")
plot(f, P, type="l", col=4, xaxs="i",
xlab="Frequency [Hz]", ylab="PSD [x^2/Hz]",
main="FFT (two-sided PSD; plotted 0..fs/2)")
lines(f, S_th_fft, lty=2, lwd=2, col=2)
legend("topright",
legend=c("estimated PSD", "theoretical PSD"),
lty=c(1,2), lwd=c(1,2), bty="n", col=c(4,2))
plot(f_sp, S_sp, type="l", col=4, xaxs="i",
xlab="Frequency [Hz]", ylab="PSD [x^2/Hz]",
main="spectrum() default (converted to per-Hz)")
lines(f_sp, S_th_sp, lty=2, lwd=2, col=2)
legend("topright",
legend=c("estimated PSD", "theoretical PSD"),
lty=c(1,2), lwd=c(1,2), bty="n", col=c(4,2))
par(op)
2.3 spectrum() で平滑化処理オプションを指定
spectrum() では計算したパワースペクトルを平滑化するオプションが用意されています。それが、spans オプションです。たとえば、
spectrum(x, spans = c(9, 9, 9))
のようにします。spans は 奇数を並べた整数ベクトルです。大きな数字にすれば、平滑化するウィンドウの幅が広がります。そして、同じ値をたくさん並べれば並べるほど、平滑化の重みが正規分布に近づきます。平滑化に使う重みは以下のRスクリプトで計算できます。
spans <- "c(9 ,9 ,9)"
w0 <- kernel("modified.daniell",m=eval(parse(text = spans)))$coef
w <- c(rev(w0[-1]),w0)
m <- (length(w)-1)/2
j <- (-m):(m)
w.max <- max(w)
plot(j,w,ylim=c(0,w.max*1.1),pch=16,col=2,las=1,yaxs="i",main=paste("spans =",spans))
arrows(j, 0, j, w, col =2,length=0,lwd=2)
spans = c(9 ,9 ,9) のときに、平滑化に使う重み係数。上に示したスクリプトで計算できます。
【Rで時系列解析】スペクトルウインドーを使ったパワースペクトルの平滑化：Rのspectrumでのspansオプションの意味 - ケィオスの時系列解析メモランダム
ここでは、数値実験の結果のみ示しておきます。平滑化すると真のパワースペクトルに近い形が現れます。
spectrum()でオプション指定なしの場合、右はspectrum()でspans=c(9, 9, 9)を指定した場合。赤の破線が真のパワースペクトル。
【R】ここで使用したRスクリプト
############################################################
# AR(2): spectrum() default vs spectrum() with spans
# + theoretical PSD overlay
# (units: x^2/Hz; NO one-sided doubling)
############################################################
set.seed(1)
# ---- sampling & AR(2) parameters ----
fs <- 10
N <- 4096
f0 <- 2.3
r <- 0.85
sigma <- 1.0
sigma2 <- sigma^2
w0 <- 2*pi*f0/fs
a1 <- 2*r*cos(w0)
a2 <- -r^2
# ---- simulate AR(2) ----
x <- as.numeric(arima.sim(model = list(ar = c(a1, a2)), n = N, sd = sigma))
# ==========================================================
# spectrum() (A) default -> convert to per-Hz
# ==========================================================
sp_def <- spectrum(x, plot = FALSE) # freq: cycles/sample (0..0.5)
f_hz <- sp_def$freq * fs # Hz
S_def <- sp_def$spec / fs # per-Hz
# ==========================================================
# spectrum() (B) spans specified -> convert to per-Hz
# spans: odd integers; larger => stronger smoothing
# ==========================================================
spans_use <- c(9, 9, 9) # <-- adjust as needed
sp_spn <- spectrum(x, spans = spans_use, plot = FALSE)
S_spn <- sp_spn$spec / fs # per-Hz
# f grid is identical to default when N/taper/detrend are same:
# f_hz <- sp_spn$freq * fs
# ==========================================================
# theoretical PSD in per-Hz (two-sided definition)
# S_f(f) = (1/fs) * sigma^2 / |1 - a1 e^{-j2π f/fs} - a2 e^{-j4π f/fs}|^2
# ==========================================================
S_theory_f <- function(f_hz){
z1 <- exp(-1i*2*pi*f_hz/fs)
z2 <- exp(-1i*4*pi*f_hz/fs)
(sigma2 / Mod(1 - a1*z1 - a2*z2)^2) / fs
}
S_th <- S_theory_f(f_hz)
# ==========================================================
# plot: left=default, right=spans
# ==========================================================
op <- par(mfrow = c(1,2), mar = c(4,4,2,1), pty = "m")
plot(f_hz, S_def, type="l", col=4, xaxs="i",
xlab="Frequency [Hz]", ylab="PSD [x^2/Hz]",
main="spectrum() default")
lines(f_hz, S_th, lty=2, lwd=2, col=2)
legend("topright",
legend=c("estimated PSD", "theoretical PSD"),
lty=c(1,2), lwd=c(1,2), bty="n", col=c(4,2))
plot(f_hz, S_spn, type="l", col=4, xaxs="i",
xlab="Frequency [Hz]", ylab="PSD [x^2/Hz]",
main=sprintf("spectrum() with spans=c(%s)",
paste(spans_use, collapse=",")))
lines(f_hz, S_th, lty=2, lwd=2, col=2)
legend("topright",
legend=c("estimated PSD", "theoretical PSD"),
lty=c(1,2), lwd=c(1,2), bty="n", col=c(4,2))
par(op)
3. まとめ
パワースペクトルの推定はとても大切な技術なので、つづきを書いていきます。
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
