Multitaper method の直観的な説明・その３

前回の続き。

有限長のデータを使う限り、スペクトルの「漏れ」は避けられず、その漏れ方は窓関数の関数形によって決まるということを見た。
これを踏まえ multitaper method (MTM) では、漏れに対する許容値をあらかじめパラメータとして与えた上で、「その許容値の範囲内で漏れが最小になるような窓関数を構成する」ことでスペクトル推定の精度を上げるというアプローチをとる。

データ長 N の時系列データ $X_t \quad (0 \le t \le N-1)$ を考えよう(簡単のためサンプリング間隔は1とする)。これをそのまま離散フーリエ変換すると、周波数解像度は $(2 N)^{-1}$ となる。言い換えれば、周波数領域では $(2 N)^{-1}$ Hz ごとにデータ点が存在する。

スペクトル漏れに対する許容値 W を、この解像度を単位として指定することにする。つまり、漏れの許容値を $W=1/(2N)$ と指定したら、それは周波数領域において、あるデータ点で表わされるスペクトルが隣のデータ点に漏れることは許容する（がそれ以上の漏れは許さない）ということで、 $W=2/(2N)$ ならさらにその隣に漏れてもよいとする、ということ。

さてここで、仮に許容値 W を完璧に実現するような窓関数が存在するとしよう。そのような窓関数のスペクトル分布は、周波数 f が $- W \lt f \lt W$ の範囲にあるときのみパワーを持ち、それ以外では 0 になっているはずだ。有限長の窓関数でそのようなスペクトル分布にできるだけ近い分布を持つものを求めるために、窓関数を N 次元のベクトル $\mathbf{w} = (w_1, w_2, ..., w_N)$ で表わし、この窓関数の周波数範囲 $- W \lt f \lt W$ でのパワー
$\int_{-W}^{W} |U(f)|^2 df \quad\quad \bigl(\quad U(f) = \sum_{t=1}^{N} w_t \exp (-2 \pi i f t)\quad \bigr)$
を最大化することで $\mathbf{w}$ を決めることを考える。

窓関数の規格化を制約条件として Lagrange の未定乗数法を使うと、
$\frac{\partial}{\partial \mathbf{w}} [ \int_{-W}^{W} |U(f)|^2 df - \lambda (|\mathbf{w}|^2 - 1) ] =0$
の解が求める窓関数となる。これは手計算で簡単に
$\sum_{t'=1}^N \frac{\sin [2 \pi W (t-t') ]}{\pi (t-t')} w_{t'} = \lambda w_t$
となることが示せるが、ここで $A_{tt'}=\frac{\sin [2 \pi W (t-t') ]}{\pi (t-t')}$ として行列 A を定義すると、この式は $A\mathbf{w}=\lambda \mathbf{w}$ という行列 A に関する固有値問題の解が、求める窓関数であることを表していることがわかる。

長くなったので続きは次回。