離散フーリエ変換
これまで、離散的な時刻tk=kΔt (k=0, ±1, ±2, …, ±∞) でサンプリングされた無限個のデータに対するフーリエ変換Gs(f)を考え、その性質について論じてきた。
しかし、実際にコンピュータ上で数値計算を行う場合、データは有限個である。このとき、フーリエ変換は連続関数とはなり得ず、離散的な周波数fnに対する離散データとなる。これを計算するものが離散フーリエ変換(Discrete Fourier transform：DFT)である。
まず、DFTの定義を与える。
Nをデータ数とする。時刻の添え字をk、周波数の添え字をnとし、ともに0, 1, …, N-1の値をとるとする。このとき、離散データgkに対する離散フーリエ変換Gnを以下のように定義する。
$$G_n=\sum_{k=0}^{N-1}g_ke^{-\frac{2\pi ink}{N}} ：離散フーリエ変換$$
$$g_k=\frac{1}{N}\sum_{k=0}^{N-1}G_n^{\frac{2\pi ink}{N}} ：離散逆フーリエ変換$$
この定義と、元々のフーリエ変換・逆フーリエ変換の定義を見比べると、その類似性がわかる。
$$G(f)=\int_{-\infty}^{\infty}g(t)e^{-2\pi ift}dt$$
$$g(t)=\int_{-\infty}^{\infty}G(f)e^{2\pi ift}df$$
ここで注意すべきことは、g(t)およびG(f)は時刻tおよび周波数fの関数として与えられるのに対し、離散フーリエ変換に登場するgkおよびGnは単に添え字kとnによって決まる量に過ぎない点である。すなわち、GnからG(f)を知るためには
- 添え字kおよびnと時刻tおよび周波数fがどのように対応しているか
- gkおよびGnとg(t)およびG(f)がどのように対応しているか
を把握しておかないと結果を解釈することができない。
離散フーリエ変換とフーリエ変換の対応
離散フーリエ変換とフーリエ変換をどのように結び付けたらよいかを考える。
離散フーリエ変換の計算機能やプログラム自体は今日では簡単に利用できるが、その結果を正しく活用するための知識が必要である。
▶時刻0～Tのデータを得た場合
サンプリング周期をΔt、データ数をNとすると、T=NΔtの関係がある。このとき得られたN個のデータは、時刻
$$t_k=k\Delta t (k=0, 1, 2, …, N-1)$$
におけるデータg(tk)となる。
周波数fについては、0≦f≦fs (=1/Δtはサンプリング周波数)がN等分されてフーリエ変換が出力される。したがって、n番目の離散周波数fnは
$$f_n=\frac{n}{N\Delta t} (n=0, 1, 2, …, N-1)$$
となる。
これらを用いて、フーリエ変換式を近似する。
フーリエ変換式の離散化
\[
\begin{align*}
G(f_n) & ≃\sum_{k=0}^{N-1}g(t_k)e^{-2\pi if_nt_k}\Delta t \\
& =\Delta t\sum_{k=0}^{N-1}g(t_k)e^{-2\pi i\frac{n}{N\Delta t}k\Delta t} \\
& =\Delta t\sum_{k=0}^{N-1}g(t_k)e^{-\frac{2\pi ink}{N}}
\end{align*}
\]
ここで、時刻tkにおけるg(tk)をDFTへの入力データgkと同一とみなすと、
$$G(f_n)=\Delta t\sum_{k=0}^{N-1}g_ke^{-\frac{2\pi ink}{N}}$$
となる。Σ以降の部分はDFTの定義式そのものである。すなわち、周波数fnにおけるフーリエ変換G(fn)とDFTのn番目の出力Gnには以下の関係式がある。
$$G(f_n)=\Delta tG_n$$
$$t_k=k\Delta t (k=0, 1, 2, …, N-1)，f_n=\frac{n}{N\Delta t} (n=0, 1, 2, …, N-1)$$
DFTの出力Gnにサンプリング周波数Δtを掛けることでフーリエ変換G(fn)となることがわかる。
逆フーリエ変換式の離散化
\[
\begin{align*}
g(t_k) & ≃\sum_{n=0}^{N-1}G(f_n)e^{2\pi if_nt_k}\Delta f \\
& =\sum_{n=0}^{N-1}G(f_n)e^{2\pi i\frac{n}{N\Delta t}k\Delta t}\frac{1}{N\Delta t} \\
& =\frac{1}{N}\sum_{k=0}^{N-1}\frac{G(f_n)}{\Delta t}e^{\frac{2\pi ink}{N}}
\end{align*}
\]
ここで、G(fn)=ΔtGnを代入すると
$$g(t_k)=\frac{1}{N}\sum_{k=0}^{N-1}\frac{G(f_n)}{\Delta t}e^{\frac{2\pi ink}{N}}$$
となるが、右辺は離散逆フーリエ変換の定義式そのものである。
$$g(t_k)=g_k$$
これより、離散逆フーリエ変換の結果gkが時刻tkにおける元データの値g(tk)になることがわかる。
パワースペクトルの計算
DFTの出力Gnから周波数fnにおけるパワースペクトルを求める。
\[
\begin{align*}
P(f_n) & =\frac{1}{T}|G(f_n)|^2 \\
& =\frac{1}{N\Delta t}|\Delta tG_n|^2 \\
& =\frac{\Delta t}{N}|G_n|^2
\end{align*}
\]
▶時刻-T/2～T/2のデータを得た場合
参考のため、フーリエ変換とDFTの対応を計算してみる。
この場合、離散的な時刻tkは
$$t_k=-\frac{T}{2}+k\Delta t=-\frac{N\Delta t}{2}+k\Delta t (k=0, 1, 2, …, N-1)$$
である点が先の場合と異なる。これを用いてフーリエ変換を離散化すると
\[
\begin{align*}
G(f_n) & ≃\sum_{k=0}^{N-1}g(t_k)e^{-2\pi if_nt_k}\Delta t \\
& =\Delta t\sum_{k=0}^{N-1}g(t_k)e^{-2\pi i\frac{n}{N\Delta t}k\left(-\frac{N\Delta t}{2}+k\Delta t\right)} \\
& =\Delta t\sum_{k=0}^{N-1}g_ke^{-\frac{2\pi ink}{N}}e^{in\pi}
\end{align*}
\]
となる。ここで、exp(inπ)=(-1)^nと書けることから
\[
\begin{align*}
G(f_n) & \Delta t(-1)^n\sum_{k=0}^{N-1}g_ke^{-\frac{2\pi ink}{N}} \\
& =\Delta t(-1)^nG_n
\end{align*}
\]
となる。
