曲がった迷路その38 - vDSPを使ったDFT

軽く独自実装することに決めた離散Fourier変換。定義式そのままのFourier変換をAppleのAccerelate frameworkにあるvDSPライブラリを使って書いてみる。

愚直DFTをvDSPで書き直す

さてこれで前々回と同じ1次元のDFTを書いてみよう。どこでvDSPを使うかと言うと一番内側のループ。

void    dftf5(DSPSplitComplex *cin, DSPSplitComplex *cout, float complex *wcoef, int length, int stride)
{
    DSPSplitComplex tempc;
    DSPSplitComplex wc;
    wc.realp = (float *)malloc(length * sizeof(float));
    wc.imagp = (float *)malloc(length * sizeof(float));
    tempc.realp = (float *)malloc(length * sizeof(float));
    tempc.imagp = (float *)malloc(length * sizeof(float));
    int i;
    
    for (i = 0 ; i < length ; i ++) {
        float   re, im;
        alignWCoef(&wc, wcoef, length, i);      //(1)
        vDSP_zvmul(cin, stride, &wc, 1, &tempc, 1, length, 1);//(2)
        vDSP_sve(tempc.realp, 1, &re, length);  //(3)
        vDSP_sve(tempc.imagp, 1, &im, length);
        cout->realp[i * stride] = re;           //(4)
        cout->imagp[i * stride] = im;
    }
    free(wc.realp);
    free(wc.imagp);
    free(tempc.realp);
    free(tempc.imagp);
}

外からexponentialの係数をもらうのは一緒。 (1)のalignWCoef()という関数はもらったexponentialの係数を、dftf2()でmoduloをとるのと同じことをしながらDSPSplitComplexの形にする。

void    alignWCoef(DSPSplitComplex *outw, float complex *cw, int length, int stride)
{
    float   *realp = outw->realp;
    float   *imagp = outw->imagp;
    int i;
    for (i = 0 ; i < length ; i ++) {
        float complex   tmp = cw[(i * stride) % length];
        *(realp ++) = crealf(tmp);
        *(imagp ++) = cimagf(tmp);
    }
}

(2)で入力ベクトルとexponential係数のかけ算をやって、(3)で足しあわせている。ベクトルの要素を足しあわせるvDSP関数は実数ベクトル用しかないのでそれを実部と虚部それぞれやって、出力ベクトルに(4)で代入している。

とても速くなるようには見えないんだけど、まあやってみよう。

ところで前回あげたvDSPのpdfドキュメントには間違いが多い。関数の名前が一文字違うだけで全然別物になるのでこういう短い名前は気をつけないといけないのに、明らかにドキュメントを起こすときにコピペのミスをいたるところでやってる。これにはまいりました。例えば今回使ったvDSP_zvmul()関数。33ページのVector-to-Vector Arithmetic Operationsの表2-6に載ってない。かわりにvDSP_zvcmul()というふたつ目のベクトルの共役をとってかけ算するという関数が載っている。僕は一文字しか違わなくて気がつかずにこの表から関数をコピペした。DFTの正変換と逆変換が逆になってこれがなんでかわかるまでずいぶん悩んだ。

トリビアルなものではあるけど間違いはこれだけではないので読む人は要注意。

これを使った2次元DFTは

void    dft2dimf3(DSPSplitComplex *cin, DSPSplitComplex *cout, int row, int column)
{
    float   *realp = (float *)malloc(row * column * sizeof(float));
    float   *imagp = (float *)malloc(row * column * sizeof(float));
    DSPSplitComplex tin, tout;
    float complex   *wcoef = allocWCoef(row, forward);
    int r, c;
    for (c = 0 ; c < column ; c ++) {
        tin.realp = cin->realp + c * row;
        tin.imagp = cin->imagp + c * row;
        tout.realp = realp + c * row;
        tout.imagp = imagp + c * row;
        dftf5(&tin, &tout, wcoef, row, 1);
    }
    free(wcoef);
    wcoef = allocWCoef(column, forward);
    for (r = 0 ; r < row ; r ++) {
        tin.realp = realp + r;
        tin.imagp = imagp + r;
        tout.realp = cout->realp + r;
        tout.imagp = cout->imagp + r;
        dftf5(&tin, &tout, wcoef, column, row);
    }
    free(wcoef);
    free(realp);
    free(imagp);
}

見た目にはどんどん遅くなるような気がしてしょうがない。