SSブログ

曲がった迷路その35 - 任意基数のFourier変換 [曲がった壁を持つ迷路の生成]

曲がった壁を持つ迷路はこないだとりあえず一段落したのでキリをつけようと思っていた。でもアプリのサイズをよく見るとfftwのライブラリがリソースを含めて自分で書いた分より大きい。ちょっとこれは恥ずかしい。ということでFFTに関して追加する。
どうでもいいけどFFTでググるとファイナルファンタジータクティクスとかいうのがトップにヒットする。あ〜あ。

きっかけ

こないだの迷路アプリではスペックルパターンを得るためにFFTをする必要が出て、fftwライブラリを使った。
fftwは
  • フリー
  • 高速
  • 移植性が高い
  • 多次元、任意点数対応
で強力なんだけど非常に大きなライブラリになっている。

ほかにもFFTのライブラリはたくさんあるけど多くがFFTのオリジナルの2のベキの点数だけに対応するものが多い。
今回のアプリは
  • 2次元
  • 任意の縦横比にしたい、従って任意点数
  • FFTは初期化で使われるだけで、アプリ全体の効率を律速することはない
であった。最初のふたつの特徴からfftwを使うことにしたけど処理効率は問われない。ということでサイズの大きなfftwライブラリのかわりに
  1. 任意点数の2次元データの外に0を追加して点数が2のベキになるようにして普通のライブラリを使う
  2. Fourier展開の式を愚直に計算する独自ルーチンを作る
ということが考えられる。

ひとつめの0詰めは音声信号の解析などではよくやるパターンである。この場合、領域の端では0を埋めなかったときの結果と異なってくる。

今回の迷路生成ではFFTで作られたスペックルパターンがトーラスにラップしている(上と下、右と左が繋がっている)とみなして端っこの処理をズルしているので、これも一緒に変更しないと使えない。それはもうめんどう。
ということで、残ったふたつめの方法を考えてみる。

これは式をそのままコードに展開すればいいだけなので簡単だけど、2次元であること、複素数を扱わなければならないことにちょっと頭を使う必要がある。

ということで、どうせならfftwの効率を体感するために
  1. C99複素数を使って式をそのまま展開
  2. MacOS XのAccelerate frameworkにあるvDSPルーチンを使ってみる
として比較してみることにする。

離散Fourier変換

FFTみたいに、連続量ではない数値配列のFourier変換(サンプリング列に対するFourier展開というべきか)を一般に、離散Fourier変換(Discrete Fourier Transform DFT)という。N個の要素を持った複素数値配列fnのDFTFk
0820eq1.png
となる。ここでiは虚数単位である。まぎらわしい。N個の点の値を得るためにN個のかけ算をするのでN2回のかけ算をすることになる。

0820eq2.png
の性質からNがたくさんの小さな素因数に分解できるときはModulo演算によって現れるexponentialの係数の種類が減って、それが規則正しく並ぶのでそれをFFTでは利用して演算回数を減らしている。

昔自分でFFTを実装して、これを考えついた人は頭いいなあと感心した記憶がある。ちなみに、当時は回折計算に使うためにMC68020を積んだHPのワークステーション用に、まずCで書いた。けっきょくそれは遅くて使い物にならなかった。そのあと68K+68881(浮動小数点コプロセッサ。「なんですかそれは?」と言われましても...)のアセンブラで、しかも68020の256バイト(!256kバイトではない!)の命令キャッシュに一番内側のループが収まるように苦労して書いた。若い頃に集中してコードを書く勉強にはなったけど、大変な労力を費やした。

Nが2のベキなら簡単で一番速くなる(かけ算回数はNlog2N)けど、Nを素因数分解したときにそれぞれの素因数が小さな値のときには同じ手を使って高速化できる。fftwライブラリはNが13までの素数の積になるときにこの手法を使っていると書いてある。しかしそのアルゴリズムはウンザリするほどめんどくさい。最終的な結果を並び替える(2のベキの場合が一番簡単で2進ビット反転順になる)汎用のアルゴリズムを考えるだけでも気が遠くなる。もちろんfftwはそれをやっている。

昔は足し算コストよりもかけ算コストの方が高かったのでかけ算の回数だけを最適化していたが、最近のCPUは同じスループットになっていたり、足し算とかけ算を同時に実行する命令(D=A×B+CがC=A×Bと同じスループットで実行できる)を持っていたりするので、逆にアルゴリズムがそういった特徴に最適化されるようにもなってきた(例えばかけ算と足し算のコストが同じ場合、基数を4にとる方が速くなる)。さらにデータキャッシュの方式やデータサイズとの大小関係でもアルゴリズムによって有利不利が起こる。fftwではそういったことをひっくるめてfftw_planを作るときに実際に計算してみて最適なアルゴリズムを選択できるようになっている。

また、Nが2のベキ以外のときオリジナルのFFTとは違った手法で最終的に2のベキのFFTに帰着させるDFTの高速化がいろいろある。

特にNが素数の場合、式-2によって、式-1のどのkに対してもexponentialの係数が0からN-1までまんべんなく現れることを利用してN-1のDFT3回に変換するRaderのアルゴリズムというのもある。これは素数の超基本的な性質の応用で、暗号の分野の屋台骨でもある。このアルゴリズムは最近知ったがいろいろなことを考えるねえ、面白い。

こういったDFTの高速化は、Nが100を超えると計算量は桁が変わるぐらい違うのでめちゃめちゃ強力である。音声のフィルタリングみたいな、長いデータに対してFFTが使えることでリアルタイム性が早い時期に確保できるようになった分野もある。

今回は高速化が目的ではないので、こういったところには突っ込まないつもり。

nice!(0)  コメント(1)  トラックバック(0) 

nice! 0

コメント 1

Viagra vs cialis vs levitra

沒有醫生的處方
when will generic cialis be available http://cialisyoues.com/ 5 mg cialis coupon printable
by Viagra vs cialis vs levitra (2018-04-14 03:47) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

献立08/20献立08/21 ブログトップ

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。