【時系列解析】自己回帰過程のパワースペクトルを解析的に求める

自己回帰過程のパワースペクトルを求めるテクニックの説明です．自己回帰過程の係数が与えられた状況から，パワースペクトルの式を導くことはとても簡単です．でも，大切なことは理屈を感じることです．

必要な基本知識

長さNの時系列 $\{x[i]\}_{i=0}^{N-1}$ のフーリエ変換とその逆変換は，

　 $\begin{eqnarray} x[i] &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \nonumber \\ X(f_k) &=& \sum_{i=0}^{N-1} x[i] \, e^{- {\rm \bf i} 2 \pi f_k i} \nonumber \end{eqnarray}$

です． $f_k = k/N$ で， ${\rm \bf i}$ は虚数単位です．

　 $\begin{eqnarray} x[i] &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \nonumber \end{eqnarray}$

は恒等式なので， $i$ に， $i -j$ を代入すると，

　 $\begin{eqnarray} x[i-j] &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k (i-j)} \nonumber \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \, e^{-{\rm \bf i} 2 \pi f_k j} \nonumber \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \, \left( e^{-{\rm \bf i} 2 \pi f_k} \right)^j \nonumber \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} \left\{X(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-j} \right\} e^{{\rm \bf i} 2 \pi f_k i} \nonumber \end{eqnarray}$

フーリエ変換の式と見比べてみると，

$x[i-j]$ のフーリエ変換は， $X(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-j}$

となってます．つまり， $x[i-j]$ の， $-j$ の部分に注目して， $\left( e^{{\rm \bf i} 2 \pi f_k} \right)$ の $-j$ 乗を， $X(f_k)$ にかけてあるだけです．この部分は，各周波数成分の波の位相が周波数 $f_k$ に依存してシフトすることを表しています． $i$ と $i-j$ では，位置がほんの少しずれているので，振幅は同じで，波の位置（位相）だけずれます．

練習すると，

$x[i-1]$ のフーリエ変換は， $X(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-1}$

$x[i-2]$ のフーリエ変換は， $X(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-2}$

$y[i]$ のフーリエ変換が $Y(f_k)$ なら， $y[i-m]$ のフーリエ変換は， $Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-m}$

最後に，長さNの時系列 $\{x[i] \}_{i=0}^{N-1}$ のフーリエ変換を， $X(f_k)$ とすれば，パワースペクトルの推定量 $\hat{S}_x (k_f)$ は，

$\displaystyle \hat{S}_x (f_k) = \frac{|X(f_k)|^2}{N}$

以上で準備は終わりです．

差分方程式の両辺をフーリエ変換してパワースペクトルを計算

ということで，上の知識を使えば差分方程式の両辺をフーリエ変換して，フーリエ変換の関係式を導けます．パワースペクトルも簡単に求められます．

例えば， $x[i]$ のフーリエ変換が $X(f_k)$ で， $y[i]$ のフーリエ変換が $Y(f_k)$ として，1次自己回帰過程のような差分方程式，

　 $y[i] = a y[i-1] + x [i]$

を考えます．この両辺をフーリエ変換すると，

　 $\begin{eqnarray} Y(f_k) &=& a Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-1} + X (f_k) \\ \left(1- a e^{-{\rm \bf i} 2 \pi f_k} \right) Y(f_k) &=& X(f_k) \end{eqnarray}$

両辺をパワースペクトルの形にするために，両辺の絶対値の２乗をとって，Nでわると，

　 $\begin{eqnarray} \frac{\left|\left(1- a e^{-{\rm \bf i} 2 \pi f_k} \right) Y(f_k)\right|^2}{N} &=& \frac{\left|X(f_k)\right|^2}{N} \\ \left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2 \frac{ \left|Y(f_k)\right|^2}{N} &=& \frac{\left|X(f_k)\right|^2}{N} \end{eqnarray}$

$x[i]$ のパワースペクトルの推定量は， $S_x(f_k) = \frac{ \left|X(f_k)\right|^2}{N}$ で，
$y[i]$ のパワースペクトルの推定量は， $S_y(f_k) = \frac{ \left|Y(f_k)\right|^2}{N}$ なので，
パワースペクトルの関係式

　 $\begin{eqnarray} \left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2 S_y(f_k) &=& S_x(f_k) \\ S_y(f_k) &=& \frac{1}{\left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2}S_x(f_k) \end{eqnarray}$

がえられます． $x [i]$ が白色ノイズ $w [i]$ （分散 $\sigma^2$ ）のときが，１次自己回帰過程

　 $y[i] = a y[i-1] + w [i]$

だったので， $w [i]$ のパワースペクトル $S_x [i] = S_w [i] = \sigma^2$ を上の式に代入すれば，

　 $\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{\left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2} \end{eqnarray}$

となります．複素数の絶対値の２乗が， $|z|^2 = z \bar{z}$ であることに注意して，オイラーの公式

　 $e^{{\rm \bf i} \theta} = \cos \theta + {\rm \bf i} \sin \theta$

を使えば，今まで何度もお世話になってきた，１次自己回帰過程のパワースペクトル

　 $\displaystyle S(f)=\frac{\sigma^2}{\left( 1- a e^{-{\rm \bf i} 2 \pi f_k} \right) \left( 1- a e^{ {\rm \bf i} 2 \pi f_k} \right)} = \frac{\sigma^2}{1+a^2-2 a \cos (2 \pi f)}$

になっていることが分かります．

白色ノイズの積分過程も計算してみる

数学的にパワースペクトルが存在するとか，しないとか，私はそんなこと気にしないので，白色ノイズの積分過程

　 $y[i] = y[i-1] + w [i]$

を考えます．両辺をフーリエ変換して，絶対値をとって，２乗して，Ｎで割れば，

　 $\begin{eqnarray} S_y(f_k) &=& \frac{1}{\left| 1- e^{-{\rm \bf i} 2 \pi f_k} \right|^2}S_w(f_k) \\ S_y(f_k) &=& \frac{\sigma^2}{2 - 2 \cos (2 \pi f)} \\ \end{eqnarray}$