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

基底

FFT の計算時間オーダーが $O(n \log n)$ になることについては DFT をサイズが半分の DFT 2つに分割する形で説明した. この様に FFT を 2 つに分けていく方法を 2 基底 FFT という. しかし FFT を行う上で,基底は 2 以外でも自由に取ることができる.

また,基底の異なる FFT をコンピュータプログラムとして動作させる場合, 演算量や演算量に対するメモリ Load/Store の割り合いなどが原因となり速度に違いが出てくる.

ここでは2基底3基底4基底の FFT を具体例とし, FFT の概念的な説明と, それらをプログラムで実装する際の演算の省略について述べる.

2基底FFT

まずは比較対象となる 2 基底 FFT の説明を行う. 時間間引きで見た FFT の形をアルゴリズムとして表すと以下のようになる. 簡単のため $n=2^p$ として計算する.

図 11: 2基底FFTのイメージ図

後の説明のため,データの配置を2次元で表現すると, 高さ $H$,幅 $W$ (各変数の値は下記擬似コード参照)の四角形を 2 個並べることで,各ステップでの計算対象位置・ 格納位置が理解しやすくなる. (図1)

なお,1 次元配列へは $A_L[k,j]=a_{2Wk+j}$ で変換できる.

  • for $L$ = 0,1,2,...,$p-1$ :
    • $H=2^L$, $W=2^{p-L-1}$
    • for $k$ = 0,1,2,...,$H-1$ :
      • for $j$ = 0,1,2,...,$W-1$ :
        • $A_{L+1}[k,j] = A_L[k,j] + \omega^{kW}A_L[k,j+W]$
        • $A_{L+1}[k+H,j] = A_L[k,j] - \omega^{kW}A_L[k,j+W]$
      • end of for
    • end of for
  • end of for

次にこの演算は(2 基底の FFT)の演算量を数えると,最内ループ内で 1 回の複素数乗算と 2 回の複素数加減算(以下,単純に加算で代表する) を行っている. さらにそれが $WH(=n/2)$ 回繰り返され,$L$ のループで $p(=\log_2n)$ 回繰り返されるので,全体で $(1/2)n\log_2 n$ 回の複素数乗算と $n\log_2n$ 回の複素数加算が,つまり $2n\log_2n$ 回の実数乗算と $3n\log_2n$ 回の実数加算が行われる.

3基底FFT

図 12: 3基底FFT

$O(n\log n)$ で FFT を行うためには,データ数 $n$ が 2 のべき乗である必要はない. 基本的な考え方を同じくすれば 3 基底の FFT を作成することが可能である. ここでは $n=3^p$ として 3 基底の FFT アルゴリズムを示す.

また,表現の単純化のため $\Omega = \exp(-2\pi i/3) = \dfrac{-1-\sqrt{3}i}{2}$ とする.

  • for $L$ = 0,1,2,...,$p-1$ :
    • $H=3^L$, $W=3^{p-L-1}$
    • for $k$ = 0,1,2,...,$H-1$ :
      • for $j$ = 0,1,2,...,$W-1$ :
        • $A_{L+1}[k,j]=A_L[k,j] + \omega^{kW}A_L[k,j+W] + \omega^{2kW}A_L[k,j+2W]$
        • $A_{L+1}[k+H,j]=A_L[k,j] + \Omega\omega^{kW}A_L[k,j+W] + \Omega^2\omega^{2kW}A_L[k,j+2W]$
        • $A_{L+1}[k+2H,j]=A_L[k,j] + \Omega^2\omega^{kW}A_L[k,j+W] + \Omega\omega^{2kW}A_L[k,j+2W]$
      • end of for
    • end of for
  • end of for

ちなみに,ここで $\Omega$ と表した部分が 2 基底 FFT では $\exp(\pi i)=-1$ となるため 2 基底 FFT の片方は $-$ になっており, 乗算を行う必要がないため演算量が大きく削減されている.

具体的な演算数は複素数乗算,複素数加算がそれぞれ $8n\log_3n$ 回ずつであるが, $\Omega^2=\overline\Omega$ であることを利用することで, 実数乗算が $\dfrac{20}{3}n\log_3 n \simeq 4.21n\log_2n$ 回,実数加算が $8n\log_3n \simeq 5.05n\log_2n$ 回となるため ほぼ同数のデータに対する 2 基底 FFT とくらべて約 2 倍の演算量となる.

