[Перевод] Трюк с тригонометрией06.08.2019 17:19

Скорее всего, вам известны следующие соотношения еще со школы:

$\sin(\alpha + \beta) = \sin\alpha \times \cos\beta + \cos\alpha \times \sin\beta \\ \cos(\alpha + \beta) = \cos\alpha \times \cos\beta - \sin\alpha \times \sin\beta$

Когда вы в детстве впервые познакомились с этой формулой, скорее всего, вашим первым чувством была боль из-за того, что эту формулу надо запомнить. Это очень плохо, потому что на самом деле вам не нужно запоминать эту формулу — она сама выводится, когда вы поворачиваете треугольник на бумаге. На самом деле, я делаю то же самое, когда записываю эту формулу. Это толкование будет очевидным к середине этой статьи. Но сейчас, чтобы оставить все веселье на потом и отодвинуть момент, когда вы скажете «Эврика!», давайте подумаем, а зачем нам вообще задумываться об этой формуле.

Тригонометрические функции sin() и cos() возможно самые популярные в компьютерной графике, поскольку они являются основой для описания любой круглой формы параметрическим способом. Среди мест их возможного применения генерация кругов или объемных объектов вращения, при вычислении преобразования Фурье, процедурная генерация волн на плоскости воды, генераторы для программного синтезатора звука, и так далее. Во всех этих случаях sin() и cos() вызываются внутри цикла, как тут:

for(int n=0; n < num; n++)
{
    const float t = 2.0f*PI*(float)n/(float)num;
    const float s = sinf(t);
    const float c = cosf(t);

    // do something with s and c
    ...
}

Мы начинаем переписывать цикл инкрементальным образом (см. код ниже), так что нам легче представить, что на итерации n данного цикла с фазой t, следующая итерация, n+1, будет считать sin() и cos() для t+f. Другими словами, у нас сосчитаны sin(t) и cos(t) и нам надо сосчитать sin(t+f) и cos(t+f):

const float f = 2.0f*PI/(float)num;
const float t = 0.0f;
for( int n=0; n < num; n++ )
{
    const float s = sinf(t);
    const float c = cosf(t);

    // do something with s and c
    ...
    t += f;
}

Неважно, каким образом мы получили t и каков ее диапазон значений (в примере выше — $[0;2\pi]$ ). Единственное, что нас беспокоит, так это то, что есть цикл, который постоянно вызывает sin() и cos() с параметром, который увеличивается в постоянных шагах (в данном случае, $\frac{2\pi}{\text{num}}$ ). Эта статья о том, как оптимизировать этот код для скорости таким образом, что одни и те же вычисления могут выполняться вообще без использования функций sin() или cos() (во внутреннем цикле), и даже более быстрой функции sincos().

Но если посмотреть на первую формулу в статье, мы увидим, что если $f = \alpha$ и $t = \beta$ , мы можем переписать это как

sin(t+f) = sin(f)*cos(t) + cos(f)*sin(t)
cos(t+f) = cos(f)*cos(t) - sin(f)*sin(t)

или, другими словами

new_s = sin(f)*old_c + cos(f)*old_s
new_c = cos(f)*old_c - sin(f)*old_s

Так как f — константа, то sin(f) и cos(f) тоже. Назовем их a и b соответственно:

new_s = b*old_c + a*old_s
new_c = a*old_c - b*old_s

Это выражение может быть напрямую использовано в нашем коде, и тогда мы получим цикл, в котором не вызывается дорогих (да вообще никаких) тригонометрических функций!

const float f = 2.0f*PI/(float)num;
const float a = cosf(f);
const float b = sinf(f);
float s = 0.0f;
float c = 1.0f;
for( int n=0; n < num; n++ )
{
    // do something with s and c
    ...

    const float ns = b*c + a*s;
    const float nc = a*c - b*s;
    c = nc;
    s = ns;
}

К настоящему моменту мы слепо играли с математикой не понимая, что же на самом деле происходит. Давайте перепишем внутренний цикл так:

$s_{n+1} = s_na + c_nb\\ c_{n+1} = c_na - s_nb$

Некоторые из вас могли заметить, что это — формула поворота объекта в двухмерном пространстве. Если вы все еще не поняли этого, возможно, матричная форма вам поможет.

