panda's tech note

デジタル信号処理:高速フーリエ変換

離散フーリエ変換・逆離散フーリエ変換 で説明した通り,離散フーリエ変換・逆離散離散フーリエ変換を適用することで,デジタル信号(離散信号)に対して時間領域と周波数領域の変換を行うことができる。以下に離散フーリエ変換および逆離散フーリエ変換の定義を再掲する。

\begin{align} X_n = \sum_{k=0}^{N-1} x_k e^{-\frac{2\pi}{N}ikn} \label{eq:dft} \\ x_k = \frac{1}{N} \sum_{n=0}^{N-1} X_n e^{\frac{2\pi}{N}ikn} \label{eq:idft} \end{align}

つまり,標本数 \( N \) の離散フーリエ変換の計算量は,単純に式 \ref{eq:dft} の計算を \( 0 \le n < N \) に対して行うと \( O(N^2) \) となる。

高速フーリエ変換とCooley-Tukey アルゴリズム

上述したように,ナイーブな手法で計算すると離散フーリエ変換の計算量は \(O(N)\) となるが,分割統治法を用いることで \( O(N \log N) \) と高速化することがでる。これは高速フーリエ変換 (FFT: Fast Fourier Transform) と呼ばれ,様々なデジタル信号処理で用いられている。ここでは高速フーリエ変換のアルゴリズムとして最も有名な Cooley–Tukey アルゴリズムについて説明する。なお,高速フーリエ変換を実装しているライブラリとしては FFTW が有名である。

Cooley-Tukey アルゴリズムでは,\( N \) を2の冪乗として分割統治法を用いる。式 \ref{eq:dft} を \( n = \{ 0, \cdots, N-1 \} \) に対する \( X_n \) をベクトルとして表し, \(\boldsymbol{x} = [x_0, \cdots, x_{N-1}]^T, \boldsymbol{X} = [X_0, \cdots, X_{N-1}]^T \) とし,\( w_N = e^{-\frac{2\pi}{N}i} \) とおくと,

\begin{align*} \boldsymbol{X} &= \begin{bmatrix} w_N^{0 \cdot 0} & \cdots & w_N^{0 \cdot (N-1)} \\ \vdots & \ddots & \\ w_N^{(N-1) \cdot 0} & & w_N^{(N-1) \cdot (N-1)} \end{bmatrix} \boldsymbol{x} \end{align*}

と書ける。ここで,

$$ W_N = \begin{bmatrix} w_N^{0 \cdot 0} & \cdots & w_N^{0 \cdot (N-1)} \\ \vdots & \ddots & \\ w_N^{(N-1) \cdot 0} & & w_N^{(N-1) \cdot (N-1)} \end{bmatrix} $$

(\(m\)行\(n\)列の要素が \( w^{mn}\) の \(N \times N\) 行列)とすると,

$$ \boldsymbol{X} = W_N \boldsymbol{x} $$

となる。\(W_N\) の偶数列と奇数列を分け,\( \boldsymbol{x_e} = [x_0, x_2, \cdots, x_{N-2}]^T, \boldsymbol{x_o} = [x_1, x_3, \cdots, x_{N-1}]^T \) とすると,

\begin{align} \boldsymbol{X} &= \begin{bmatrix} w_N^{0 \cdot 0} & \cdots & w_N^{0 \cdot (N-2)} & w_N^{0 \cdot 1} & \cdots & w_N^{0 \cdot (N-1)} \\ \vdots & \ddots & & \vdots & \ddots & \\ w_N^{(N-1) \cdot 0} & \cdots & w_N^{(N-1) \cdot (N-2)} & w_N^{(N-1) \cdot 1} & \cdots & w_N^{(N-1) \cdot (N-1)} \end{bmatrix} \begin{bmatrix} \boldsymbol{x_e} \\ \boldsymbol{x_o} \end{bmatrix} \label{eq:fft1} \end{align}

\( w_{N/2} = w_{N}^2 \) より,

\begin{align*} \begin{bmatrix} w_N^{0 \cdot 0} & \cdots & w_N^{0 \cdot (N-2)} \\ \vdots & \ddots & \\ w_N^{(N-1) \cdot 0} & \cdots & w_N^{(N-1) \cdot (N-2)} \end{bmatrix} &= \begin{bmatrix} w_N^{0 \cdot 0} & \cdots & w_N^{0 \cdot (N-2)} \\ \vdots & \ddots & \\ w_N^{(N/2-1) \cdot 0} & \cdots & w_N^{(N/2-1) \cdot (N-2)} \\ w_N^{(N/2) \cdot 0} & \cdots & w_N^{(N/2) \cdot (N-2)} \\ \vdots & \ddots & \\ w_N^{(N-1) \cdot 0} & \cdots & w_N^{(N-1) \cdot (N-2)} \end{bmatrix} \\ &= \begin{bmatrix} w_{N/2}^{(0 \cdot 0)/2} & \cdots & w_{N/2}^{(0 \cdot (N-2))/2} \\ \vdots & \ddots & \\ w_{N/2}^{((N/2-1) \cdot 0)/2} & \cdots & w_{N/2}^{((N/2-1) \cdot (N-2))} \\ w_{N/2}^{((N/2) \cdot 0)/2} & \cdots & w_{N/2}^{((N/2) \cdot (N-2))/2} \\ \vdots & \ddots & \\ w_{N/2}^{((N-1) \cdot 0)/2} & \cdots & w_{N/2}^{((N-1) \cdot (N-2))/2} \end{bmatrix} \\ &= \begin{bmatrix} W_{N/2} \\ W_{N/2} \end{bmatrix} \end{align*}