4基底FFT

4 基底 FFT についても 3 基底 FFT と同様に構成することが可能であるが, 4 = 2×2 であることを利用するとさらに計算量を減らすことができる. つまり,2 基底の 2 段ループアンローリングという形で考え, $L$ のループについて 0,2,4,… という動きができるようにする.

図 13: 4基底FFTのイメージ図

基本的には図 3 に示した概念図のうち,$L+1$ の段階でメモリへの読み書きを減らし, また,回転子の乗算を纏めて行い演算回数も減らすことで速度を速めている. 以降の説明では $n=2^p$ (ただし $p$ は偶数)とする。

  • for $L=0,2,4,\cdots p-2$ :
    • $H=2^L$,$W=2^{p-L-2}$
    • for $k=0,1,2,\cdots,H-1$ :
      • for $j=0,1,2,\cdots,W-1$ :
        • この部分
      • end of for
    • end of for
  • end of for

だけを取り出して表示する. 2 基底 FFT の考え方を単純に適用すると

$L\rightarrow L+1$ の段階

  • $A_{L+1}[k ,j ] = A_L[k,j ] + \omega^{2kW} A_L[k,j+2W]$
  • $A_{L+1}[k+H,j ] = A_L[k,j ] - \omega^{2kW} A_L[k,j+2W]$
  • $A_{L+1}[k ,j+W] = A_L[k,j+W] + \omega^{2kW} A_L[k,j+3W]$
  • $A_{L+1}[k+H,j+W] = A_L[k,j+W] - \omega^{2kW} A_L[k,j+3W]$

$L+1\rightarrow L+2$ の段階

  • $A_{L+2}[k ,j] = A_{L+1}[k ,j] + \omega^{kW} A_{L+1}[k ,j+W]$
  • $A_{L+2}[k+ H,j] = A_{L+1}[k+H,j] - i\omega^{kW} A_{L+1}[k+H,j+W]$
  • $A_{L+2}[k+2H,j] = A_{L+1}[k ,j] - \omega^{kW} A_{L+1}[k ,j+W]$
  • $A_{L+2}[k+3H,j] = A_{L+1}[k+H,j] + i\omega^{kW} A_{L+1}[k+H,j+W]$

となるが,後半で行うはずの $\omega^{kW}$ の乗算を前半へ持って行くことで

$L\rightarrow L+1$ の段階

  • $a = A_L[k,j ] + \omega^{2kW} A_L[k,j+2W]$
  • $b = A_L[k,j ] - \omega^{2kW} A_L[k,j+2W]$
  • $c = \omega^{kW}A_L[k,j+W] + \omega^{3kW} A_L[k,j+3W]$
  • $d = \omega^{kW}A_L[k,j+W] - \omega^{3kW} A_L[k,j+3W]$

$L+1\rightarrow L+2$ の段階

  • $A_{L+2}[k ,j] = a + c$
  • $A_{L+2}[k+ H,j] = b - id$
  • $A_{L+2}[k+2H,j] = a - c$
  • $A_{L+2}[k+3H,j] = b + id$

とシンプルになる. 基本的に $P$ 基底 FFT に対する $P^m$ 基底 FFT も同様の変形を行うことで実現できる. 4 基底 FFT での $\omega^{3kW}$ をかける部分の様に $P$ 基底では $m$ 回に分けて回転子をかける (4 基底 では $\omega^{kW}$,$\omega^{2kW}$ の 2 回)必要があるものを $P^m$ 基底では 1 回の乗算に抑えることで演算量を下げることになる.

また,4 基底 FFT では $\Omega=-i$ ($\Omega$ については 3 基底 FFT 参照)であるために, 本来なら演算すべき($i$ をかける)はずが, 計算対象(実部と虚部)を取り換えるだけですむので実数乗算・実数加算の回数が減り, より効率的になる. この点から FORTRAN などの複素数型を扱える言語でも,複素数値を実数値×2 個として扱う方が計算効率がいい.

具体的な数で表すと,単純に 2 基底 FFT を展開しただけのものでは $2n\log_2n$2 の実数乗算と$3n\log_2n$ の実数加算が必要だったのに対し, 纏め直した形では $\dfrac{3}{2}n\log_2n$ の実数乗算と $\dfrac{11}{4}n\log_2n$ の実数加算で済むため, 実数演算量が 85% に削減される.

