Введение

Всем привет, я работаю программистом 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 байтного типа удовлетворит любой задаче.

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

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

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


  1. alan008
    27.05.2022 00:18
    +5

    Стоило начать с того, что такое ЦОС.


    1. AWE64
      27.05.2022 00:57
      +3

      Цифровая обработка сигналов, надо полагать.


      1. alan008
        27.05.2022 11:10
        +2

        А как до этого догадаться тем, кто этим никогда не занимался :-)


        1. SmileOn Автор
          27.05.2022 12:32
          +3

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


          1. byman
            27.05.2022 13:15
            +1

            хорошая статья. Ждем следующую часть :)


      1. omxela
        27.05.2022 22:20

        малопонятное чтиво

        Чтиво как раз понятное, но ЦОС дешифруется неоднозначно, например, Цифровая образовательная среда (как сказал мне Гугл). И я долго думал. На самом деле это очень распространенный глюк в статьях: автор считает, что читатель сидит внутри его головы. Читатель же, видя такое, может подумать, что его не уважают (и это именно так, имхо). Ну чего стоило добавить расшифровочку в самом начале статьи? А если бы автор писал статью в настоящий журнал, его бы редактор вежливо попросил убрать аббревиатуру из заголовка, а в тексте - расшифровать.


        1. byman
          28.05.2022 10:29
          +2

          Ну чего стоило добавить расшифровочку в самом начале статьи?

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


    1. TheCalligrapher
      28.05.2022 12:19

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

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


  1. serjeant
    27.05.2022 09:13

    Спасибо за статью! А как называется фирма, если не секрет?


    1. SmileOn Автор
      27.05.2022 12:33

      Пожалуйста, надеюсь, статья кому-нибудь поможет. Фирма называется Ольвия.


  1. FGV
    27.05.2022 10:03

    О, да это же "моя любимая" фиксированная запятая. В свое время работал с 32 битным значением с ценой мл. разряда 360/2^32 ну или PI/2^31 (кому как больше нравится).

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


  1. LittleAlien
    27.05.2022 12:39

    Подозреваю, что подобных трюков много в статьях по старым псевдо-3D играм, наверняка в Doom 1,2 именно так углы и хранились. И таблицы для тригонометрии, разумеется.


  1. AndrewSu
    27.05.2022 20:16
    +2

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

    Возможно, вам пригодится circular mean.


    1. iShrimp
      28.05.2022 10:47
      +1

      Спасибо за ссылку! В Википедии есть и более обширный материал по круговой статистике: https://en.wikipedia.org/wiki/Directional_statistics


  1. AVDerov
    28.05.2022 10:57
    +2

    Спасибо за статью! Обычно работаешь во float и такие простые алгоритмы проходят мимо. Статья хорошо написана, кратко. Не все статьи должны разжевывать всё для всех.


  1. iShrimp
    28.05.2022 11:00
    +1

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

    Осторожнее со взятием абсолютного значения знакового числа. Результат abs() может внезапно оказаться отрицательным, если аргумент равен INT_MIN.

    Хотя, если разность не хранится, а сразу возводится в квадрат, это не важно. Но тогда и abs() не нужен.