曲がった迷路その37 - 複素数の四則演算をvDSPにやらせる [曲がった壁を持つ迷路の生成]
曲がった壁を持つ迷路のとりあえずのアプリのFourier変換にfftwライブラリを使ったけど、fftwはデカ過ぎて全体の5割を超えてしまった。容量的に本末転倒なので、FFTを自分で実装することにした。
vDSPを使ったDFT
ここまではなにも頭を使っていない。ということでvDSPを使ってみる。vDSPはAccelerate frameworkに含まれていてOSの一部としてすべてのMacにインストールされているライブラリである。あまり難しいことはできないけど、ライブラリの内部では可能ならAltiVecやSSEのベクタユニットを呼ぶようになっているという(ちなみに、「Taking Advantage of the Accelerate framework」にあるvDSPへのリンクは日本語だけでなく英語版でも切れている)。
vDSPのルーチンは例えばAppleのガイド(リンク先はpdf。html版が見つからなくなっている)から引き写すと
void vDSP_vmul( float *input_1, /* input vector 1 */ SInt32 stride_1, /* address stride for input vector 1 */ float *input_2, /* input vector 2 */ SInt32 stride_2, /* address stride for input vector 2 */ float *result, /* output vector */ SInt32 strideResult, /* address stride for output vector */ UInt32 size /* real output count */ );となっていて、これはfloatのベクトルinput_1とinput_2のかけ算の結果をresultに入れて返す、というもの。sizeにベクトルの長さ、stride_1、stride_2、strideResultに要素の間隔を入れる。
これを
float a[1024], b[1024], c[1024]; vDSP_vmul( a, 1, b, 1, c, 1, 1024 );というふうに使う。
vDSPにはこういった関数がたくさん定義されていて、このスタイルでfloat用、double用、そして一部複素数用、それから2のベキの長さの1次元と2次元のFFTも含まれている。
今回はFFTではなく、複素数のベクトルの演算を使ってみる。 vDSPではC99の複素数ではなく独自の定義をしている。vDSPであつかう複素数にはふた通りあって
struct DSPComplex { float real; float imag; }; typedef struct DSPComplex DSPComplex;というC99複素数と互換性のある(キャストすれば使える)ものと
struct DSPSplitComplex { float * realp; float * imagp; }; typedef struct DSPSplitComplex DSPSplitComplex;というふうに実部と虚部が別の配列になっているのとがある。ここにあげたのはfloat幅の複素数だけどもちろんdouble幅のもまったく同じ形式で定義されている。vDSPで計算できるのは別々の配列になったDSPSplitComplex型だけで、C99複素数と互換性のある型の配列とは相互に変換する関数が用意されている。
これがまたややこしくてバグの原因になる。けど、まあしかたがない。
ちなみに、fftwも
#include <complex.h> #include <fftw3.h>としてC99複素数のヘッダを先に読み込むとデータのプロトタイプ宣言はC99複素数になる。中身は同じだけど、C99の実数-複素数間の型変換や虚数単位の「I」とかが使えて記述がすっきりする。
2009-08-25 22:11
nice!(0)
コメント(0)
トラックバック(0)
コメント 0