部分再帰的FFT
実際、再帰的アルゴリズムに基づいてFFTを書くと、とてもシンプルなプログラムになります。アルゴリズムの基本を覚えている人ならば、「ははーん、クイックソートね」と気付いたかもしれません。クイック・ソート同様FFTの計算量もO(n)=N*log(N)になります*1。クイックソートを開発したC.A.R.Hoarは最初FORTRANで苦労して書いたそうです。で、新言語Algolのセミナーがあったので同僚と聞きに行った後、クイックソートの再帰版をAlgolで書き、「こりゃいい！」と感動したそうです。
閑話休題。
FFTの再帰版はシンプルで見通しが良いのですが、いかんせん効率が悪い*2ため、普段は使われません。一般には3重ループを使った「フラット」なアルゴリズムが使われます。
このフラットなアルゴリズムはN点の入力データをlog(N)回スキャンします。スキャンは連続的なものからステップ的なものまで変化しながら行われます。問題は、Nが大きい場合です。Nが小さい時はいいのですが、数万を超えると最新のx86と言えどもL1キャッシュに入らなくなります。スキャンは局所性が少なく、1回のスキャンにおいて、あるデータが読み書きされる回数はそれぞれ1回です。つまり、キャッシュの再利用はまったくなされません。スキャンがアドレスに対して連続的ならばまだ良いのですが、ステップ的になるとバースト・フィルが悪い方向に働き、性能ががた落ちになります。
このような場合、部分再帰的にFFTを書きなおすと劇的に性能が向上します。N点の部分再帰的FFTは次のような構造を持ちます。
最初にNを再帰分割するメリットは、L2サイズが不明であっても常に高い性能をえられることです。再帰分割によってL2キャッシュにデータがすっぽり入る大きさまで配列を分割すると、その先は常にL2キャッシュがヒットします。そしてサイズがKまで小さくなると、あとはL1上で残りのFFTが実行されます。現在のx86ではKは512か1024くらいに取れば良いでしょう。
再帰分割サイズがL2に入りきらない領域であっても、アクセスはアドレスに対して連続的なのでキャッシュのバースト・フィルが生かされます。
なお、フラットなFFTアルゴリズムの後半で内側の二つのループのネストを入れ替えるとループ・オーバーヘッドを減らすことができます。
部分再帰的FFTは、ムーアの法則の恩恵をたっぷり受けることのできる、優れたアルゴリズムです。
