落胆がらくた街

多分三日坊主で終わると思うんですけど

2PIでループする位相変数をUINTに正規化すると美味い

http://www.musicdsp.org/archive.php?classid=1
これの Fast SIN approximation for usage in e.g. additive synthesizers を試してみた。

float fastsin(UINT32 phase)
{
    const float frf3 = -1.0f / 6.0f;
    const float frf5 = 1.0f / 120.0f;
    const float frf7 = -1.0f / 5040.0f;
    const float frf9 = 1.0f / 362880.0f;
    const float f0pi5 = 1.570796327f;
    float x, x2, asin;
    UINT32 tmp = 0x3f800000 | ( phase >> 7 );
    if (phase & 0x40000000)
        tmp ^= 0x007fffff;
    x = ( *( (float*)&tmp ) - 1.0f ) * f0pi5;
    x2 = x * x;
    asin = ( ( ( ( frf9 * x2 + frf7 ) * x2 + frf5 ) * x2 + frf3 ) * x2 + 1.0f ) * x;
    return ( phase & 0x80000000 ) ? -asin : asin;
}

検証コード↓

int main() {
    timerInit();

    float delta_f = 2 * PI * 1000 / SRATE;
    float phase_f = 0;
    UINT32 delta_i = 1000.0 / SRATE * UINT32_MAX;
    UINT32 phase_i = 0;
    

    //特に工夫無く、位相phaseと角速度deltaをfloatにしてmath.hのsinfを呼ぶ
    double t = 0;
    float test = 0;

    for (size_t i = 0; i < SRATE; i++)
    {
        timerStart();

        test += sinf(phase_f);
        phase_f += delta_f;

        t += timerStop();

        if (phase_f >= 2 * PI)phase_f -= 2 * PI; //2PIで位相がループするのはとりあえず時間計測の外とする
    }

    cout << "sinf(float)::" << "total:" << t << "[ms]\n";

    //UINT32正規化版
    t = 0;
    float test2 = 0;

    for (size_t i = 0; i < SRATE; i++)
    {
        timerStart();

        test2 += fastsin(phase_i);
        phase_i += delta_i;

        t += timerStop();
    }

    cout << "sinf(UINT32)::" << "total:" << t << "[ms]\n";

    cout << "test:" << test << ", test2:" << test2 << endl; //何らかの形で変数を参照してやらないと最適化で省かれる

}

結果。

sinf(float)::total:3.06791[ms]
sinf(UINT32)::total:0.292572[ms]
test:0.0168861, test2:6.52373e-05

UINT32、速い(確信)。
ちなみに計測環境はWin8.1、cpu:intel 4690k、DRAM 16GBでコンパイル環境はVisual studio 2015のVisual C++、x64 Releaseビルドで最適化レベル最大。
sinの出力結果の総和が微妙に異なってるのは精度の問題なんだが、別に精度計測のソース書いたら1サンプル毎の精度は十分だったので割愛。


位相をUINT32に正規化するメリットを考えてみる。

  • 検証コードでは時間計測から省いたけど、真面目に2PIでループさせようとする場合1サンプル毎に位相が2PI超えたら戻す処理がいる。UINTならオーバーフローするので不要。
  • 位相を元にテーブルを引きたいみたいな時、位相変数の上何桁と下何桁を分けるみたいな処理が高速で出来る float->intキャストよりずっとはやい!
  • fastsinが速い
  • 位相変数そのものの計算もビット演算が使える分そこそこ速い。特に頻出の2n倍が高速なのは嬉しい

とかこんな感じかな。結構嬉しいと思うんだけど、逆にデメリットはどんなんあるんだろう。誰か教えてください。