Comments 93
PS Если не гнаться за точностью, то есть элементарный Alpha max plus beta min algorithm.
Вы сами изобрели CORDIC или просто попытались изложить его для непосвященного слушателя?
Изобрел сам. Ну и в публикации об этом есть. Ну и реализацию выложил, хотя уже еще ее оптимизировал.
Это верно для соответствия 1:1 при расчете классическим «пифагоровым» методом.
При меньших таблицах — да тоже работает, но точность в некоторых значениях начинает убывать.
… пришел к простому империческому правилу.Лучше писать эмпирическому, конечно если Вы не имели в виду имперскому.
Извините, функция Ctrl/Enter (чтобы послать сообщение об опечатке лично автору) в комментариях не работает.
Ну, строго говоряя, этому изобретению тоже сто лет в обед.
Я не на что не претендую. Вы, слишком прямо ко всему отнеслись. Я описал свой опыт, который выглядит как история, под названием «как я отказался...».
Под изобретением я имею ввиду не что-то фундаментальное, но способность собрать новый дом из старых кирпичей.
Например, человек, который придумал автомагнитолу, всего-лишь придумал, как вставить радио в автомобиль, хотя и радио и автомобиль изобрели до него. И он — изобретатель. А есть человек, который придумал прикрепить ластик на обратную сторону карандаша, и это тоже изобретение.
Про такие вещи рассказывают на какой дисциплине? Вычмат?
Э… Мне казалось, что LUT каждый, кто пишет обработку сигналов, должен левой пяткой делать. Конкретно в вашем подходе – вы даже чуток усложнили выкладки, тригонометрия вообще ни к чему, переход к косинусу и арктангенсу не требовался: просто переходим от функции двух переменных sqrt(y²+x²) к произведению двух функций от одной переменной x*sqrt(1+t²), где t=y/x, и LUT у нас в кармане. Не уверен, что это оптимальное решение в вашем случае, но рабочее, факт.
И личная просьба: прогоните текст через спеллчекер. Хотя бы s/длинна/длина/g. Буду очень благода...
А вы пробовали вычислять квадратный корень через разложение в ряд Тейлора ?
Наверняка не скажу, но навскидку там получается примерно вдвое меньше инструкций + не нужен LUT.
Кроме того, вполне допускаю, что есть целочисленные варианты алгоритма. Сам метод Ньютона вполне себе целочисленный. Нужно будет только грамотно прикинуть первое приближение. На худой конец использовать для него тот же LUT.
И ещё момент. Если всё это применялось для ускорения БПФ, не уж-то получилось быстрее, чем предложенная самой STM реализация для cortex-ядер?
Нужно будет только грамотно прикинуть первое приближение.Обычно применяется подсчет ведущих нулей.
И ещё момент. Если всё это применялось для ускорения БПФ, не уж-то получилось быстрее, чем предложенная самой STM реализация для cortex-ядер?У меня это применялось не для БПФ, но хочу попробовать.
В ядре «быстрого» обратного корня находится хак, связанный с бинарным представлением floating-point т.е. алгоритму нужно fp-представление, причем в ieee754 формате
А лучше ли это, чем быстрый инверсный квадратный корень?
Мне этот метод не понравился как «апроксимативный». Да, это неплохое решение, но, как уже написали ниже, не для 107й STMки.
… + не нужен LUT.
Про LUT тоже много думал, в основном, как генерировать. Чтобы можно было развернуть в ОЗУ.
Третьим будете (кто эту ошибку заметил). Кстати, Вы про сообщение об ошибке автору через выделение+Ctrl-Enter в курсе?
Заметил, кстати, что в последнее время ну о-о-очень многие этим грешат («длинной», имею в виду).
Заметил, кстати, что в последнее время ну о-о-очень многие этим грешат («длинной», имею в виду).
Лично у меня это было я знаю почему. Одно время переписывался с коллегой, который все время писал «антена». Мне это настолько резало глаз, что прямо коробило. Поправить коллегу я не решался, он был заметно старше и опытней в теме.
Но в моей голове возникла гиперкомпенсация по поводу двойного «н».
многие этим грешат («длинной», имею в виду)
Я, например, так слышу. Знаю, как пишется, но если бы писал на слух — писал бы "длинна".
А почему не запилить таблицей? Если это STM32 то памяти должно хватить
Между точками таблицы линейная или квадратичная интерполяция
На порядок быстрее бы работало.
Похоже на дихотомический алгоритм…Да, половина. Как раз чтобы LUТ уменьшить.
А почему не запилить таблицей?Если все запилить таблицей, то надо 32Мб. Суть алгоритма в небольшом LUT.
Между точками таблицы линейная или квадратичная интерполяцияМежду точками нет интерполяции. Просто количество точек больше необходимого.
При всём уважении, из частного Вы брали арктангенс, а не «арккосинус». Извините, бросается в глаза
Здорово. Но было бы круто ещё бенчмарки посмотреть. Любопытства ради
A = SQR(X^2+Y^2) ≈ Max(X, Y) + Min(X, Y)/2, другими словами задача сводится к определению максимального из них, сдвигу и сумме — да, с погрешностью, зато халява.
Подобрав другие коэффициенты при Max и Min можно поиграться с точностью результата. Под спойлером приводится зависимость погрешности от угла для различных коэффициентов
кажется 20 точек мне вполне хватало для поворотов.
В целом я смотрю, когда (сейчас) кому-то нужен алгоритм, то его передирают начисто с какого-нибудь сайта и запиливают не разбираясь в сути и принципе. Работает? Да! Что еще нужно. Тормозит? Возьмем процессор по мощнее!
Пока Кули и Тьюки не придумали БПФ, все использовали ДПФ. И т. д. :)
;r5 = X
;r6 = Y
;Magnitude
r2 = 0x7FF
r3 = 0x7FF
сравнить r5 и r2, результат записывается в регистр флагов
r5 -= r2 при условии, что установлен флаг "больше"
r2 -= r5 при условии, что установлен флаг "меньше или равно"
r5 = r2 при условии, что установлен флаг "меньше или равно"
;Magnitude_X ;X in r5
сравнить r6 и r3, результат записывается в регистр флагов
r6 -= r3 при условии, что установлен флаг "больше"
r3 -= r6 при условии, что установлен флаг "меньше или равно"
r6 = r3 при условии, что установлен флаг "меньше или равно"
;Magnitude_Y ;Y in r6
и так далее.
Справедливости ради, кусок выдранный из общей прошивки рабочего изделия я выдрал и выложил сюда. Для наглядности и понимания сути всего алгоритма в целом.
Например, в начале указаны значения r2 и r3 как 0x7FF. В реальной прошивке они не указываются явно, а являются входными парамертрами, которые формирует другая часть адаптивной корректировки по медиане.
Предложенный мной алгоритм лишь рыба, которая создана для понимания.
Оптимизировать можно много чего, но надо ориентироваться на задачу, которую решаем.
может дешевле более мощный проц взять?
Сколько времени потратили и какой тираж применения?
С другой стороны, за цену процессора вы платите один раз и используете тоже один раз, а за цену НИОКР платите один раз, а используете многократно...
может дешевле более мощный проц взять?Проц достаточно производителен, просто я люблю все оптимизировать. Иногда могу заморочиться, что программа работает быстрее, если переставить местами несколько инструкций.
Сколько времени потратили и какой тираж применения?Сверх обычного, потратил сутки-полторы не более. Тираж, сейчас опытная партия изделий (после карантина) 10 шт. Потом испытания и сертификация. Затем первичная поставка на 100+ изделий по доп. соглашению. А там посмотрим.
Суть в заделе на будущее, такой алгоритм много где пригодится. Кроме того, хоть я и пишу на АСМе, я уже не пишу много с нуля, а копипастю из собственных проектов, или напрямую или макробиблиотеками.
А если серьёзно, спасибо за содержательную статью.
Навряд ли на столько. Всё-таки людей, которые могут в одного разработать ASIC, во много раз меньше людей, собравших метеостанцию на ардуине, и стоят они пропорционально встречаемости, так что такая разработка отбивалась бы на стотысячных тиражах. В дешёвые и популярные FPGA я скорее поверю.
DARPA — это такой дедушка-инкубатор стартапов, там "выстреливает" хорошо если процентов 10, да и отдельные заявки могут висеть десятилетиями, например, подводный лидар.
С другой стороны очевидно, что запрос на такой софт есть, раз на это готовы потратить почти 200 миллионов долларов на исследования и разработку, а значит рано или поздно что-то появится. Возможно не этот проект, возможно это будет платная среда.
Я бы предположил, что компании-производители чипов очень заинтересованы в уменьшении влияния дефицита специалистов и экспертов на процесс и в увеличении пропускной способности, потому что это очень сильно тормозит наращивание оборотов их бизнеса. Так что с моей точки зрения это только вопрос времени. Как и с любыми САПРами.
Обещают превратить разработку ASIC из HDL в процесс подобный компиляции программы из исходников.Как говорит Топа из utopia show, «хочется верить, но верится с трудом» (с).
У нас на работе заказывали ASIC, так это какие-то неприличные деньги и полгода какого-то промежуточного перепроектирования.
Или можно так: sqrt(x^2+y^2)~=x+a*y^2+b*y^4; вроде, точность д.б. нормальной.
sin(float x)
работал на 2-3% быстрее, чем arr[(int) float x]
. У меня возникло подозрение, что вычисление синуса как-то хитро оптимизировано (набором инструкций SSE, например), что даже прямое обращение в массив не смогло выиграть. Эту тему я дальше копать не стал, она так и осталась для меня тайной. Может, какой-то флаг при компиляции был отключен (target был Release в Borland Developer Studio C++). Это было в 2010-м году на Intel Core 2 Duo. Если кто знает, в чём секрет, поделитесь =)Мне тоже приходилось делать БПФ...
… я пробовал вычисление синуса заменить на таблицу...
Когда я писал на x86, то я использовал FPU.
Но БПФ лучше писать «бабочковую» версию, там есть небольшой геморрой с поворотными коэффициентами, но в целом алгоритм хороший. А еще лучше, если количество точек кратно четырем, то делать 4-кратную бабочку, там еще больше выигрыш.
… делать БПФ (на миллионе значений)...Зачем так много? Для какой цели?
Зачем так много? Для какой цели?Корреляционный анализ высокочастотных сигналов. Два сигнала частотой 500КГц длительностью 1,05 секунд, разрядность 16 бит, если не ошибаюсь. Длительность брал, чтобы было степенью двойки, алгоритм БПФ, с бабочкой. Я ради интереса пробовал несколько реализаций, ковырялся, и на каком-то этапе проверял один и тот же алгоритм, заменяя в нём только вычисление синуса на массив — и выигрыш синуса стал для меня главным открытием.
Вот тут умные люди пишут, что
Lookup tables never really work on modern processors — random access to memory is a significant bottleneck these days. The only real way to improve this is to use vector processing (SIMD / SSE) but that does depend on the algorithm being used.То есть вычисление синуса оптимизировано и оно ожидаемо быстрее, чем простой прямой доступ к памяти (в отличие от оптимизированного доступа через векторы SIMD, SSE). Проверять, конечно, не буду, поверю им на слово.
Младшее значение 0x800h, что означает 2048. Cтаршее 0xB50, и оно означает 2896.
Старшие 4 бита не используются. 12...15 биты всегда 0.
Бит 11 всегда единица.
Значения в пределе от 2048+(0...848).
Используя вместо 16бит 8бит можно получить таблицу в 2 раза меньше.
Поиск байта, умножить на три сдвигом и/или сложением + 2048
На сколько упадет точность при входных целых числах?
На самом деле 848/3=283, что более чем 256, значит надо делить на 4, а это 212 значений.
Точность падает, до расхождения в 3 целых МЗР.
Правильнее в алгоритме (пункт 5) умножать Y не на 2048, а на 256. Тогда после деления Y/X
получаемый адрес будет состоять из 8 бит. И можно уже использовать «байтную таблицу», она получается уже 256 байт, вместо 4кБ.
Но точность будет на 2 МЗР отличаться.
Вполне можно сделать, но в моей задаче, нужна была точность до целой единицы.
Ну и 4кБ даже для AVR не много, тем более для STM32. :)
Если Вы посмотрите на таблицу в моей статье, то увидите, что в начале таблицы много одинаковых значений, потом, почти в два раза меньше, значений на единицу больше и т.д. А к концу таблицы значения идут подряд.
Да, таблицу можно ужимать без потери, но это вызовет замедление алгоритма.
Как здесь писали, есть компромиссы между объемом данных и скоростью. Я описал случай с максимальной скоростью, при среднем объеме таблицы.
Не ясно, с какой скоростью меняется сигнал на входе. Скорость АЦП. Допустимое запаздывание по обработке. Точность результата и его размерность.
Я вам сразу указал, что можно сократить таблицу на четверть за счет не используемых битов.
Не сталкивался с данным процессором, можно ли там использовать адресацию со смещение по битам L MW[AR1,P#0.0]
Я вам сразу указал, что можно сократить таблицу на четверть за счет не используемых битов.Нет, такой адресации нет.
Не сталкивался с данным процессором, можно ли там использовать адресацию со смещение по битам L MW[AR1,P#0.0]
Да, я знаю, что LUT можно сократить. Я знаю примерно 4,5 способа сделать так, чтобы при этом не терялась точность. Я знаю, как ускорить алгоритм.
Дело в другом.
Дело не в конкретной реализации, это «рыба», понимаете?
Дело в том, чтобы показать, что решить задачу можно иначе, чем «принято» это делать. И алгоритм здесь — пример такого решения.
То, что Вы предлагаете, не меняет алгоритма, а меняет суть его реализации в программе.
Проведите эксперимент — попросите любого знакомого программиста решить задачу нахождения длины вектора. И Вы увидите, Вам решат через корень. Машинально. Конечно при условии, что мой пост не читали.
Зачем так все усложнять?Кто усложняет? По моему все максимально упрощено.
Если вернуться к началу задачи то есть некая точка к координатами X Y, причем модуль значения координаты помещается в 11 бит.Верно. Обратите внимание, что модуль берется в первом и втором действии приведенного алгоритма.
Итого максимальная длина отрезка 2048...Нет. Максимальная длина 2048*sqrt(2) примерно = 2895. А 2048 — это максимальная длина проекции на ось.
… и угол ( если отбросить знак) от 0 до 90 градусов.Верно. Однако 3 действием я свел угол к 45.
Так как длина отрезка величина дискретная, то угол тоже дискретная величина.Да. И длина вектора и угол величины дискретные. Но нет, дискретность угла не вытекает из дискретности длины. Скорее и то, и другое есть следствия одной причины.
Количество значений, которое может принять угол, равно длине четверти окружности в «пикселях». Для отрезка длиной 2048 получаем 3217 значений. Приблизительно 0,028 градуса каждое значение. Погрешность «пол деления» итого 0,042 градуса.Допустим.
Погрешность АЦП пусть даже один разряд.Погрешность АЦП не нужно рассматривать. Вектор может быть получен усреднением множества точек с АЦП, или проходить через какой-нибудь фильтр.
Итого 0,05 градуса или менее 0,1%. Но если рассматривать отрезок длиной, скажем 205 ( 10% от максимального уровня) то тогда будет всего 322 значения. Погрешность все равно менее 1%. Далее выбираем допустимую погрешность и вычисляем допустимую дискретность угла. По ней вычисляем «длину отрезка». Имея длину отрезка можно получить его координаты при условии что он совпадает с тем, угол которого надо измерить. А дальше имея координаты конечной точки это отрезка по табличке вычисляем угол.Вы сейчас описали обратный метод разработки. Примерно все это я проделал на этапе сбора сведений (еще до написания алгоритма) для решения этой задачи.
Более того, я просчитал все возможные варианты. Ну то есть буквально — все погрешности для всех X и всех Y при всех длинах таблиц. Написал программу, которая мне «пробрутфорсила» это все. И сделал вывод, как я и писал, импирический, что размер LUT должен быть равен максимальной длине проекции на ось. То есть для 8-битных векторов нужна таблица 256 значений. И т.д.
И просчитывать все возможные значения с учетом разрядности АЦП не пришлось бы.
Очень грубый пример — предположим в зависимости от угла вектора надо получить амплитудно модулируемый сигнал.Что значит «в зависимости от угла вектора»? Как угол влияет на амплитуду? Угол — это фаза.
Если исполнительное устройство оперирует в процентах и минимально различимое изменение 1%, то дискретность угла в 3.6 градуса уже достаточна.Да, но 1% при амплитудной модуляции это около 37 дБ. Это очень маленький динам. диапазон.
Я указал способ при котором есть алгоритм, который дает результат такой-же как классический метод (тот что через корень).
Зачем мне «ухудьшать» точность, подгоняя ее под какие-то «определенные характеристики»?
И использование 12 разрядного ЦАП не имеет никакого смысла в этом случае.В той задаче, которую я решал 12 разрядный ЦАП был выбран не с проста. И предварительные расчеты, и НИР и ОКР показали необходимость именно 12 бит. 10-ти битного было мало.
Но, причем здесь алгоритм?
Моя статья об алгоритме, а не о изделии. Можете обратить внимание, что кроме однократного упоминания ЦАП 12-бит, там нет больше ни слова. И то, про 12 бит говориться только в контексте, для объяснения причины использования во входных данных чисел без знака и виртуального нуля.
… Только усложнение конструкции, алгоритма и удорожание компонентной базы. Прежде всего надо понять диапазон выходных значений и необходимую точность их. Все остальное потом.Так и было сделано на этапе НИР. Но Вы, почему-то, говорите о недостаточной продуманности самого изделия, несмотря на то, что не знаете, что за изделие, какие параметры у него должны быть, что вышло в итоге, о коком-то удорожании (непонятно по сравнению с чем).
Еще раз, речь в статье не про изделие, а про алгоритм, который поможет заменить метод.
Вы привели алгоритм вычисления угла вектора ( вектора у Вас на картинках)...Нет, я привел алгоритм вычисления длины вектора, а не угла.
я предложил рассмотреть случай, что значение этого угла используется для амплитудной модуляции — в зависимости от значения угла меняется амплитуда некоего сигнала, просто для примера.Хорошо, но опять-же в моем алгоритме угол как раз выкидывается, его величина не важна. Если бы нужен был бы угол, был бы другой алгоритм.
37дБ, по известной формуле, это изменение уровня в 70.79 раза 1% это изменение уровня в 0,99 раз ( за исключением изменения уровня от 0).Хорошо, 37дБ я написал, потому-что прикинул в уме. И навскидку ошибся на 3 дБ. Давайте укажем 39.96688дБ. :)
1% изменения уровня это всего 0,0864 дБ… «Разы» = ( x+ x/100*1)/x = 1,01 в децибелах всего 0,0864
Вы писали:
Если исполнительное устройство оперирует в процентах и минимально различимое изменение 1%… то это означает, что «устройство» может видеть всего 100 значений, а значит его динам. диапазон около 40 дБ.
При всем уважении, по-моему Вы ударились в демагогию.
Что именно Вы хотите мне сказать?
Как я отказался от вычисления квадратного корня