円周率.jp (http://xn--w6q13e505b.jp/method/fft/rft.html)

実数FFT

DFT は複素数の計算として提唱され,高速化に関する研究がなされているが, 実際には入力値として実数を用いることが一般的である. 入力値が実数であるということは入力に関する虚部に関するデータが必要無い. そのため,理論的には変数の領域が半分で済む上に演算量も半分程度に抑えられるはずであり,それは実現可能な事実である. 本ページではそのような実数 FFT を計算する上で必要な事を記述する.

ここでは多数桁に応用することを前提にしたものを含めた実数 FFT について記述していく.

共役複素数を利用した実数FFT

入力値が実数であるということは $a_j=\overline{a_j}$ である.このことから

\[ A_k = \sum_{j=0}^{n-1} a_j \omega^{jk} = \overline{ \sum_{j=0}^{n-1} \overline{a_j \omega{jk}} } = \overline{ \sum_{j=0}^{n-1} a_j \omega^{-jk}} = \overline{ A_{n-k} } \]

ということが分かる. $k \in (0, n/2)$ について $A_k$ の実部を ${\tt [k]}$ の位置に, 虚部を ${\tt [n-k]}$ の位置に格納することで,実質的に $n-2$ の複素数を同数の実数領域に収めることができる. 一方,$k=0$,$n/2$ の場合には $A_k$ が実数となり虚部を保持する必要がないので, 全体としてちょうど $n$ の配列で保持できる.

\[ a_0 = {\tt a[0]}, a_{n/2} = {\tt a[n/2]}, a_k = {\tt a[k]} + i{\tt a[n-k]}, a_{n-k} = {\tt a[k]} - i{\tt a[n-k]} \quad ({\rm for}\ 0\lt k\lt n/2) \]

複素数FFTを利用した実数FFT

\[ x_j = a_{2j} + ia_{2j+1} \quad ({\rm for}\ 0\leq j \lt n/2) \tag{1} \]

という形の複素数列 $\{x_j\}$ の DFT 結果 $\{X_k\} = \left\{ \sum_{j=0}^{n/2-1} x_j\omega^{2jk} \right\}$ を利用して

\[ \begin{eqnarray} A_k &=& X_k &-& \dfrac12(X_k - \overline{X_{n/2-k}}) (1 + i\omega^k)\\ \overline{A_{n/2-k}}&=& \overline{X_{n/2-k}} &+& \dfrac12(X_k - \overline{X_{n/2-k}}) (1 + i\omega^k)\\ \end{eqnarray} \] \[ A_0 = \Re (X_0) + \Im (X_0), A_{n/2} = \Re(X_0) - Im(X_0) \]

で計算できる. 上の 2 式は括弧内が共通部分しているので, 並行して計算を行うと効率よく演算でき, $k=n/4$ のときは $1+i\omega^{n/4}=2$ なので $X_{n/4}=\overline{A_{n/4}}$, $X_{3n/4} = \overline{A_{3n/4}}$ となる.

同様に変形することで,逆演算は

\[ \begin{eqnarray} X_k &=& A_k &-& \dfrac12 (A_k - \overline{A_{n/2-k}}) (1+i\omega^{-k})\\ \overline{X_{n/2-k}} &=& \overline{A_{n/2-k}} &+& \dfrac12 (A_k - \overline{A_{n/2-k}}) (1+i\omega^{-k}) \end{eqnarray} \] \[ X_0 = \dfrac12 (A_0 + A_{n/2}) + \dfrac{i}{2} (A_0 - A_{n/2}) \]

で計算できる. この変換を行った後に要素数 $n/2$ の複素数 IFFT を行うと式 (1) の形で実数 IFFT の結果が出る.

計算結果については上記 「共役複素数を利用した実数FFT」 に記載した性質を利用すれば,$0\leq k \leq n/4$ の計算と $n$ の記憶領域ですむ上に,1 つの複素数を連続した領域に保持するのが容易である.

導出方法

この式ができるまでの導出を簡単ながら行ってみる.

まず $X_k$ の定義を展開する.

\[ \begin{eqnarray} X_k &=& \sum_{j=0}^{n/2-1} x_j \omega^{2jk} &=& \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} + i\sum_{j=0}^{n/2-1} a_{2j+1}\omega^{2jk}\\ X_{n/2-k} &=& \sum_{j=0}^{n/2-1} \overline{x_j \omega^{2jk}} &=& \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} - i\sum_{j=0}^{n/2-1} a_{2j+1}\omega^{2jk} \end{eqnarray} \] \[ \therefore \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} = \frac12 (X_k + \overline{X_{n/2-k}}), \quad \sum_{j=0}^{n/2-1} a_{2j+1}\omega^{2jk} = -\frac{i}{2} (X_k -\overline{X_{n/2-k}}) \]

これを利用して元の DFT の定義式を若干変形すると

\[ \begin{eqnarray} A_k &=& \sum_{j=0}^{n/2-1} a_{jk}\omega^{2jk} + \omega^k \sum_{j=0}^{n/2-1} a_{2j+1} \omega^{2jk} &=& X_k - \dfrac{1}{2} (X_k - \overline{X_{n/2-k}})(1+i\omega^k)\\ \overline{A_{n/2-k}} &=& \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} - \omega^k \sum_{j=0}^{n/2-1} a_{2j+1} \omega^{2jk} &=& \overline{X_{n/2-k}} + \dfrac{1}{2} (X_k - \overline{X_{n/2-k}}) (1 + i\omega^k) \end{eqnarray} \]

から計算できる. なお,$A_0$,$A_{n/2}$ の導出については,定義式から自明である.