【Rでスペクトル解析２】直線トレンドが作る「見かけの1/fスペクトル」

時系列解析においてパワースペクトルは、確率過程の特性を調べるための基本的な道具と言えます。確かに、「純粋な」確率過程から生成された時系列サンプルに対してスペクトル解析を行うのであれば、この理解で問題ありません。しかし、実際のデータにはそうでない構造、たとえば直線トレンドが含まれていることが多く、その場合には確率過程として単純に読み解くことができなくなります。
　このブログでは「1/f ゆらぎ」を頻繁に扱っていますが、実は直線トレンドを含んだ時系列に対して前処理を行わずにパワースペクトルを推定すると、パワースペクトル $S(f)$ の低周波数側に $S(f) \sim 1/f^ 2$ という構造が現れることがあります。 $S(f) \sim 1/f^ 2$ という形だけを見ると、「ブラウン運動のような拡散的な時系列だ」と誤って解釈してしまうかもしれません。
　そこで今回は、確率過程のパワースペクトルを正しく推定・解釈するためには、直線などのトレンド成分を事前に除去することが不可欠であるという点について解説します。

有限長さの直線的時系列のパワースペクトル。左は直線的に増加する時系列、右はそのパワースペクトル。

1. Rでパワースペクトルを計算してみる

　ここでは、単純な直線関数 $x(t)=t$ をサンプリング間隔 $dt$ (dt) で $N$ (N) 点抽出した場合のパワースペクトルをRで計算してみます．
　まず、有限長の直線を次のように作ります。観測時間は $T=N \,dt$ で、時間軸は $t _ n = n \, dt$ ( $n=0, 1, \dots, N-1$ )です。信号は $x(t _ n) = t _ n$ とします。パワースペクトルの推定は、最も基本的な方法として、FFTの振幅二乗に基づくピリオドグラム（periodogram）を使います。
　以下のスクリプトでは、

長さ N の有限長データを作る
FFTで周波数成分を計算する
片側パワースペクトル $S(f)$ を作る
log–log プロットで描く
傾きを比べるために $1/f^ 2$ のグラフを重ねる

という流れで計算しています。以下のスクリプトをコピーして、みなさんもRで実行してみてください。

############################################################
# ノイズなし：有限区間 [0, T] の直線 x(t)=t
# カラー表示・指数表記軸つきパワースペクトル
############################################################

# =========================
# 0) 基本設定
# =========================
N  <- 4096
dt <- 1.0
T  <- N * dt

t <- (0:(N-1)) * dt
x <- t    # ノイズなし直線

# =========================
# 1) パワースペクトル
# =========================
X  <- fft(x)
P2 <- (Mod(X)^2) * (dt / N)

n_half <- N / 2
f  <- (0:n_half) / (N * dt)
S1 <- P2[1:(n_half + 1)]
S1[2:n_half] <- 2 * S1[2:n_half]

# log 軸のため 0 Hz を除外
idx <- which(f > 0)

# =========================
# 2) 軸用の目盛（10^n）
# =========================
fx_ticks <- 10^(floor(log10(min(f[idx]))) : ceiling(log10(max(f[idx]))))
fy_ticks <- 10^(floor(log10(min(S1[idx]))) : ceiling(log10(max(S1[idx]))))

# =========================
# 3) プロット
# =========================
par(mar = c(5, 5, 3, 1), las = 1)

plot(f[idx], S1[idx],
     log  = "xy",
     type = "l",
     col  = "steelblue",
     lwd  = 4,
     xaxt = "n",
     yaxt = "n",
     xlab = "Frequency f",
     ylab = "Power spectrum S(f)",
     main = "Power spectrum of a finite-length linear trend")

# x 軸（10^n 表記）
axis(1, at = fx_ticks,
     labels = parse(text = paste0("10^", log10(fx_ticks))))

# y 軸（10^n 表記）
axis(2, at = fy_ticks,
     labels = parse(text = paste0("10^", log10(fy_ticks))))

box()

# =========================
# 4) 参照線：1/f^2
# =========================
k0 <- 20
f0 <- f[idx][k0]
S0 <- S1[idx][k0]
ref <- S0 * (f[idx] / f0)^(-2)

lines(f[idx], ref,
      col = "firebrick",
      lwd = 2,
      lty = 2)

legend("topright",
       legend = c("S(f) (linear trend)", "reference: 1/f^2"),
       col    = c("steelblue", "firebrick"),
       lwd    = 2,
       lty    = c(1, 2),
       bty    = "n")

　上の例では直線の“パワースペクトル”を計算しましたが、ホワイトノイズに直線トレンドが加わった時系列のパワースペクトルを計算する場合は、上のスクリプトの x を以下のように変更してください。

