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


$\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);
    }
}

Комментарии (27)


  1. Refridgerator
    06.08.2019 15:29
    -1

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


  1. GCU
    06.08.2019 15:39

    хотя в других реализациях эти вызовы будут

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


    1. developerxyz Автор
      06.08.2019 16:16
      +1

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


      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. IgorPie
        06.08.2019 16:18

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


        1. slopestyler
          07.08.2019 10:17

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


        1. Aquahawk
          07.08.2019 11:51

          Ох, не первых лабах по симуляции показывают как уменьшение дельты приводит к взрыву. А всё из-за флотов и ошибки в точности, при маленькой дельте, ошибка одной итерации начинает приближаться к самой дельте. Пробуйте.


          1. IgorPie
            07.08.2019 16:37

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

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


        1. mayorovp
          07.08.2019 13:15

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


      1. pvvv
        06.08.2019 19:54

        Гражданин Герцель вам 5 умножений во внутреннем цикле поможет сэкономить, с осциллятором второго порядка, вместо честного поворота умножением на матрицу:
        en.wikipedia.org/wiki/Goertzel_algorithm#Power-spectrum_terms
        и который вполне себе устойчивый и без накопления ошибок


        1. panteleymonov
          07.08.2019 00:47

          Если посмотреть с дурой стороны, я бы назвал это не накоплением ошибок, а дискретизацией. На самом деле в самой волне, при генерации ее поворотом, ни каких ошибок и отклонений нету, есть только несоответствие заданной частоте и оно довольно мизерное, но не генерируемой.


        1. BD9
          07.08.2019 17:06

          Всё-таки Гёрцель.


      1. GCU
        06.08.2019 20:01
        +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 );

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


        1. Source
          07.08.2019 19:22

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


  1. IgorPie
    06.08.2019 15:59

    Как бы, для кого — открытие, а для кого — серые будни, известные еще с 8-биток.
    Добавить еще пару действий и получится резонансный фильтр 2го порядка, LP/HP/BP.
    Имхо, 128х полосовыми фильтрами можно было бы менее ресурсоемко сделать подобие FFT, хватая по 1 сэмплу.


    1. panteleymonov
      07.08.2019 00:53

      Быстрый FFT уже изобретен.


      1. mad_god
        07.08.2019 10:54

        быстрый быстрый?


  1. MSC6502
    06.08.2019 21:49
    +1

    А что автор так боится синусов косинусов? Таблица значений на четверть периода и все. Скорость для 1024-точечного Фурье на 386 процессоре была достаточна для визуализации процесса виброиспытаний некоторого изделия в реальном масштабе времени.


    1. DarkWolf13
      07.08.2019 01:08

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


  1. azudem
    07.08.2019 06:31

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


    1. screwer
      07.08.2019 08:37
      +1

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


      1. IgorPie
        07.08.2019 16:38

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


      1. KvanTTT
        09.08.2019 01:56

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


    1. Refridgerator
      07.08.2019 08:46

      И ещё вот (что кстати на PC-платформе не имеет смысла, потому что там есть FPU с 80-битной точностью для промежуточных вычислений).


    1. IgorPie
      07.08.2019 16:38

      если числа меньше 1, что характерно для DSP, я бы почитал чего есть быстрого.


  1. Refridgerator
    07.08.2019 07:06
    +1

    Посмотрите, как организован внутренний цикл, тот, что выполняется 4096 раз: ни одного вызова функций sin() или cos(), хотя в других реализациях эти вызовы будут.
    В других реализациях этих вызовов тоже нет — там используются предварительно подсчитанные табличные данные. Более того — трюк с последовательным поворотом вектора для вычисления FFT хорошо известен и в частности используется в книге «Numerical Recipes» (страница 612 в третьем издании).


  1. Zoolander
    07.08.2019 08:32
    +2

    Продолжайте писать. Не все начинали тут с восьмибиток и реайлтайма.

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


  1. Kron7
    07.08.2019 11:26

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