vDSPのフーリエ変換がどのくらい速いのか検証してみた

#include <Accelerate/Accelerate.h>

static void dft_manual() {
    const auto N = 64;
    const auto xr = std::make_unique<float[]>(N);
    const auto xi = std::make_unique<float[]>(N);
    const auto yr = std::make_unique<float[]>(N);
    const auto yi = std::make_unique<float[]>(N);
    auto x_re = xr.get();
    auto x_im = xi.get();
    auto y_re = yr.get();
    auto y_im = yi.get();

    const auto t1 = std::chrono::system_clock::now();

    for (int m = 0; m < 4096; ++m) {
        // input
        for (int i = 0; i < N; ++i) {
            x_re[i] = (float)i / (float)N;
            x_im[i] = 0;
        }

        for (int k = 0; k < N; ++k) {
            y_re[k] = 0;
            y_im[k] = 0;
            for (int i = 0; i < N; ++i) {
                const auto w_re = cos(2 * M_PI * k * i / N);
                const auto w_im = - sin(2 * M_PI * k * i / N);
                y_re[k] += (x_re[i] * w_re - x_im[i] * w_im);
                y_im[k] += (x_re[i] * w_im + x_im[i] * w_re);
            }
        }
    }

    const auto t2 = std::chrono::system_clock::now();

    double elapsed = std::chrono::duration_cast<std::chrono::microseconds>(t2-t1).count();
    printf("%f\n", elapsed);

    // output
    for (int k = 0; k < N; ++k) {
        printf("%d %f + %f j\n", k, y_re[k], y_im[k]);
    }
}

static void dft_vdsp() {
    const auto N = 64;
    const auto xr = std::make_unique<float[]>(N);
    const auto xi = std::make_unique<float[]>(N);
    const auto yr = std::make_unique<float[]>(N);
    const auto yi = std::make_unique<float[]>(N);
    DSPSplitComplex x;
    x.realp = xr.get();
    x.imagp = xi.get();
    DSPSplitComplex y;
    y.realp = yr.get();
    y.imagp = yi.get();

    for (int i = 0; i < N; ++i) {
        x.realp[i] = (float)i / (float)N;
        x.imagp[i] = 0;
    }

    // 64 = 2 ** 6
    const auto length = 6;
    auto fftSetup = vDSP_create_fftsetup(length, FFT_RADIX2);

    const auto t1 = std::chrono::system_clock::now();

    for (int m = 0; m < 4096; ++m) {
        vDSP_fft_zop(fftSetup, &x, 1, &y, 1, length, FFT_FORWARD);
    }

    const auto t2 = std::chrono::system_clock::now();

    double elapsed = std::chrono::duration_cast<std::chrono::microseconds>(t2-t1).count();
    printf("%f\n", elapsed);

    // output
    for (int k = 0; k < N; ++k) {
        printf("%d %f + %f j\n", k, y.realp[k], y.imagp[k]);
    }

    vDSP_destroy_fftsetup(fftSetup);
}

// M1 manual -O0 254142us
// M1 manual -Os 149266us
// M1 vdsp 252us

manualというのはまったく最適化せずに離散フーリエ変換(Discrete Fourier Transform)のアルゴリズムを実装したものです。
vDSP_fft_zop()は高速フーリエ変換(Fast Fourier Transform)で高速化アルゴリズムになっていることもあり桁が3つも違いますね。
任意の要素数のフーリエ変換はN²回の乗算とN(N-1)回の加算が必要ですが、
要素数が2の冪数の場合は効率のよいアルゴリズムが存在しN/2 * log2(N-1)回の乗算とN * log2(N)回の加算で済みます。

• 実高速フーリエ変換(real fast fourier transform)

ここでもやっていますが、フーリエ変換する際、入力データの虚数部分は全て0を設定します。
N個の入力データを前半と後半に分けて、N/2個の複素数にして入力し、出力を調整すれば、同じ結果が半分の長さの変換で結果が得られることが知られています。
これを実高速フーリエ変換と呼びます。
この場合はvDSP_fft_zop()ではなくvDSP_fft_zrop()を使用します。

N個の実値をN/2個の複数に変換するのがvDSP_ctoz()です。