%
% An Introduction to Atmospheric Gravity Waves (2002) 
% 1st edition
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%             Style  Setting             %%%%%%%%
% フォント: 12point (最大), 片面印刷
\documentclass[a4j,12pt,oneside]{ltjreport}

\input{preemble}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%             Text Start                 %%%%%%%%
\begin{document}

\chapter{スペクトル解析}

\section{直接法}

離散的で有限な区間で定義された関数$ x(t_i) $を考える.
離散時間間隔を$ \Delta t $, データの個数を$ N $とすると,
基本周期$ T=N\Delta t $となる.
独立変数は$ t_i\equiv(i-1)\Delta t $の点で定義される.
これのフーリエ級数展開は
\begin{align}
	x(t_i)
	=&\frac{a_0}{2}
	+\sum_{n=1}^{m}\qty(
	a_n\cos\frac{2\pi}{T}nt_i
	+b_n\sin\frac{2\pi}{T}nt_i)
	\notag
	\\
	&(1\le m \le N/2),
	\\
	a_n
	=&\frac{2}{T}\int_{0}^{T}x(t_i)\cos\frac{2\pi}{T}nt_idt
	\notag
	\\
	=&\frac{2\Delta t}{T}\sum_{i=1}^{N}x(t_i)\cos\frac{2\pi}{T}nt_i,
	\\
	b_n
	=&\frac{2}{T}\int_{0}^{T}x(t_i)\sin\frac{2\pi}{T}nt_idt
	\notag
	\\
	=&\frac{2\Delta t}{T}\sum_{i=1}^{N}x(t_i)\sin\frac{2\pi}{T}nt_i
\end{align}
となる.
ただし, $ m $を$ N/2 $まで取ったときは$ b_{N/2}\equiv 0 $である.
ここから解像できる最小の周期は$ 2\Delta t $であることがわかる.
この周波数$ f_N\equiv 1/(2\Delta t) $をナイキスト周波数という.
ある周波数$ f_n=1/n\Delta t $の変動の強さは振幅
$ X(f_n)=\sqrt{a_n^2+b_n^2} $で表される.
その周波数の持つ変動の強さを表したものを\textbf{スペクトル(spectrum)}
という.
一般的には取扱が簡単な$ |X(f_n)| $を二乗した量,
すなわち
\begin{equation}\label{def_spect.eq}
	|X(f_n)|^2=a_n^2+b_n^2
\end{equation}
がスペクトルとして用いられる場合が多い.
直接法はデータを直接フーリエ変換し, スペクトルを求める方法である.
現在では計算の高速化のために高速フーリエ変換(FFT)
を用いて行われる場合が多い.
また誤差を減少させるために平滑化を行ったりテーパーと呼ばれる
重みをつけたりする場合もある.

\section{Blackman-Tukey 法}

Blackman-Tukey 法は自己共分散関数のフーリエ変換によって
パワースペクトルを推定する方法である.
「自己共分散関数のフーリエ変換はスペクトルであり,
スペクトルの逆フーリエ変換は自己共分散関数である」
という Wiener-Khintchine の関係に基づいている.
この関係から, 自己共分散関数, あるいはその規格化された関数
である自己相関関数で示される情報はすべてスペクトル中に
存在しているといえる.

まず無限のデータがあると考える. これを$ x(t_i) $とする.
このとき自己共分散関数は以下のように定義される:
\begin{equation}\label{kyoubunsan-func}
	C(\tau)\equiv
	\lim\limits_{N\rightarrow\infty}
	\frac{1}{N}
	\sum_{i=1}^{N}x(t_i)x(t_i+\tau).
\end{equation}
このとき$ C(0) $は定義から$ x(t_i) $の分散であり,
$ r(\tau)\equiv C(\tau)/C(0) $はラグ相関である.
自己共分散関数の定義から, $ C(\tau)=C(-\tau) $,
すなわち偶関数であること, $ |C(\tau)|\le C(0) $であることは自明である.

