フーリエ変換再{n}入門?

はじめに

先日(といっても1週間くらい前のことですが)、

実験で取得したデータをFFTで変換し、その結果得られた周波数領域のデータの各点の「周波数軸」方向の間隔って、どう計算すればよいかわかりますか?

という趣旨の質問を受けたのですが、即答することができず、次の日に回答するようなこととなってしまいました。お恥ずかしい限りです。( ;´･ω･`)

FFT,DFT,DCTなどのフーリエ変換方面は今の仕事では必ずしも専門ではないのですが、道具として使う可能性があるところで仕事をしていることと、質問されて翌日に答えても役に立つレベルであるということは、

「少し勉強(おさらい)しておけば、何かしらの役には立つんじゃね?」

と思ったので、ちょっと勉強してみることにしました。

そういうことなので、間違いとかそこそこあるかもしれませんが、適当に補って読んでいただけると幸いです。

最初の質問に対する回答

まず、最初の質問に対する回答から考えてみます。

FFTの入力とする観測データ(データはすべて実数値であるとします。)を得るために要した時間を$T$とします。また、観測データのサンプリングの間隔を$N$とすると、この観測データのサンプリング周波数$f_s$は
\[
f_s = \frac{N}{T}
\]
になります。

この観測データに対してFFTを実行した結果得られるデータは(複素共役になっているナイキスト周波数より後の部分は負の周波数のデータであり、対応する正の周波数のデータが存在しているので、ここでは除外して考えると)、$\displaystyle\frac{N}{2}+1$点の周波数領域のスペクトルデータになります(“+1″はナイキスト周波数成分に相当します)。また、データの示すスペクトルデータは周波数$0$(直流成分)から$\displaystyle\frac{f_s}{2}$(ナイキスト周波数成分)までのデータなので、周波数領域におけるスペクトルデータの間隔(周波数分解能)$f_r$は、
\[
f_r = \frac{f_s/2}{N/2} = \frac{f_s}{N} = \frac{1}{T}
\]
となります。

冷静に考えればそれほど難しい話ではないんですが、きれいさっぱり忘れてましたね。

~~(急に質問が来たからということにしておこう…)。~~

ちょっと気になったので、計算してみる。

ここでは、具体的なフーリエ変換の計算を行わないまま記事を書いてきましたが、いったん考えたり調べたりし始めると

「そういえば、正弦波のDFTってどうやって求めるんだっけ? ってか、なんで$\delta$関数になるんだっけ?」

とか、気になり始めたりして夜も眠れなくなることは以前ほどはなくなった(ベッドをシモンズのやつに交換したので、寝つきが格段に良くなりました。(｀・ω・´))のですが、それでも導出方法くらいは思い出しておきたかったので、他のサイトの記事を参考にしつつ、計算してみることにしました。

まず、余弦関数のDFTを計算します。

まず、周波数$f=\displaystyle\frac{f_0}{T}\left(0 < f_0 \le \displaystyle\frac{N}{2}, f_0 \in \mathbb{N} \right)$の余弦関数(のサンプリングデータ)
\[
x_t = \cos \left(\displaystyle\frac{2\pi j f_0 t}{N}\right)
\]
($j$は虚数単位)に対してDFTを実行してみます。上式の右辺に登場する$f_0$は実際の周波数と周波数分解能の比を表します。また、$t$は観測データの観測間隔を1とした相対的な経過時刻です。前節での議論の通り、観測データを得るために要した時間は周波数分解能の逆数になっているために打ち消し合うので、上記の右辺には直接現れることはありません。

(↑の議論は忘れがちになるところです。)

すると、変換の結果得られるスペクトルデータ$X_f$は以下の式で表されます。
\[
X_f = \sum^{N-1}_{t=0}\,x_t\,e^{-\frac{2\pi jft}{N}}
\]

$x_t = \displaystyle\frac{1}{2}\left(e^{\frac{2\pi j f_0 t}{N}}+e^{-\frac{2\pi j f_0 t}{N}}\right)$と変形してから上式に代入すると、

\[
X_f = \sum^{N-1}_{t=0}\,\displaystyle\frac{1}{2}\left\{e^{\frac{2\pi j t}{N}(f_0-f)}+e^{\frac{2\pi j t}{N}(-f_0-f)}\right\}
\]

上式の右辺第2項に$e^{-2\pi j t}$($t$が整数なので1になります。)をかけると、

\begin{eqnarray}
X_f &=& \sum^{N-1}_{t=0}\,\displaystyle\frac{1}{2}\left\{e^{\frac{2\pi j t}{N}(f_0-f)}+e^{-2\pi j t}\,e^{\frac{2\pi j t}{N}(-f_0-f)}\right\} \nonumber \\
{} &=& \sum^{N-1}_{t=0}\,\displaystyle\frac{1}{2}\left\{e^{\frac{2\pi j t}{N}(f_0-f)}+e^{-\frac{2\pi j t}{N}(N-f_0-f)}\right\} \nonumber \\
{} &=& \sum^{N-1}_{t=0}\,\displaystyle\frac{1}{2}e^{\frac{2\pi j t}{N}(f_0-f)}+\sum^{N-1}_{t=0}\,\displaystyle\frac{1}{2}e^{-\frac{2\pi j t}{N}(N-f_0-f)} \label{eq:fftfirsteq}
\end{eqnarray}

と式変形できます。

(\ref{eq:fftfirsteq})式右辺の各項はともに等比級数になっているので、まず第1項について計算してみます。

第1項は$f = f_0$か否かで場合分けが必要で、$f = f_0$の場合には冪の部分が0になりますので、$e^{\frac{2\pi j t}{N}(f_0-f)} = 1$になります。また、$f \neq f_0$の場合には、

\begin{eqnarray}
\sum^{N-1}_{t=0}\,e^{\frac{2\pi j t}{N}(f_0-f)} &=& \displaystyle\frac{1-e^{2\pi j(f_0-f)}}{1-e^{\frac{2\pi j t}{N}(f_0-f)}}
\end{eqnarray}

と計算できます。$f_0$及び$f$は整数なので、上式右辺の分子の値は0となります。よって、上式左辺の値も0となりますので、(\ref{eq:fftfirsteq})式右辺の第1項はクロネッカーのデルタ$\delta_{t,t^{\prime}}$を用いて、

\begin{eqnarray}
\sum^{N-1}_{t=0}\,e^{\frac{2\pi j t}{N}(f_0-f)} = N\delta_{f_0,f} \label{eq:deltaf0f}
\end{eqnarray}

となります。(\ref{eq:fftfirsteq})式右辺第2項についても$N-f_0$及び$f$について同様の議論ができますので、

\begin{eqnarray}
\sum^{N-1}_{t=0}\,e^{\frac{2\pi j t}{N}(N-f_0-f)} = N\delta_{N-f_0,f} \label{eq:deltaNf0f}
\end{eqnarray}

となります。(\ref{eq:deltaf0f}),(\ref{eq:deltaNf0f})式を(\ref{eq:fftfirsteq})式の右辺に代入すると、

\begin{eqnarray}
X_f &=& \displaystyle\frac{1}{2}(N\delta_{f_0,f} + N\delta_{N-f_0,f}) \\
{} &=& \displaystyle\frac{N}{2}(\delta_{f_0,f} + \delta_{N-f_0,f})
\end{eqnarray}

となります。正の周波数と負の周波数の成分が出てきますが、$N \to \infty$とすると、デルタ関数っぽくなりますね。

正弦関数のDFT等を計算します。

次に、
\[
x_t = \sin \left(\displaystyle\frac{2\pi j f_0 t}{N}\right)
\]

を計算してみます。これは、余弦関数の場合と同様に

\begin{eqnarray}
x_t &=& \sin \left(\displaystyle\frac{2\pi j f_0 t}{N}\right) \\
{} &=& \displaystyle\frac{1}{2j}\left(e^{\frac{2\pi j f_0 t}{N}}-e^{-\frac{2\pi j f_0 t}{N}}\right)
\end{eqnarray}

と表せることを利用して、

\begin{eqnarray}
X_f &=& \sum^{N-1}_{t=0}\,\displaystyle\frac{1}{2j}(\delta_{f_0,f} – \delta_{N-f_0,f}) \\
{} &=& \displaystyle\frac{N}{2j}(\delta_{f_0,f} – \delta_{N-f_0,f})
\end{eqnarray}

と計算できます。2個のデルタ関数っぽい関数の重ね合わせで表現されるのは余弦関数の場合と同様ですが、正の周波数成分は余弦関数の正の周波数成分を複素平面上で$- \displaystyle\frac{\pi}{2}$だけ、負の周波数成分は余弦関数の負の周波数成分を複素平面上で$\displaystyle\frac{\pi}{2}$だけ回転したものになっています。

もうちょっと一般に、位相(初期位相)$\phi$の正弦関数

\[
x_t = \sin \left(\displaystyle\frac{2\pi j f_0 t}{N}+\phi\right)
\]

のDFTは

\begin{eqnarray}
X_f&=&\displaystyle\frac{N}{2j}(\delta_{f_0,f}e^{j\phi} – \delta_{N-f_0,f}e^{-j\phi})
\end{eqnarray}

になります。

まとめ

ここまで超簡単ではありますが、フーリエ変換、特にDFTを中心におさらいしたことを書きました。

最近はfftwみたいな優秀なライブラリがあったりするので、自分でFFTのプログラムを作らなくてもデータを食わせれば計算結果を得ることはできます。ただ、計算結果の見方や計算結果に対する解析の方法くらいは頭に入っていないと、また次に急に質問が来た時に対応できないので、おさらいのタイミングとしては悪くはなかったのではないかと思います。