$\left(\begin{matrix}s_{n+1} \\ c_{n+1}\end{matrix}\right) = \left(\begin{matrix}a & b \\ -b & a\end{matrix}\right) \cdot \left(\begin{matrix}s_{n} \\ c_{n}\end{matrix}\right)$

В самом деле, sin(t) и cos(t) можно сгруппировать в вектор длины 1 и отрисовать как точку на экране. Назовем этот вектор x. Тогда, $x = \{\sin\beta, \cos\beta\}$ . Значит, векторная форма выражения — $x_{n+1} = Rx_n$ , где $R = \left(\begin{matrix}a&b\\-b&a\end{matrix}\right)$ . Мы видим, что наш цикл выполняет небольшой двухмерный поворот каждую итерацию так, что x вращается по кругу во время выполнения цикла. Каждую итерацию x вращается на $\frac{2\pi}{\text{num}}$ радиан.
Итак, в основном,

$\sin(\alpha + \beta) = \sin\alpha \times \cos\beta + \cos\alpha \times \sin\beta \\ \cos(\alpha + \beta) = \cos\alpha \times \cos\beta - \sin\alpha \times \sin\beta$

это формула движения точки $x = \{\sin\beta, \cos\beta\}$ по окружности с шагом в $\beta$ радиан. Чтобы это сделать, мы построим одну из двух осей поворота, к примеру, $u = \{\cos\beta, \sin\beta\}$ . Первый компонент поворота — проекция $x$ на $u$ . Так как $x$ и $u$ нормализованы (имеют длину 1), проекция — их скалярное произведение. Следовательно, $s = x\cdot u = \sin\alpha\cdot\cos\beta + \cos\alpha\cdot\sin\beta$ , и конечно второй компонент — антипроекция, которую можно найти, спроецировав на перпендикулярную ось, $v$ . Мы можем создать этот вектор, развернув координаты $u$ и изменить знак на противоположный у первой координаты: $c = x\cdot v = \cos\alpha\cdot\cos\beta + \sin\alpha\cdot\sin\beta$

Обычно вы должны иметь возможность выполнять эти крошечные вращения снова и снова. В самом деле, $|R| = \left|\begin{matrix}a&b\\-b&a\end{matrix}\right| = a^2 + b^2 = \sin^2\alpha + \cos^2\alpha = 1$ , что означает, что матрица $R$ не увеличивает и не уменьшает пространство, к которому применена, что значит, что $x$ будет двигаться по идеальной окружности. Однако из-за того, что компьютеры не точны, $x$ будет двигаться по спирали и в конце концов совпадет с центром окружности вращения. У меня не возникало таких проблем, но я думаю, что они могут возникнуть при очень больших num, т.е. маленьких углах поворота.

В Kindercrasher, 4096-байтной демке из 2008 (скриншот на КДПВ), группа сфер пульсирует под музыку. Для этого я сосчитал преобразование Фурье звука. Существуют алгоритмы, делающие это в реальном времени, к примеру, FFT. Однако, мой код должен вместиться в несколько килобайт, и я решил пойти иным путем. Вместо реализации FFT, я написал DFT по его простому определению. Вы можете проверить это в википедии.

$X_k = \sum_{n=0}^{N-1}{x_ne^{-\frac{2\pi i}{N}kn}}\quad k=0,1,\ldots,N-1$

Моя функция также принимает 16-битный звуковой стерео буфер, x, и возвращает первые 128 частот звукового спектра звука y. Посмотрите, как организован внутренний цикл, тот, что выполняется 4096 раз: ни одного вызова функций sin() или cos(), хотя в других реализациях эти вызовы будут.

#include 
void iqDFT12( float *y, const short *x )
{
    for( int i=0; i<128; i++ )
    {
        const float wi = (float)i*(2.0f*3.1415927f/4096.0f);
        const float sii = sinf( wi );
        const float coi = cosf( wi );

        float co = 1.0f;
        float si = 0.0f;
        float acco = 0.0f;
        float acsi = 0.0f;
        for( int j=0; j<4096; j++ )
        {
            const float f = (float)(x[2*j+0]+x[2*j+1]);
            const float oco = co;
            acco += co*f; co = co*coi -  si*sii;
            acsi += si*f; si = si*coi + oco*sii;
        }
        y[i] = sqrtf(acco*acco+acsi*acsi)*(1.0f/32767.0f);
    }
}