Трюк с тригонометрией

Original author: Inigo Quilez
  • Translation

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


$\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 = \{\cos\alpha, \sin\alpha\}$ по окружности с шагом в $\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 <math.h>
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);
    }
}

Only registered users can participate in poll. Log in, please.

Альтернативное голосование

  • 90.3%+1253
  • 9.6%-127
Support the author
Share post

Similar posts

AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 28

    +7
    Статья выглядит так, как будто автор не знаком с комплексными числами.
      +1
      хотя в других реализациях эти вызовы будут

      Скорее всего не во внутреннем цикле — там просто будут использоваться табличные значения.
      Кстати насколько разошлись значения синуса и косинуса от реальных на последнем шаге?
        +3

        Оригинальная статья умалчивает, но я провел собственное расследование: ссылка
        Результаты:


        BEFORE: Sin = +0.00000000000000000000 (expected 0), Cos = +1.00000000000000000000 (expected 1)
        AFTER:  Sin = -0.00000000000000035196 (expected 0), Cos = +1.00000000000642160630 (expected 1)
        ERROR:  Sin = -0.00000000000000035196 (expected 0), Cos = +0.00000000000642160630 (expected 0)
        ERROR:  Sin = -3.51959e-16 (expected 0), Cos = +6.42161e-12 (expected 0)

        То есть, на 1 миллион итераций точность около 10-12 знаков.
        Для сравнения:
        Точность float (32): 7-8 знаков
        Точность double (64): 15-17 знаков
        Точность long double (80): 18-19 знаков
        Я использовал long double.

          +1
          От шага же зависит. Чем меньше шаг, тем меньше врёт.
          Для генератора сигналов на бюджетных МК — вполне годно. Там вытянуть амплитуду можно в аналоговой части
            0

            Там нет монотонной зависимости, сначала точность растёт, но всё же с определённого момента начинает падать.

              0
              Ох, не первых лабах по симуляции показывают как уменьшение дельты приводит к взрыву. А всё из-за флотов и ошибки в точности, при маленькой дельте, ошибка одной итерации начинает приближаться к самой дельте. Пробуйте.
                0
                да пробовал, еще на фиксированной точке для аудио нужд с частотой колебаний от 10Гц, до нескольких кГц. Т.е. дельта не сферическая, а вполне себе злободневная.
                Пользоваться, в принципе, можно. Еще и короче (и шустрее) библиотечных вариантов получается. С фиксированной точкой дружит.

                А, например, arm_sin_f32 считает синус по табличке с шагом чуть ли не в 1 градус и кубической интерполяции. К этой библиотеке тоже есть вопросы по точности.
                0

                Чем меньше шаг — тем больше самих шагов.

                0
                Гражданин Герцель вам 5 умножений во внутреннем цикле поможет сэкономить, с осциллятором второго порядка, вместо честного поворота умножением на матрицу:
                en.wikipedia.org/wiki/Goertzel_algorithm#Power-spectrum_terms
                и который вполне себе устойчивый и без накопления ошибок
                  0
                  Если посмотреть с дурой стороны, я бы назвал это не накоплением ошибок, а дискретизацией. На самом деле в самой волне, при генерации ее поворотом, ни каких ошибок и отклонений нету, есть только несоответствие заданной частоте и оно довольно мизерное, но не генерируемой.
                    0
                    Всё-таки Гёрцель.
                      0
                      Ссылаемся на Википепию? Хорошо…
                      В русскоязычной литературе нет устоявшегося варианта транскрипции фамилии автора алгоритма. Распространены варианты «Алгоритм Герцеля», «Алгоритм Гертцеля», «Алгоритм Горцеля» и другие.

                      Ваша же ссылка
                    +1
                    Это отличный результат.
                    Тогда не понятно, почему автор вызывает внутри внешнего цикла:
                    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 );

                    ведь по сути там такой-же трюк подходит
                      0

                      Эта часть отрабатывает в 4096 раз реже, так что экономия будет уже не столь существенна.

                  +2
                  Как бы, для кого — открытие, а для кого — серые будни, известные еще с 8-биток.
                  Добавить еще пару действий и получится резонансный фильтр 2го порядка, LP/HP/BP.
                  Имхо, 128х полосовыми фильтрами можно было бы менее ресурсоемко сделать подобие FFT, хватая по 1 сэмплу.
                    0
                    Быстрый FFT уже изобретен.
                      +1
                      быстрый быстрый?
                  • UFO just landed and posted this here
                      0
                      иногда синусы и косинусы даже табличные либо не могут обеспечить видимую точность (и такое бывает), иногда оин не доступны в компиляторе, не работаю потому что их забыли реализовать (спасибо менеджерам спешившим продать новый новую железку), и просто желание программиста без конкретных явных причин, только какие-то лично-субъективные
                      +3

                      Вот недавно на хабре была статья, как числа округлять. Сегодня: как применять школьную формулу тригонометрии. Это все сложно. Напишите пожалуйста, как извлекать квадратный корень числа с плавающей точкой в C++. Это будет очень актуально для хабра.

                        +2

                        Корень тоже был, см. 0x5F3759DF

                          0
                          1/на корень, если что
                            0

                            Таких констант можно кучу напридумывать для любого показателя степени от -1 до 1, о чем в статье «Магическая константа» 0x5f3759df и говорится.

                            0
                            И ещё вот (что кстати на PC-платформе не имеет смысла, потому что там есть FPU с 80-битной точностью для промежуточных вычислений).
                              0
                              если числа меньше 1, что характерно для DSP, я бы почитал чего есть быстрого.
                              +1
                              Посмотрите, как организован внутренний цикл, тот, что выполняется 4096 раз: ни одного вызова функций sin() или cos(), хотя в других реализациях эти вызовы будут.
                              В других реализациях этих вызовов тоже нет — там используются предварительно подсчитанные табличные данные. Более того — трюк с последовательным поворотом вектора для вычисления FFT хорошо известен и в частности используется в книге «Numerical Recipes» (страница 612 в третьем издании).
                                +3
                                Продолжайте писать. Не все начинали тут с восьмибиток и реайлтайма.

                                Статья полезная. Тут рядом где-то статья была, как человек графики для Телеграма оптимизиировал — думаю, ему бы понравилось.
                                  0
                                  Если тригонометрические функции по какой-то причине так неприятны, то почему бы не перейти к показательным функциям с помощью формулы Эйлера?

                                  Only users with full accounts can post comments. Log in, please.