となる。また,

\begin{align*} \begin{bmatrix} w_N^{0 \cdot 1} & \cdots & w_N^{0 \cdot (N-1)} \\ \vdots & \ddots & \\ w_N^{(N-1) \cdot 1} & \cdots & w_N^{(N-1) \cdot (N-1)} \end{bmatrix} &= \begin{bmatrix} w_{N}^0 w_N^{0 \cdot 0} & \cdots & w_{N}^{0} w_N^{0 \cdot (N-2)} \\ \vdots & \ddots & \\ w_{N}^{N/2-1} w_N^{(N/2-1) \cdot 0} & \cdots & w_{N}^{N/2-1} w_N^{(N/2-1) \cdot (N-2)} \\ w_{N}^{N/2} w_N^{(N/2) \cdot 0} & \cdots & w_{N}^{N/2} w_N^{(N/2) \cdot (N-2)} \\ \vdots & \ddots & \\ w_{N}^{N-1} w_N^{(N-1) \cdot 0} & \cdots & w_{N}^{N-1} w_N^{(N-1) \cdot (N-2)} \end{bmatrix} \\ &= \begin{bmatrix} w_{N}^0 & & \boldsymbol{0} \\ & \ddots & \\ \boldsymbol{0} & & w_{N}^{N/2-1}\\ w_{N}^{N/2} & & \boldsymbol{0} \\ & \ddots & \\ \boldsymbol{0} & & w_{N}^{N-1}\\ \end{bmatrix} \begin{bmatrix} W_{N/2} \\ W_{N/2} \end{bmatrix} \end{align*}

となる。ここで,

\begin{align*} C_N &= \begin{bmatrix} w_{N}^0 & & \boldsymbol{0} \\ & \ddots & \\ \boldsymbol{0} & & w_{N}^{N/2-1} \\ \end{bmatrix} \quad (\in \mathbb{C}^{N/2 \times N/2}) \\ D_N &= \begin{bmatrix} w_{N}^{N/2} & & \boldsymbol{0} \\ & \ddots & \\ \boldsymbol{0} & & w_{N}^{N-1} \end{bmatrix} \quad (\in \mathbb{C}^{N/2 \times N/2}) \end{align*}

とおくと,式 \ref{eq:fft1} は,\(\boldsymbol{1}\) を単位行列,\(\boldsymbol{0}\) をゼロ行列とすると,

\begin{align} \boldsymbol{X} &= \begin{bmatrix} W_{N/2} & C_N W_{N/2} \\ W_{N/2} & D_N W_{N/2} \end{bmatrix} \begin{bmatrix} \boldsymbol{x_e} \\ \boldsymbol{x_o} \end{bmatrix} \notag \\ &= \begin{bmatrix} \boldsymbol{1} & C_{N} \\ \boldsymbol{1} & D_{N} \end{bmatrix} \begin{bmatrix} W_{N/2} & \boldsymbol{0} \\ \boldsymbol{0} & W_{N/2} \end{bmatrix} \begin{bmatrix} \boldsymbol{x_e} \\ \boldsymbol{x_o} \end{bmatrix} \notag \end{align}

と表せる。ここで,

$$ \begin{bmatrix} W_{N/2} & \boldsymbol{0} \\ \boldsymbol{0} & W_{N/2} \end{bmatrix} \begin{bmatrix} \boldsymbol{x_e} \\ \boldsymbol{x_o} \end{bmatrix} = \begin{bmatrix} W_{N/2} \boldsymbol{x_e} \\ W_{N/2} \boldsymbol{x_o} \end{bmatrix} $$

となる。\( W_{N/2} \boldsymbol{x_e}, W_{N/2} \boldsymbol{x_o}\) はそれぞれ,\( \frac{N}{2} \) のサイズの離散フーリエ変換である。また,\( C_N, D_N \) は対角行列なので,

$$ \boldsymbol{X} = \begin{bmatrix} \boldsymbol{1} & C_{N} \\ \boldsymbol{1} & D_{N} \end{bmatrix} \begin{bmatrix} W_{N/2} \boldsymbol{x_e} \\ W_{N/2} \boldsymbol{x_o} \end{bmatrix} = \begin{bmatrix} (W_{N/2} \boldsymbol{x_e}) + C_N (W_{N/2} \boldsymbol{x_o}) \\ (W_{N/2} \boldsymbol{x_e}) + D_N (W_{N/2} \boldsymbol{x_o}) \end{bmatrix} $$

の演算は \( W_{N/2} \boldsymbol{x_e}, W_{N/2} \boldsymbol{x_o}\) が与えられたとき \( O(N) \) の計算量で行える。これを再帰的に繰り返すと \(N=1\) のサイズの離散フーリエ変換に帰着できる(ただし,多くの実装では十分小さいサイズの離散フーリエ変換は \(O(N^2)\) の離散フーリエ変換を行い,\(N=1\) まで分割しない)。このように分割統治法により,半分のサイズの高速フーリエ変換を再帰的に計算することで,\( O(N \log N) \) の高速フーリエ変換が実現できる。