ディジタル信号処理｜周波数解析(FFT)とは

今回のテーマはディジタル信号処理です🐜
その中でも、周波数解析またはFFTの要点について整理したいと思います。¹

もちろん、継続的に加筆修正していきます！

周波数解析
フーリエ変換｜連続時間領域と連続周波数領域
離散時間フーリエ変換｜離散時間領域と連続周波数領域
離散フーリエ変換｜離散時間領域と離散周波数領域
高速フーリエ変換(FFT)
フーリエ変換のまとめ
実装時のポイント
付録｜複素数オイラーの公式
おわりに
参考文献

周波数解析

周波数解析（スペクトル解析）とは、信号にどんな周波数の周期信号が含まれているのかを解析することである。

信号にフーリエ変換を施して、時間領域から周波数領域に変換して解析する。

フーリエ変換｜連続時間領域と連続周波数領域

フーリエ変換は、信号を連続時間領域→連続周波数領域に変換する。 もちろん、フーリエ逆変換または逆フーリエ変換は、信号を連続周波数領域→連続時間領域に変換する。

フーリエ変換の定義式を以下に示す。

$\displaystyle F(\Omega) = \int^{\infty}_{-\infty} f(t) e^{-j \Omega t} dt$

フーリエ逆変換の定義式を以下に示す。

$\displaystyle f(t) = \frac{1}{2 \pi} \int^{\infty}_{-\infty} F(\Omega) e^{ j \Omega t} d \Omega$

ただし、 $T$ を周期、 $f$ を周波数としたとき、 $\Omega = 2 \pi f = \frac{2 \pi}{T}$ が成り立つ。この $\Omega$ は（非正規化）角周波数を表す連続変数であり、 $F(\Omega$ )は $f(t)$ に含まれる角周波数 $\Omega$ の振動成分の量を示す。ちなみに、 $|F(\Omega)|$ は $f(t)$ の振幅スペクトル、 $\angle F(\Omega)$ は $f(t)$ の位相スペクトル、 $|F(\Omega)|^2$ は $f(t)$ のパワースペクトルと呼ぶ。

これを線形・時不変システムで図示すると以下の通り。

なお、畳み込み式は次式で定義される。

$\displaystyle y(t) = \int^{\infty}_{-\infty} x(\tau) h(t-\tau) d\tau = x(t) * h(t)$

しかし、計算機（ディジタルな領域）で解析するには、離散化が必要であることがわかる。

離散時間フーリエ変換｜離散時間領域と連続周波数領域

離散時間フーリエ変換(Dicrete-Time Fourier Transform; DTFT)は、信号を離散時間領域→連続周波数領域に変換する。 もちろん、離散時間フーリエ逆変換は、信号を連続周波数領域→離散時間領域に変換する。

離散時間フーリエ変換の定義式を以下に示す。

$\displaystyle F(\omega) = \sum^{\infty}_{n=-\infty} f[n] e^{-j \omega n}$

離散時間フーリエ逆変換の定義式を以下に示す。

$\displaystyle f[t] = \frac{1}{2 \pi} \int^{\pi}_{-\pi} F(\omega) e^{ j \omega n} d \omega$

この $\omega$ は正規化角周波数を表す連続変数であり、 $F(\omega)$ は $f[n]$ に含まれる正規化角周波数 $\omega$ の振動成分の量を表す。また、 $F(\omega)$ は周期 $2\pi$ の周期関数で、通常 $-\pi \lt \omega \lt \pi$ の範囲で考える。

なお、サンプリング周期を $T_s$ としたとき、非正規化角周波数 $\Omega$ と正規化角周波数 $\omega$ の間に以下の式が成り立つ。

$\displaystyle \omega = \Omega T_s$

これを線形・時不変システムで図示すると以下の通り。

しかし、計算機で解析するには、無限和の解消と周波数領域の離散化が必要であることがわかる。

離散フーリエ変換｜離散時間領域と離散周波数領域

離散フーリエ変換(Discrete Fourier Transform; DFT)は、信号を離散時間領域→離散周波数領域に変換する。 もちろん、離散フーリエ逆変換は、信号を離散周波数領域→離散時間領域に変換する。

離散フーリエ変換の定義式を以下に示す。

$\displaystyle F[k] = \sum^{N-1}_{n=0} f[n] e^{-j \frac{2\pi}{N} kn}$

離散フーリエ逆変換の定義式を以下に示す。

$\displaystyle f[n] = \frac{1}{N} \sum^{N-1}_{k=0} F[k] e^{j \frac{2\pi}{N} kn}$

$f[n]$ は離散時間 $n$ （整数）で定義された周期 $N$ の周期関数であり、 $F[k]$ も整数 $k$ で定義される周期 $N$ の周期関数である。 $k$ は周波数のインデックスを表す整数で、 $F[k]$ は $f[n]$ に含まれる正規化角周波数 $\frac{2\pi k}{N}$ [rad]の振動成分の量を表す。

これを線形・時不変システムで図示すると以下の通り。

なお、畳み込み式は次式で定義される。

$\displaystyle y[n] = \sum^{\infty}_{k=-\infty} x[k] h[n-k] = x[t] * h[t]$

高速フーリエ変換(FFT)

結局、周波数解析で求めたいのは $F(\Omega)$ または $|F(\Omega)|$ であるが、計算機で解析するために様々な近似を経て、 $F[k]$ を計算して解析を行う。

しかし、DFTでも計算時間が $O(N^{2})$ であり、現実的ではない。

そこで一般には、この計算時間が $O(N \log_2 N)$ になる高速フーリエ変換(Fast Fourier Transform; FFT)を利用する。

このFFTでは、そのアルゴリズムや計算機の性質から $N=1024$ など2の乗数の点数で計算を行うと高速にできる。

アルゴリズムの詳細はここでは議論しない。

フーリエ変換のまとめ

ここまで3種のフーリエ変換について概要を述べたが、これらにフーリエ級数展開を加えて、離散性と周期性に注目して以下に整理する。

なお、フーリエ級数の定義式を以下に示す。ただし、ここでは議論しない。

$\displaystyle F_k = \frac{1}{T_0} \int^{T_0 /2}_{-T_0 /2} f(t) e^{-j \Omega_0 kt} dt$

$\displaystyle f(t) = \sum^{\infty}_{k=-\infty} F_k e^{j \Omega_0 kt}$

実装時のポイント

Pythonの信号処理ライブラリで有名なSciPyを使う場合、主に以下の2つの関数を使います。

scipy.fftpack.fft
中身はDFTのFFTです。戻り値は $F[k]$ ですので、複素数であることに注意しましょう。
また、一般に周波数解析でy軸にとるのは $|F[k]|$ なので、適宜パワースペクトルに変換します。
scipy.fftpack.fftfreq
scipy.fftpack.fftの戻り値に対応する、離散化された周波数を計算する関数です。
必要に応じて正規化各周波数→非正規化角周波数に変換します。