次にパワースペクトル密度関数$ S(f) $を定義する.
それは先程述べたように自己共分散関数のフーリエ変換として定義される:
\begin{align}
	C(\tau)
	&=\int_{-\infty}^{\infty}
		S(f)e^{i2\pi f\tau}df
	=2\int_{0}^{\infty}
		S(f)\cos(2\pi f\tau)df,
	\label{CtoS}
\\
	S(f)
	&=\int_{-\infty}^{\infty}
		C(\tau)e^{-i2\pi f\tau}d\tau
	=2\int_{0}^{\infty}
		C(\tau)\cos(2\pi f\tau)d\tau.
	\label{StoC}
\end{align}
このように自己共分散関数とパワースペクトルは互いにフーリエ変換の
関係にある.
ただし普通に言うパワースペクトル密度関数は
この$ S(f) $を2倍して, 正の周波数のみで定義された片側スペクトル
\begin{equation}\label{def_power}
	G(f)=2S(f)
	\;\;\;\;\;
	(0\le f<\infty)
\end{equation}
で定義されるものである.
$ G(f) $を使うと, \eqref{CtoS}, \eqref{StoC}は次のようになる:
\begin{align}
C(\tau)
&=\int_{0}^{\infty}
G(f)\cos(2\pi f\tau)df,
\label{CtoG}
\\
G(f)
&=4\int_{0}^{\infty}
C(\tau)\cos(2\pi f\tau)d\tau
\;\;\;\;\;
(0\le f<\infty).
\label{GtoC}
\end{align}
また
\begin{equation}\label{bunsan_C0}
	C(0)=\int_{0}^{\infty}G(f)df
\end{equation}
となるから, パワースペクトル密度関数の積分が分散に等しい.

次にデータレコードが有限である場合を考える.
すなわちデータを$ x_i\ (i=1,2,\cdots,N) $として
時間間隔を$ \Delta t $とする.
ラグが$ l(\ge 0) $のとき, 自己共分散関数は
\begin{equation}\label{kyoubunsan_risanka}
	\hat{C}(l)
	=\frac{1}{N-l}
	\sum_{i=1}^{N-l}x_i x_{i+l}
\end{equation}
と書ける.
$ l $の定義域は, $ l $に負を許して$ 0\le |l|\le M< N/2 $
であり($ -M\le l \le M $), 最大周期は$ 2M\Delta t $である.
次に\eqref{GtoC}を離散化すると
\begin{align}\label{BTpower}
	\hat{G}(k)&=4\qty[
			\frac{1}{2}\hat{C}(0)
			+\sum_{l=1}^{M-1}\hat{C}(l)
			 \cos\qty(2\pi\frac{k}{2M\Delta t}l\Delta t)
			+\frac{1}{2}\hat{C}(M)
			\cos\qty(2\pi\frac{k}{2M\Delta t}M\Delta t)
		]\Delta t
	\notag
	\\
	&=4\qty[
			\frac{1}{2}\hat{C}(0)
			+\sum_{l=1}^{M-1}\hat{C}(l)
			\cos\qty(\frac{\pi kl}{M})
			+\frac{1}{2}\hat{C}(M)
			\cos\qty(\pi k)
		]\Delta t
	\notag
	\\
	&=4\qty[
		\frac{1}{2}\hat{C}(0)
		+\sum_{l=1}^{M-1}\hat{C}(l)
		\cos\qty(\frac{\pi kl}{M})
		+\frac{1}{2}\hat{C}(M)(-1)^k
	]\Delta t
	\notag
	\\
	&=4\Delta t\qty[
		\sum_{l=0}^{M}\hat{C}(l)
		\cos\qty(\frac{\pi kl}{M})\cdot\delta_l]
	\notag
	\\
	&\mathrm{for}\ k=0,1,2,\cdots,M
\end{align}
ここで
\begin{align*}
	\delta_l
	=\begin{cases}
	1 & \mathrm{for}\ l=1,2,\cdots,M-1
	\\
	1/2 & \mathrm{for}\ l=0,M
	\end{cases}