eps <- rnorm(N)
sig <- 100
x <- t + sig*eps    # 直線トレンド ＋ 白色ノイズ

この x は下の図のような時系列になります。

sig の値を変化させると、パワースペクトルは下の図のようになります。みなさんも、sig の値を増やして、パワースペクトルがどのようになるか体験してください。

直線トレンドを含むホワイトノイズのパワースペクトル。ノイズの強さを変化させた場合。赤破線は $1/f^ 2$ の傾き。

直線トレンドを含むホワイトノイズのパワースペクトル。ノイズの強さを変化させた場合。赤破線は $1/f^ 2$ の傾き。

2. 解析的にパワースペクトルを見積ってみる

　ここからは、先ほど数値的に確認した「直線トレンドが低周波側で $1/f^ 2$ のように見える」現象を、解析的に整理します。

2.1 大前提：直線は非定常で、厳密な意味のパワースペクトルは定義できません

　まず大前提を確認します。直線

$\displaystyle{ x(t)=a t }$

のような信号は、確定信号（決定論的信号）であり、確率過程ではありません。また、平均値が時間とともに変化するため定常でもありません。したがって、定常確率過程に対して成り立つ Wiener–Khinchin 定理（自己相関関数のフーリエ変換としてパワースペクトル密度 $S(f)$ を定義する枠組み）における、厳密な意味でのパワースペクトル密度は存在しません。

2.2 仮定：有限区間 $[0,T]$ の連続信号 $x(t) =a t$

　以下の連続信号を考えます。

$\displaystyle{ x(t)=a t,\qquad 0\le t\le T. }$

ここでは、「有限区間フーリエ変換」を

$\displaystyle{ F(f)=\int_{0}^{T} a t\, e^{-i2\pi f t}\,dt }$

で計算します。 $T \to \infty$ では、 $x$ が発散するので、有限区間でないとこの積分は計算できません。

2.3 フーリエ変換の計算

　定数 $a$ を外に出して

$\displaystyle{ F(f)=a\int_{0}^{T} t\, e^{-i2\pi f t}\,dt }$

を計算します。部分積分を用いて計算すると、

$\displaystyle{ \begin{aligned} F(f)&=\int_{0}^{T} a\, t\, e^{-i2\pi f t}\,dt \\ &= a \left( \left[\frac{t\,e^{-i2\pi f t}}{-i2\pi f}\right]_{t=0}^{t=T} -\int_{0}^{T}\frac{e^{-i2\pi f t}}{-i2\pi f}\,dt \right)\\ &=a\left( \frac{T e^{-i2\pi f T}}{-i\,2\pi f} +\frac{1-e^{-i 2\pi f T}}{(2\pi f)^2} \right). \end{aligned} }$

がえられます。

2.4 有限区間スペクトル

　ここでは便宜的に

$\displaystyle{ S(f)=|F(f)|^2 }$

を「パワースペクトルに相当する量」として扱います。区間幅 $T$ で割っていないため、これは「パワースペクトル密度」ではなく、有限区間信号のフーリエエネルギー分布のようなものです。

2.5 低周波側のスケーリング

　低周波数 $f\to 0$ では $\displaystyle{ 2\pi f T \ll 1 }$ が成り立つため、

$\displaystyle{ e^{-i2\pi f T} =1-i2\pi f T-\frac{(2\pi f T)^2}{2}+O(f^3) }$

と展開できます。これを上の $F(f)$ に代入して整理すると、支配的な項は

$\displaystyle{ F(f)\sim a\,\frac{T}{2\pi f} \qquad (f\to 0) }$

となります。したがって

$\displaystyle{ \begin{aligned} S(f)&=|F(f)|^2\\ &\sim a^2\frac{T^2}{(2\pi f)^2}\\ &\propto \frac{1}{f^2} \qquad (|2\pi f T| \ll 1). \end{aligned} }$

となります。このように、有限区間 $[0,T$ ] の直線トレンドをそのままフーリエ変換し、正規化せずに振幅二乗を取ると、低周波側に $1/f^ 2$ 型の構造が現れる**ことが分かります。

3. まとめ：直線トレンドはパワースペクトル推定の前に必ず除け

　今回は、直線トレンドを含む時系列に対してパワースペクトルをそのまま計算すると、低周波数側に $1/f^ 2$ 型の構造が現れてしまうことを、数値実験と解析計算の両面から確認しました。この結果が示しているのは、時系列データに対してパワースペクトルを推定する際には、直線トレンドや低次多項式トレンドといった確率的でない成分を事前に除去することが不可欠であるという点です。これらの前処理を行って初めて、パワースペクトルは「確率過程としての揺らぎの性質」を正しく反映したものとして解釈できるようになります。直線トレンドはパワースペクトル推定の「邪魔者」であり、解析の前段階で必ず取り除くべき成分です。

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