vDSPによる2の基数のFFTの実装

2の基数のFFTはすでに独自に実装を終えている。しかしこの実装はまったく効率を考えていないし、一方でMacOS Xには古い（と言うとネガティブに聞こえるかもしれないので枯れた、と言うべきか）vDSPライブラリがあってこれはCPUのベクタユニットに最適化さているという。どのくらい枯れた最適化が有効なのかベンチマークの意味でこれを組み込むことにする。

vDSPによる2の基数のFFT

vDSPにも2の基数のFFTがあり、fftwのサイトでのベンチマークでは、ちょっと古いデータだけどPowerPCでfftwと同程度に速いという結果になっている。

これを期待して2の基数のFFTはこちらも使えるようにする。

vDSPでは

• 実、および複素数
• 1次元、2次元
• float配列、double配列
• in-place、out-of-place

のすべての組み合わせと、out-of-place型の3×2ⁿと5×2ⁿの1次元複素FFTがある。

前にも書いたが複素配列の場合はsplit型と言う実数部と虚数部が別の配列になったものしか受けない。そのため普通の（interleavedと呼んでいる）配列と相互変換する関数も用意されている。

2次元FFTもあるけど、例えば縦は2のベキだけど横は素数とか言う場合にvDSPが使えなくなるので、今回は1次元のルーチンだけを使って2次元化は単に繰り返すだけにする。今回使うのは

• 1次元
• 複素数
• float幅
• in-place

のものだけ。この関数は

void vDSP_fft_zip (FFTSetup setup,
   DSPSplitComplex * signal,
   vDSP_Stride stride,
   vDSP_Length log2n,
   FFTDirection direction);

というもの。FFTSetupというのはfftwのplan構造体のようなもので、回転因子や中間結果のバッファなんかを確保しているらしいが、いつものOpaque型で中身はわからない。signalが泣き別れ型の複素配列で、関数から返ってくると結果が上書きされている。他は名前の通り。

これを使うには

1. FFTSetup構造体を作る
2. 入力配列をDSPSplitComplex型に変換する
3. vDSP_fft_zip()を呼ぶ
4. 結果の配列をDSPSplitComplex型から変換する
5. FFTSetup構造体を破棄する

ということになる。

これを使って、前の抽象クラスのサブクラスとして実装する。

@interface ASvDSPFFT1D : ASDFT1D {
    unsigned int    log2n;
    FFTSetup        setup;
    DSPSplitComplex splitArray;
    FFTDirection    direction;
    float           *inverseBuffer;
}
- (id)initWithLog2OfPoints:(unsigned int)log2OfLength andDirection:(int)dir;
@end

こういうふう。 実装のまず初期化は

- (id)initWithLog2OfPoints:(unsigned int)log2OfLength andDirection:(int)dir
{
    log2n = log2OfLength;
    self = [super init];
    direction = dir > 0 ? FFT_FORWARD : FFT_INVERSE;
    splitArray.realp = (float *)malloc((1 << log2n) * sizeof(float));
    splitArray.imagp = (float *)malloc((1 << log2n) * sizeof(float));
    setup = vDSP_create_fftsetup(log2n, FFT_RADIX2);
    if (direction == FFT_INVERSE)
        inverseBuffer = (float *)malloc((1 << log2n) * sizeof(float));
    else
        inverseBuffer = NULL;
    return self;
}

読んで字のごとく。

FFTの実行は

- (void)performForComplexArray:(float complex *)dat withStride:(int)stride
{
    unsigned int    N = 1 << log2n;
    vDSP_ctoz((const DSPComplex *)dat, stride * 2, &splitArray, 1, N);
    vDSP_fft_zip(setup, &splitArray, 1, log2n, direction);
    if (direction == FFT_INVERSE)
        [self unnormalizeBack];
    vDSP_ztoc((const DSPSplitComplex *)(&splitArray), 1,
                                     (DSPComplex *)dat, stride * 2, N);
}

というもの。配列を変換して、FFTを呼んで、戻して、ということをしているだけ。

注意が必要なのは配列変換のvDSP_ctoz()とvDSP_ztoc()のstrideで、普通の配列側はfloat二つ分なので2倍しないといけない。そんなの当然だから関数内部で調整してくれればいいのに、そうなっていない。こんなところに落とし穴があるなんて思わなかったのでドキュメントを見ずにそのまま使ってたらハマった。

FFTの逆変換のときのunnormalizeBackというメソッドはfftwと互換をとるためのもので、vDSPのFFTは逆変換のときに規格化しているとマニュアルに書いてある。

それをやめるもので

- (void)unnormalizeBack
{
    unsigned int    N = 1 << log2n;
    float           mul = 1.0f / N;
    vDSP_vsmul((const float *)(&(splitArray.realp)), 1,
                                            &mul, inverseBuffer, 1, N);
    memcpy(splitArray.realp, inverseBuffer, N * sizeof(float));
    vDSP_vsmul((const float *)(&(splitArray.imagp)), 1,
                                             &mul, inverseBuffer, 1, N);
    memcpy(splitArray.imagp, inverseBuffer, N * sizeof(float));
}

というように要素全部をN倍している。非常に間抜けだけどしかたがない。

しかし実際これで計算してみると逆変換では値が大きくなってしまう。試しにこのメソッドをやめるとfftwの結果と一致した。つまりvDSPのマニュアル には逆変換側はNで全要素を割り算しているとなっているが、vDSPは正規化しないライブラリであるということだ。まったくもってひどい。