Скорее всего, вам известны следующие соотношения еще со школы:
Когда вы в детстве впервые познакомились с этой формулой, скорее всего, вашим первым чувством была боль из-за того, что эту формулу надо запомнить. Это очень плохо, потому что на самом деле вам не нужно запоминать эту формулу — она сама выводится, когда вы поворачиваете треугольник на бумаге. На самом деле, я делаю то же самое, когда записываю эту формулу. Это толкование будет очевидным к середине этой статьи. Но сейчас, чтобы оставить все веселье на потом и отодвинуть момент, когда вы скажете "Эврика!", давайте подумаем, а зачем нам вообще задумываться об этой формуле.
Вступление
Тригонометрические функции 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
и каков ее диапазон значений (в примере выше — ). Единственное, что нас беспокоит, так это то, что есть цикл, который постоянно вызывает sin()
и cos()
с параметром, который увеличивается в постоянных шагах (в данном случае, ). Эта статья о том, как оптимизировать этот код для скорости таким образом, что одни и те же вычисления могут выполняться вообще без использования функций sin()
или cos()
(во внутреннем цикле), и даже более быстрой функции sincos()
.
Но если посмотреть на первую формулу в статье, мы увидим, что если и , мы можем переписать это как
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;
}
Толкование
К настоящему моменту мы слепо играли с математикой не понимая, что же на самом деле происходит. Давайте перепишем внутренний цикл так:
Некоторые из вас могли заметить, что это — формула поворота объекта в двухмерном пространстве. Если вы все еще не поняли этого, возможно, матричная форма вам поможет.
В самом деле, sin(t)
и cos(t)
можно сгруппировать в вектор длины 1 и отрисовать как точку на экране. Назовем этот вектор x
. Тогда, . Значит, векторная форма выражения — , где . Мы видим, что наш цикл выполняет небольшой двухмерный поворот каждую итерацию так, что x
вращается по кругу во время выполнения цикла. Каждую итерацию x
вращается на радиан.
Итак, в основном,
это формула движения точки по окружности с шагом в радиан. Чтобы это сделать, мы построим одну из двух осей поворота, к примеру, . Первый компонент поворота — проекция на . Так как и нормализованы (имеют длину 1), проекция — их скалярное произведение. Следовательно, , и конечно второй компонент — антипроекция, которую можно найти, спроецировав на перпендикулярную ось, . Мы можем создать этот вектор, развернув координаты и изменить знак на противоположный у первой координаты:
Примечания
Обычно вы должны иметь возможность выполнять эти крошечные вращения снова и снова. В самом деле, , что означает, что матрица не увеличивает и не уменьшает пространство, к которому применена, что значит, что будет двигаться по идеальной окружности. Однако из-за того, что компьютеры не точны, будет двигаться по спирали и в конце концов совпадет с центром окружности вращения. У меня не возникало таких проблем, но я думаю, что они могут возникнуть при очень больших num
, т.е. маленьких углах поворота.
Пример
В Kindercrasher, 4096-байтной демке из 2008 (скриншот на КДПВ), группа сфер пульсирует под музыку. Для этого я сосчитал преобразование Фурье звука. Существуют алгоритмы, делающие это в реальном времени, к примеру, FFT. Однако, мой код должен вместиться в несколько килобайт, и я решил пойти иным путем. Вместо реализации FFT, я написал DFT по его простому определению. Вы можете проверить это в википедии.
Моя функция также принимает 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)
GCU
06.08.2019 15:39хотя в других реализациях эти вызовы будут
Скорее всего не во внутреннем цикле — там просто будут использоваться табличные значения.
Кстати насколько разошлись значения синуса и косинуса от реальных на последнем шаге?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.IgorPie
06.08.2019 16:18От шага же зависит. Чем меньше шаг, тем меньше врёт.
Для генератора сигналов на бюджетных МК — вполне годно. Там вытянуть амплитуду можно в аналоговой частиslopestyler
07.08.2019 10:17Там нет монотонной зависимости, сначала точность растёт, но всё же с определённого момента начинает падать.
Aquahawk
07.08.2019 11:51Ох, не первых лабах по симуляции показывают как уменьшение дельты приводит к взрыву. А всё из-за флотов и ошибки в точности, при маленькой дельте, ошибка одной итерации начинает приближаться к самой дельте. Пробуйте.
IgorPie
07.08.2019 16:37да пробовал, еще на фиксированной точке для аудио нужд с частотой колебаний от 10Гц, до нескольких кГц. Т.е. дельта не сферическая, а вполне себе злободневная.
Пользоваться, в принципе, можно. Еще и короче (и шустрее) библиотечных вариантов получается. С фиксированной точкой дружит.
А, например, arm_sin_f32 считает синус по табличке с шагом чуть ли не в 1 градус и кубической интерполяции. К этой библиотеке тоже есть вопросы по точности.
pvvv
06.08.2019 19:54Гражданин Герцель вам 5 умножений во внутреннем цикле поможет сэкономить, с осциллятором второго порядка, вместо честного поворота умножением на матрицу:
en.wikipedia.org/wiki/Goertzel_algorithm#Power-spectrum_terms
и который вполне себе устойчивый и без накопления ошибокpanteleymonov
07.08.2019 00:47Если посмотреть с дурой стороны, я бы назвал это не накоплением ошибок, а дискретизацией. На самом деле в самой волне, при генерации ее поворотом, ни каких ошибок и отклонений нету, есть только несоответствие заданной частоте и оно довольно мизерное, но не генерируемой.
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 );
ведь по сути там такой-же трюк подходитSource
07.08.2019 19:22Эта часть отрабатывает в 4096 раз реже, так что экономия будет уже не столь существенна.
IgorPie
06.08.2019 15:59Как бы, для кого — открытие, а для кого — серые будни, известные еще с 8-биток.
Добавить еще пару действий и получится резонансный фильтр 2го порядка, LP/HP/BP.
Имхо, 128х полосовыми фильтрами можно было бы менее ресурсоемко сделать подобие FFT, хватая по 1 сэмплу.
MSC6502
06.08.2019 21:49+1А что автор так боится синусов косинусов? Таблица значений на четверть периода и все. Скорость для 1024-точечного Фурье на 386 процессоре была достаточна для визуализации процесса виброиспытаний некоторого изделия в реальном масштабе времени.
DarkWolf13
07.08.2019 01:08иногда синусы и косинусы даже табличные либо не могут обеспечить видимую точность (и такое бывает), иногда оин не доступны в компиляторе, не работаю потому что их забыли реализовать (спасибо менеджерам спешившим продать новый новую железку), и просто желание программиста без конкретных явных причин, только какие-то лично-субъективные
azudem
07.08.2019 06:31Вот недавно на хабре была статья, как числа округлять. Сегодня: как применять школьную формулу тригонометрии. Это все сложно. Напишите пожалуйста, как извлекать квадратный корень числа с плавающей точкой в C++. Это будет очень актуально для хабра.
screwer
07.08.2019 08:37+1Корень тоже был, см. 0x5F3759DF
KvanTTT
09.08.2019 01:56Таких констант можно кучу напридумывать для любого показателя степени от -1 до 1, о чем в статье «Магическая константа» 0x5f3759df и говорится.
Refridgerator
07.08.2019 08:46И ещё вот (что кстати на PC-платформе не имеет смысла, потому что там есть FPU с 80-битной точностью для промежуточных вычислений).
IgorPie
07.08.2019 16:38если числа меньше 1, что характерно для DSP, я бы почитал чего есть быстрого.
Refridgerator
07.08.2019 07:06+1Посмотрите, как организован внутренний цикл, тот, что выполняется 4096 раз: ни одного вызова функций sin() или cos(), хотя в других реализациях эти вызовы будут.
В других реализациях этих вызовов тоже нет — там используются предварительно подсчитанные табличные данные. Более того — трюк с последовательным поворотом вектора для вычисления FFT хорошо известен и в частности используется в книге «Numerical Recipes» (страница 612 в третьем издании).
Zoolander
07.08.2019 08:32+2Продолжайте писать. Не все начинали тут с восьмибиток и реайлтайма.
Статья полезная. Тут рядом где-то статья была, как человек графики для Телеграма оптимизиировал — думаю, ему бы понравилось.
Kron7
07.08.2019 11:26Если тригонометрические функции по какой-то причине так неприятны, то почему бы не перейти к показательным функциям с помощью формулы Эйлера?
Refridgerator
Статья выглядит так, как будто автор не знаком с комплексными числами.