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倍が高速なのは嬉しい
とかこんな感じかな。結構嬉しいと思うんだけど、逆にデメリットはどんなんあるんだろう。誰か教えてください。