\end{align*}
である.
$ \Delta f $は最大周期が$ 2M\Delta t $であることより
$ \Delta f =1/(2M\Delta t)=f_N/M $を使っている.
スペクトルは周波数$ f_k=k/(2M\Delta t)=kf_N/M\ 
\mathrm{for}\ k=0,1,2,\cdots,M$の$ M+1 $個で見積もっていることになる.
\eqref{BTpower}を見ると$ \hat{G}(k) $は多くの成分の加重平均である
事がわかる.
ここからこの方法によるパワースペクトルは安定したスペクトルの
推定を与えることがわかる.

実際にデータ処理を行う場合は適切な\textbf{ラグ・ウィンドウ(lag window)}
またはスペクトル・ウィンドウを用いる必要がある.
有限長さ$ P $のデータからの共分散関数の計算は
$ |\tau| \le \tau_m < P $の範囲で打ち切られる.
これは無限ラグの自己共分散関数に
\begin{align}
W_0(\tau)
=\begin{cases}
1 & (|\tau| \le \tau_m )
\\
0 & (|\tau| > \tau_m )
\end{cases}
\end{align}
のような箱型の関数を掛けていることと等価である.
これをラグ・ウィンドウと呼ぶ.
スペクトル・ウィンドウはこのラグ・ウィンドウのフーリエ変換である.
有限のラグで得られた共分散関数を$ C_0(\tau)\ (|\tau| \le \tau_m ) $
とし, $ \tau $を無限まで取った共分散関数を
$ C(\tau) $で表すとき,
これを式で表現すると
\begin{equation}\label{C_C_0W.eq}
C_0(\tau)=W_0(\tau)C(\tau)
\end{equation}
となる.
これをフーリエ変換するとスペクトルが得られるが, 
手続きは前と同じである. すなわち
$ C_0(\tau),\ W_0(\tau),\ C(\tau) $のフーリエ変換を
それぞれ$ S_0(f),\ Q_0(f),\ S(f) $とすると,
\eqref{C_C_0W.eq}のフーリエ変換は
\begin{equation}\label{FT_C0WC.eq}
	S_0(f)=
	\int_{-\infty}^{\infty}
	Q_0(f-f')S(f')df'
\end{equation}
となる.
この$ Q_0(f) $がラグ・ウィンドウのフーリエ変換なので,
ここでのスペクトル・ウィンドウである.
$ S_0(f) $は真のスペクトル$ S(f) $とスペクトル・ウィンドウ
$ Q_0(f) $のたたみ込みとなっている.
したがって$ Q_0(f) $は$ f=0 $付近でのみ値を持ち,
ほかはゼロになっていることが望ましい.
特に連続的なデータはデルタ関数,
離散的なデータでは
\begin{align}
Q_{00}(f)
=\begin{cases}
1/\Delta f=\tau_m & (|f| \le 0.5\Delta f )
\\
0 & (|f| > 0.5\Delta f )
\end{cases}
\end{align}
であれば全く歪みがない.
しかし実際の$ Q_0(f) $は
\begin{equation}\label{Q_0.eq}
	Q_0(f)
	=\int_{-\infty}^{\infty}
	W_0(\tau)e^{-i2\pi f\tau}d\tau
	=\int_{-\tau_m}^{\tau_m}
	e^{-i2\pi f\tau}d\tau
	=2\tau_m\frac{\sin(2\pi f\tau_m)}{2\pi f \tau_m}
\end{equation}
であり, 長い裾部を持ち, その一部が負の値となる.
これは真のスペクトルを大きく歪ませ, 強いスペクトルのピークの周辺に
負のスペクトルを与えることになる.

このようなことから, 望ましいラグ・ウィンドウが色々と提案されており,
Hanning ウィンドウや Hamming ウィンドウがよく使われる.
$ |\tau| $が大きくなるに連れて徐々に減衰する関数となっている.
すなわち適当なウェイトを自己共分散関数にかけて
平滑化することに対応している.


\end{document}

