vDSPについてのメモ - 2

会社がお盆休みで、横浜に今朝帰ってきた。のぼりの新幹線は空いていた。こないだのvDSPについてのメモの続き。

vDSPでなにができるのか

vDSPのReferenceでは1次元配列をvector、2次元配列をmatrix、配列ではない一つの数をscalarと呼んでいる。

vDSPの関数はだいたいこんな呼び方をする。

void vDSP_vadd (const float input1[],
   vDSP_Stride stride1,
   const float input2[],
   vDSP_Stride stride2,
   float result[],
   vDSP_Stride strideResult,
   vDSP_Length size);

これはふたつの不動小数点の配列input1とinput2を要素ごとに足してその結果を別の配列resultに入れる、という関数である。要素の数は最後の引数のsizeで指定する。vDSP_Strideは、パックされて詰った要素を一通り全部ではなく、いくつかおきに計算したい場合を考えている。これはCの2次元配列で縦方向に計算したい場合などを考慮している。要素が詰まっていてすべての要素をはじめから全部処理するにはvDSP_Strideに1を指定する。他の、例えば2を指定するとひとつおきになる。

vDSPの関数

vDSPにはごちゃごちゃといろんな関数があって探すのが面倒だけど、大まかに

• 入力がvectorで出力もvector（さっきの配列の要素同士の足し算など）
• 入力がvectorで出力はscalar（内積演算など）
• 入力がvectorとscalarの両方（vectorの定数倍など）
• 1次元、および2次元のFFT
• 1次元の畳み込みおよび自己相関、ウィンドウイング
• 実-複素配列変換
• matrixの積

ができる。この中ではFFTのルーチンを昔、僕は使っていたけど、要素数に制限がないのにパフォーマンスはvDSPとあまり変わらないのでfftwの方に移行してしまった。他を使うのは今回初めて。

演算ごとに関数があってわかりにくいが、名前はだいたい

vDSP_{次元}{演算}{精度}

のようになっている。 次元はvectorではv、matrixではm、scalarではsがつけられている。演算は例えば足し算はadd、かけ算はmul、以下同様でほかに足してから別のをかけるとかではam（addしてmul）などとなる。精度は配列の要素がdoubleの場合はD、32ビット整数の場合はi、などがつく。単精度不動小数点の場合は省略される。

たとえば

vDSP_vamD()

はふたつの倍精度の配列を要素ごとに足してそのあと別の配列を要素ごとにかける関数である。 基本的な四則演算には複素数に対応したものもあってそのときは次元を表す文字の前にzがつく。例えば

vDSP_zvadd()

は複素数のふたつの配列を足す関数である。ただし、vDSPで使える複素数はDSPSplitComplexという形式に限られる。これは

struct DSPSplitComplex {
  float *             realp;
  float *             imagp;
};
typedef struct DSPSplitComplex          DSPSplitComplex;

のように定義されて実数部と虚数部が別々のところにパックされた配列になっている。普通の（interleaved、実数部と虚数部が組になった）複素数の配列とDSPSplitComplexと相互変換する

vDSP_ctoz()
vDSP_ztoc()
vDSP_ctozD()
vDSP_ztocD()

が用意されている。実装上の理由からこうなっているんだと思うけど、大きな配列を並び替えるときは場所塞ぎでちょっとした計算（ふたつのベクトルを足すだけとか）ではオーバーヘッドが馬鹿にならない。

何ができるかはだいたいわかった。vDSP Reference Collectionなんていう、いかにも散らばってたのをかき集めました、というようなドキュメントがぽつ、とある。まあ、これを見れば全体はわかる。それではベンチマーク行ってみよう。