同様の展開をすることで 8 基底 FFT,16 基底 FFT …などを作ることができ, 8 基底 FFT では $1/\sqrt{2}$ を共有することで演算量を減らすことができる. しかし,16 以上の基底では回転子の乗算を纏めることで削減される演算数が少ない一方で, 作成すべきプログラムが非常に煩雑になるためコストパフォーマンスは低い.

$r$ 基底FFT

3 基底 FFT と同様にすることで,一般の基底 $r$ について $r$ 基底 FFT を構成することが可能である. より一般的に,データ数 $n$ が

\[ n = \prod_{k=0}^{m-1} r_k = r_0 \cdot r_1 \cdot \cdots r_{m-1} \]

で表される時の FFT が構成できる. $r_k$ は 2 以上の自然数であればどのような順序で並んでいても構わない. このとき,$r_k$ の順序に基づいて以下のようなアルゴリズムが構成される. なお,第 $L$ 回目のループで $A_L$ は $H_L\times r_LW_L$ の行列である.

  • $H_0=1$,$W_0=n$
  • for $L=0,1,2,\cdots,m-1$ :
    • $W_{L+1} = W_L \div r_L$
    • for $k=0,1,2,\cdots,H_L-1$ :
      • for $j=0,1,2,\cdots,W_L-1$ :
        • for $t=0,1,2,\cdots,r_L-1$ :
          • $A_{L+1}[k+tH_L,j] = {\displaystyle \sum_{s=0}^{r_L-1}} \Omega^{st} \omega^{skW_{L+1}} A_L[k,j+sW_{L+1}]$
        • end of for
      • end of for
    • end of for
    • $H_{L+1} = r_LH_L$
  • end of for

ここでは 2, 3, 4 基底と同じように書いたが,$k$,$j$,$t$,$s$($\sum$ の部分) のループについては互いに依存関係が無いので任意に入れ替えても問題ない. また 3 基底 FFT でも述べたように,$\Omega^{(r-s)t} = \overline{\Omega^{st}}$ が成り立つことを利用することで演算量が下がる.

また、上記の 4 基底 FFT のように いくつかの基底については各基数に依存した最適化を行うことで さらに演算数を削減することができる. 例えば 4 基底では $\Omega=i$ との乗算を実質的に省くことができた. 他の例としては 5 基底 FFT においては $\cos(2\pi/5)=\dfrac{\sqrt{5}-1}{4}$, $\cos(2\pi/5)=-\dfrac{\sqrt{5}+1}{4}$ であることが, 8 基底 FFT (細かくは 8 を因数に含む基底の FFT)では 4 基底 FFT と同じ手法の他, $\cos(2\pi/8)=\sin(2\pi/8)=\dfrac{1}{\sqrt{2}}$ であることが利用できる.

演算量

ここまでにいくつかの具体的な基底について, 各基底の特徴を利用した演算効率化を紹介したが, ここではそのような効率化を行わない前提で一般的な FFT の演算量を求める. まずデータ数 $n$ が

\[ n = \prod_{k=1}^{m} p_k^{e_k} \]

であるとき,1 つ 1 つの基底に対する演算量は表 1 のようになる. (具体的手段については $r$ 基底 FFT を参照.)

表1.p 基底 FFT の演算回数
複素数演算回数 実数演算回数
乗算回数 $(p_k-1) ne_k$ $4(p_k-1) ne_k$
加算回数 $(p_k-1)ne_k$ $4(p_k-1)ne_k$

そのため,全体としての演算回数は

\[ {\rm (実数乗算回数)} = {\rm (実数加算回数)} = 4n \sum_{k=1}^{m} (p_k-1) e_k \]

となる.

また,最近のプロセッサは演算速度に対するメモリの Load/Store の時間的コストが高い. そのため,キャッシュメモリに入りきらないサイズのデータが演算対象になる場合には 「単位演算量当たりの Load/Store 量」が重要になる.

$r$ 基底 FFT について考えると,データ量が $O(r)$,演算回数が $O(r^2)$ であるため,基本的には基底 $r$ が大きい方が有利である. しかし,本質的にはレジスタの数や, メモリ上の値を直接演算に用いる命令があるか, といった事柄が関わってくるので一概に大きい基底の方が良いとは限らない.