Введение

Всем привет, я работаю программистом DSP контроллеров в фирме, которая занимается разработкой и производством радаров. В начале моего пути, при решении задач по интеграции алгоритмов ЦОС, приходилось по долгу искать пути оптимизации вычислений.

То, что я собираюсь здесь изложить скорее всего будет интересно новичкам. Многие способы давным давно известны и реализованы в библиотеках. Однако, в свое время я потратил много часов для поиска разрозненной информации и понимания всех тонкостей.

Но довольно слов, давайте к делу.

Углы

Операций с углами в ЦОС много. Все начинается с нахождения фазы комплексного числа, далее их нужно складывать, вычитать, усреднять, находить дисперсию и т.д. При этом, нужно постоянно пользоваться функцией которая отсекает лишние периоды, ведь тригонометрический круг описывает только 360°.

Нормировка

Решение, которое решает практически все проблемы угловой математики и избавляет от кучи ненужных проверок, было подсмотрено мной в библиотеке для одного из западных DSP чипов.

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

Итак, угол хранится в типе uint8, где 0 это 0°, а 255 это почти 360°. Таким образом, 90° это 64, 180° это 128, а 270° это 192.

В чем же профит?
После сложения или вычитания больше не нужно отбрасывать лишние периоды и дергать функцию wrapTo360 (привет пользователям MatLab). При переполнении все отбросится само.

Например, если мы к 270° прибавим 180°, то во флотовой арифметике получим 450°, после нужно проверить больше ли результат 360° и если да, то вычесть лишние 360°. В результате получится 90°.

Теперь посчитаем тоже самое с нашей нормировкой: сложим 192 и 128, из-за переполнения uint8 сразу получим 64, что как раз-таки и соответствует 90°.

Благодаря свойствам дополнительного кода, приятным бонусом будет быстрый перевод в человеческие градусы или радианы в конце работы ЦОС алгоритма, при этом, изменение нотации на (-180°..+180°) никак не повлияет на производительность, просто одна ассемблерная инструкция поменяется на другую:

#define TO_ANGLE_COEF (360.0/256)

uint8 angle = 192; //270°

//0..360

float angle_deg_360 = angle * TO_ANGLE_COEF; //270

//-180..+180

float angle_deg_180 = (sint8) angle * TO_ANGLE_COEF; //-90

Очевидно, что единственным минусом данного подхода будет ошибка, которая зависит от выбранной нормировки. Например, в случае нормирования на 256 значений uint8 дискретность будет 360/256 = 1,4 градуса. Но если выбрать двухбайтный тип, то дискретность станет всего 360/65536 = 0.00549 градуса.

Усреднение

Усреднение углов в ЦОС используется для снижения влияния шумов.

Допустим, мы хотим усреднить 90° и 100°. (90 + 100)/2 = 95. То есть, 95° это середина дуги, между 90° и 100°.

Особенность усреднения углов заключается в том, что когда они "дребезжат" вокруг точки разрыва, классическая формула дает результат, в прямом смысле, противоположный. Например, результат усреднения двух близких углов 3° и 355°: (3°+355°)/2 = 179°. Но по смыслу понятно, что результат должен быть 359. То есть среднее двух углов это не просто середина дуги, а середина кратчайшей дуги.

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

#define P_A(s1, s2) (uint8)((s1) - (sint8)((s1) - (s2))/2) //average two phases

Выглядит не понятно, но присмотритесь: s1 - s2 это длина дуги. Мы не знаем какая это дуга - длинная или короткая, но нам нужна именно короткая. Какая дуга получится - зависит от направления вычитания.

Длинная дуга, это дуга больше 180°, то есть больше половины длины тригонометрической окружности. Когда такая длинна дуги преобразуется в знаковый тип, то значение >= 180° (напоминаю, в предложенной нормировке >= 128) станет отрицательным.

Например, s1 = 271°, s2 = 1°, разница будет 270°, что больше 180°. При преобразовании к знаковому типу она превратится в -90°, поделится пополам -45° и отнимется от s1. Так как вычитаемое отрицательное, то в результате получится сложение и результат будет 271°+45° = 316°, что корректно, в отличии от результата классической формулы.

Если же в результате разница будет меньше 180°, то эта формула также корректно сработает. Поменяем s1 и s2 местами, теперь s1 = 1°, s2 = 271°, s1 - s2 получится 90° (у нас же автоматический wrap, еще не забыли?). Преобразование к знаковому типу ничего не поменяет, и теперь 90°/2 отнимется от 1°, с автоматическим wrap'ом получим те же 316°.

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

Дисперсия

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

Для вычисления дисперсии, в первую очередь, необходимо среднее значение, его рассмотрели на предыдущем шаге.

В формуле дисперсии для углов, основная особенность такая же как и в вычислении среднего. Нужно находить разницу между средним и текущим, то есть найти длину дуги, а точнее кратчайшей дуги. Затем возводить в квадрат, накапливать сумму, поделить на n (или n-1, смотря что вам нужно) и затем извлечь корень. Здесь я опишу только поиск кратчайшей дуги, остальное - тривиальная задача:

(uint8)abs((sint8)(ave - sn))

Как и усреднение, такая конструкция не будет ошибаться в точках разрыва. Так что пользуйтесь.

Заключение

Данные приемы значительно ускоряют работу встраиваемых систем, однако могут оказаться полезными и на ПК в приложениях где считается большое количество углов. Поднятие точности до 4 байтного типа удовлетворит любой задаче.

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

Всем спасибо за внимание.