Pebble Coding

ソフトウェアエンジニアによるIT技術、数学の備忘録

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つも違いますね。
任意の要素数フーリエ変換はN2回の乗算とN(N-1)回の加算が必要ですが、
素数が2の冪数の場合は効率のよいアルゴリズムが存在しN/2 * log2(N-1)回の乗算とN * log2(N)回の加算で済みます。

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

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

vDSPの内積計算がどのくらい速いのか検証してみた

vDSPの内積計算がどのくらい速いのか検証してみました。

#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#endif

int main(int argc, const char * argv[]) {
    float* a = new float[4096];
    float* b = new float[4096];
    float* c = new float;

    for (int i = 0; i < 4096; ++i) {
        a[i] = (float)i / 4096.f;
        b[i] = (float)i / 4096.f;
    }

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

    for (int i = 0; i < 4096; ++i) {
#if 0
         vDSP_dotpr(a, 1, b, 1, c, 4096);
#elif SSE
        // omit
#elif AVX
        // omit
#else
        *c = 0;
        float* p = a;
        float* q = b;
        const float* end = p + 4096;
        while (p != end) {
            *c += *(p++) * *(q++);
        }
#endif
    }

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

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

    // M1 -ONone(最適化なし) : 58482 us
    // M1 -Os(最適化) : 6023 us
    // M1 vDSP: 975 us
    // Intel /Od(最適化なし) : 30764 us
    // Intel /Ox(最適化) : 7521 us
    // Intel SSE: 3578 us
    // Intel AVX256: 851 us

    printf("%f\n", *c);

    delete [] a;
    delete [] b;
    delete c;

    return 0;
}

Intelは第12世代IntelCPUで動作するWindowsマシンを使いました。
M1はM1 Pro CPUです。
vDSPはIntelCPUでAVX命令を使ったときと同じくらいの速度が出るということがわかりました。
追加で内積以外の結果もみてみました。

M1 vDSP_conv: 123ms
M1 -Os: 2182ms

std::partition_point

C++11 で導入されたstd::partition_pointはラムダ式を渡すときに
lower_bound,upper_boundに比べて引数が一つ少なくて済むというメリットがある。

const std::vector<int> v = { -50, -20, -10, -7, -1, -2, 21, 22, 26, 27 };

const int val1 = 21;

// lower_bound
// 21以上の値が現れる最初の位置
auto it1 = std::partition_point(v.begin(), v.end(),
[val1](const auto& obj) { return obj < val1; });
printf("lower_bound %d\n", *(--it1));
// lower_bound -2

// upper_bound
// 21より大きな値が最初に現れる位置
auto it2 = std::partition_point(v.begin(), v.end(),
[val1](const auto& obj) { return obj <= val1; });
printf("upper_bound %d\n", *(--it2));
// upper_bound 21

const int val2 = 20;

// lower_bound
// 20以上の値が最初に現れる最初の位置
auto it3 = std::partition_point(v.begin(), v.end(),
[val2](const auto& obj) { return obj < val2; });
printf("lower_bound %d\n",  *(--it3));
// lower_bound -2

// upper_bound
// 20より大きな値が最初に現れる位置
auto it4 = std::partition_point(v.begin(), v.end(),
[val2](const auto& obj) { return obj <= val2; });
printf("upper_bound %d\n", *(--it4));
// upper_bound -2

256bit(32バイト)にアラインされたメモリを取得する(C++17)

C++17以上ではnewでアラインされたメモリを簡単に取得できるようになりました。

    for (int i = 0; i < 10000; ++i) {
        float* a = new (std::align_val_t{32}) float;
        
        auto val = reinterpret_cast<size_t>(a);
        if ((val % 32) != 0) {
            assert(false);
        }
    }

余談ですが64bit環境のWindows VC++ではnewは16バイト、 macOS10.6以上の64bit環境のmallocでは16バイトにアラインされるようです。

https://docs.microsoft.com/ja-jp/cpp/c-runtime-library/reference/malloc?view=msvc-170

https://developer.apple.com/library/archive/documentation/Cocoa/Conceptual/Cocoa64BitGuide/OptimizingMemoryPerformance/OptimizingMemoryPerformance.html

C++20で気になる機能

  • 2 つの値の中点を計算する関数が追加された。
  • 浮動小数点型のatomic操作関数が追加された。
  • 円周率の定義が使えるようになった。
  • コンテナにcontains()関数が追加。
  • 2 の累乗数に関する関数が追加された。
  • std::make_shared()に配列サポートが追加された。
  • 定数式での仮想関数呼び出しが可能.
struct IValue {
    virtual int value() const = 0;
};

struct Value : public IValue {
    constexpr int value() const override {
        return 2022;
    }
};


int main(int argc, const char * argv[]) {

    constexpr Value v;
    static_assert(v.value() == 2022);
}
  • std::format
  • コンセプト swiftのprotocol whereに近いことができるようになった。
  • ranges (今の所MSVCのみ対応,clang,GCC非対応)