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

実数FFT

歴史的に見れば DFT は複素数の計算として提唱され高速化に関する研究がなされているが、 実際には入力値として実数を用いることが一般的である。 入力値が実数であるということは虚部が 0 であるという前提が成り立つので それを利用したメモリ削減や高速化も提唱されている。 本ページではそのような実数 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} } \]

ということが分かる。 これを利用して $0 \lt k \lt n/2$ について $A_k$ の実部を ${\tt [k]}$ の位置に、虚部を ${\tt [n-k]}$ の位置に格納すれば実質的に $n-2$ 個の複素数を同数の実数領域に収めることができる。 $k=0, n/2$ については $A_0=\sum a_j$、$A_{n/2}=\sum (-1)^ja_j$ と実数になり虚部を保持する必要がないので、 全体としてちょうど $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) \]

欠点としては共役の関係にある値の位置が逆順に並んでいるので インデックス管理が大変になることと、例外になる $k=0, n/2$ の処理を分けて作る必要があることが挙げられる。

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

実数 DFT で変換したいデータ $\{a_j\}$ の中身を 2 個ずつつなげて \[ 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) \tag{a}\\ \overline{A_{n/2-k}}&=& \overline{X_{n/2-k}} &+& \dfrac12(X_k - \overline{X_{n/2-k}}) (1 + i\omega^k) \tag{b}\\ \end{eqnarray} \] \[ A_0 = \Re (X_0) + \Im (X_0),\ A_{n/2} = \Re(X_0) - Im(X_0),\ A_{n/4} = \overline{X_{n/4}} \]

で計算できる。 特に (a) (b) の 2 式は括弧内が共通部分しているので並行して計算を行うと効率よく演算できる。 なお $n/2\lt k\lt n$ の $A_k$ については複素共役にあるデータ から復元できるので特に計算しない。

同様に変形することで逆演算は \[ \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}),\ X_{n/4} = \overline{A_{n/4}} \]

で計算できる。 この変換を行った後に $\{X_k\}$ の IFFT を行うと式 (1) の形で実数 IFFT の結果が出る。

導出方法

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

まず $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}$ の導出については、定義式から自明である。