【時系列解析】サンプリング間隔がΔtの離散時系列のパワースペクトル

今回は，サンプリング間隔が1ではなく， $\Delta t$ のときのパワースペクトルについての説明です．スペクトル解析や時系列解析についての教科書では，サンプリング間隔が1の場合のみ扱っているものが多いので，迷うことがあると思います．

サンプリング間隔 $\Delta t$ の長さ $N$ の時系列
　 $\left\{x(0), x(\Delta t), x(2 \Delta t), \cdots, x\left((N-1)\Delta t \right) \right\}$
のパワースペクトル密度の推定量 $\hat{S}(f)$ （ピリオドグラム）は，
　時刻を $t_n = n \Delta t$ ，時間の長さを $T = N \Delta t$ ，周波数を $\displaystyle f_k = \frac{k}{N\Delta t}$
として，
　 $\displaystyle \hat{S}(f_k) = \frac{\Delta t }{N} \left|\sum_{n=0}^{N-1} x(t_n) \, e^{- i 2 \pi f_k t_n} \right|^2$
です．この式の，絶対値の中身は，
　 $\displaystyle \sum_{n=0}^{N-1} x (t_n) \, e^{- i 2 \pi f_k t_n} = \sum_{n=0}^{N-1} x (t_n) \, e^{- i 2 \pi \frac{k}{N} n}$
となるので，サンプリング間隔のことは考えずに，離散フーリエ変換したものと同じです．

つまり，実際の計算では，サンプリング間隔が1のときと同じ方法でまずパワースペクトル密度を推定し，その後，周波数を $\Delta t$ で割り，パワースペクトル密度に $\Delta t$ をかければよいのです．

また，サンプリング周波数を $f_{\rm samp}$ とすれば， $f_{\rm samp} = \frac{1}{\Delta t}$ なので， $f_{\rm samp}$ で考えるときは，周波数に $f_{\rm samp}$ をかけ，パワースペクトル密度を $f_{\rm samp}$ でわってください．

上の $\hat{S}(f)$ は，有限の長さの連続信号 $\left\{ x(t) \right\}_{0 \le t \le T}$ のパワースペクトル密度の推定量
　 $\hat{S}(f) = \frac{1}{T} \left|\int_{0}^{T} x(t) \, e^{- i 2 \pi f t} dt \right|^2$
を $\Delta t$ 幅の部分区間の和として近似することで導くことができます．
　 $\begin{eqnarray} \hat{S}(f_k) &=& \frac{1}{N \Delta t} \left|\sum_{n=0}^{N-1} x(t_n) \, e^{- i 2 \pi f_k t_n} \Delta t \right|^2 \nonumber \\ &=& \frac{\Delta t^2}{N \Delta t} \left|\sum_{n=0}^{N-1} x(t_n) \, e^{- i 2 \pi f_k t_n} \right|^2 \nonumber \\ &=& \frac{\Delta t }{N} \left|\sum_{n=0}^{N-1} x(t_n) \, e^{- i 2 \pi f_k t_n} \right|^2 \end{eqnarray}$

Rのspectrumを使う場合

Rでパワースペクトル密度を計算するコマンドspectrumを使う場合，サンプリング間隔をdtとすれば（以下は，dtが0.1 sの例です），

dt <- 0.1
x <- spectrum(x,plot=FALSE)
x$freq <- x$freq/dt
x$spec <- x$spec*dt

のようにします．このとき，x$freqが周波数で，x$specはパワースペクトル密度です．ただし，spectrumでは，直線トレンドを除いたり，窓関数をかけたりする処理がされるので，パワースペクトル密度の全面積は，上の定義の値と一致しません．

サンプリング周波数をf.sampとすれば（以下は，f.sampが10 Hzの例です），

f.samp <- 10
x <- spectrum(x,plot=FALSE)
x$freq <- x$freq*f.samp
x$spec <- x$spec/f.samp

です．

サンプリング間隔が $\Delta t$ のときの離散フーリエ変換

離散フーリエ変換についても， $\Delta t$ の扱いを説明します．ポイントは，フーリエ変換のときは， $\Delta t$ をかけ，逆フーリエ変換のときは， $\Delta t$ でわるということです．

離散時系列
　 $\left\{x(0), x(\Delta t), x(2 \Delta t), \cdots, x\left((N-1)\Delta t \right) \right\}$
について， $x_n=x(n \Delta t)$ ， $t_n = n \Delta t$ ， $\displaystyle f_k = \frac{k}{N\Delta t}$ とします．

まず，サンプリング間隔が1の，離散時系列$\left\{x_n \right\}_{n=0}^{N-1}$の離散フーリエ変換と逆変換は，それぞれ，
$\begin{eqnarray} X_k &=& \sum_{n=0}^{N-1} x_n \, e^{- i 2 \pi \frac{k}{N}n} \label{X_k} \\ x_n &=& \frac{1}{N} \sum_{k=0}^{N-1} X_k \, e^{ i 2 \pi \frac{k}{N}n} \label{x_k} \end{eqnarray}$
であることを念頭においてください．

連続信号の離散化として，サンプリング間隔 $\Delta t$ の場合考えると，
　 $\begin{eqnarray} X \left(f_k \right) &=& \sum_{n=0}^{N-1} x (t_n) \, e^{- i 2 \pi f_k t_n} \Delta t \\ &=& \Delta t \sum_{n=0}^{N-1} x (n \Delta t) \, e^{- i 2 \pi \frac{k}{N \Delta t} (n \Delta t)} \nonumber \\ &=& \Delta t \sum_{n=0}^{N-1} x_n \, e^{- i 2 \pi \frac{k}{N} n} = \Delta t \, X_k \nonumber \\ x (t_n) &=& \sum_{k=0}^{N-1} X \left(f_k \right) \, e^{ i 2 \pi f_k t_n } \Delta f \\ &=& \frac{1}{N \Delta t} \sum_{k=0}^{N-1} X_k \, e^{ i 2 \pi \frac{k}{N} n} \end{eqnarray}$
となりますので，フーリエ変換のときは， $\Delta t$ をかけ，逆フーリエ変換のときは， $\Delta t$ でわるという関係になることが分かります．