Комментарии 33
Рекомендую юзать Cordic. Просто и надёжно.
Вы имеете ввиду готовую библиотеку, или алгоритм?
Много лет назад, когда работал в АО "Концерн "Созвездие", мы использовали данный алгоритм на FPGA, для синтеза синуса. Не моя часть, но я тоже ее касался, косвенно. Это были еще времена AT91RM9200, IXP465 и PXA270.
Здесь больше речь об оптимизации вычислений, т.к. при разложении в ряд тейлора, использовал все 32 бита, как дробную часть числа. Т.е. двоичная точка стоит прямо перед числом. А результат: 32-битное число. По этому после умножения, можно использовать только старшую часть результата, без сдвига. Но, потом, подумал, и сложил еще и младшую
Какое преимущество в данном случае даст Cordic?
Не сразу понял о чём идёт речь. Думал что используется метод неподвижной точки f(x) = x. И удивился, разве можно найти синус и косинус таким образом? Корень то легко, а вот тригономертические функции как!? Потом понял, что речь идёт про формат записи числа.
Можно подробнее про корень? Я не занимался ранее этим вопросом. Использовал готовые бибилиотеки. Отлаженные, стандартные.
Просто, был у меня в жизни период, когда я прошел так называемую "работу на работе". Это научило меня нескольким вещам.
Во-первых не доверять сторонним реалзациям, а если доверять, то покрывать их тестами, на 100% как BlackBox. Впрочем, лушче бы иметь открытую реализацию. Функции sin(x), cos(x) разложены в ряд тейлора, с корректно подобранными коэффициентами. Точности вычислений я знаю. Потенициальная ошибка может быть только работе процессора. Отчасти по этому выбран Тейлор. К тому же он FixedPoint.
Удалось достигнуть точности порядка 1LSB, при 32-битных вычислениях. Результат то все равно, фактически, 31-бит, т.к. еще есть знак.
Квадратный корень вычислял используя четыре аппроксимации. Первая,- линейная (всегда ниже оригинального значения). Вторая, апроксимация пораболы. Но пораболу следовало сместить, по этому аппроксимировал квадратично ошибку смещения пораболы по оси X. И наконец, апроксимация остаточной ошибки, которая дала точность результа +/-0.15. Последнюю, просто на глаз сделал, по графику. Поленился корни уравнения считать. А затем использовал итеративное приближение. Но это уже 7-12 бит дробной части из ее 16бит. Операцию деления не использую вообще. Ни при вычислении квадратного корня, ни при вычислении корня квадратного.
А что за метод неподвижной точки?
Как известно, каноническая форма любого уравнения это f(x) = 0. Здесь неявно предполагается, что число 0 входит в область опредения функции f(x), что в общем виде не так, более того эта функция вообще может быть определена не на числовом множестве. Тогда можно записать уравнение в другом виде f(x) = x. Это более универсальная форма, т.к. даже если нуля в области определения нет, то хоть какой то x уж точно существует)).
При этом можно заметить, что если f(x) = x, то f(f(x))=x и вообще f(f(f(...(x)...)))= x. Поэтому точка x называется неподвижной точной (fixed point), т.к. применение к ней функции f не изменяет саму x. А значит рекурсивное применение функции f самой к себе автоматически находит решение уравнения f(x)=x (либо не находит - т.к. в численных методах задача либо сходится, либо расходится))).
Как с помощью этого вычислить квадратный корень? По определению квадратный корень числа a это такое x, что x^2=a и x>=0. Запишем x^2=a в виде f(x) = x. f(x) = a/x = x. Стартуем от любой точки, например 1 и рекурсивно применяем f(x) к самой себе и... В результате выходит бесконечный цикл, в котором два значения x1=1 и x2=2 повторяются, прыгая вокруг правильного ответа. Решение расходится, печаль, бывает.
Но очевидно, что правильный ответ находится между x1 и x2, поэтому перепишем f(x) в таком виде f(x)=(x + a/x)/2. Видно, что f(sqrt(a))=sqrt(a), значит такая форма записи тоже удовлетворяет условию задачи. Попробуйте теперь рекурсивно вызывать f(x)=(x + a/x)/2, начиная с любого положительного x - сами увидите, что задача сходится очень быстро.
Я в эксельке накидал для наглядности. Уже на 4й итерации получился ответ с максимально возможной точностью:
К тому же в отличии от формулы Тейлора этот метод не накапливает ошибку вичислений. Т.к. каждый раз на вход f поступает неточное значение неподвижной точки x, и каждый раз f вычисляет её чуть точнее, тем самым автоматически подавляя ошибку.
Не уловил ход мысли. Подумаю еще.
Судя по последней формуле, Вы пришли к классическому вычислению корня квадратного, используя формулу применяемую в newlibc (__ieee754_sqrt(x)), т.е. формулу Герона. Ее еще в Египте применяли для этих целей. Дробление прямоугольника с длинной равной площади и высотой равной единице (начальное условие, в этом случае не 1.5, а 2 должно быть для корня из двух.
Однако.... При вычислении корня, таким образом, требуется операция деления. Там в, знаменателе стоит Y, а не степень двойки, а значит двоичным сдвигом уже не отделаешься. В моем случае, это деление еще и 64-битное. Деление,- достаточно дорогая операция, которая доступна не везде. В процессорах, не везде, а в FPGA и подавно.
Кроме того, Вы изначально очень удачно выбрали начальное условие, это непременно повлияло на количество итераций.
CORDIC- не не слышали?)
Да, Вы не первый кто упоминает Cordic. Обычно, это на слуху у ПЛИС-оводов, и тех кто занимается ЦОС. Сколько времени это займет на STM32?
У меня нет итеративной части при вычислении тригонометрии, она вычисляется всегда за одно и то же фиксированное время. Я понимаю, что быстрая математика на МК обычно использует таблиные методы с линейной апроксимацией. Это если говорить о ЦОС. Если говорить о моей задаче, то мне нужна определенная абсолютная точность. Это для станка CNC (дома стоит такой, фрезерный).
Если прочитать комментарии от первого до четвертого, то CORDIC прозвучит два раза. Ой, уже три.
А почему не использовать табличные значения, сгеренированные на этапе компиляции?
Так же хранить четверть диапазона и оставшиеся вычислять.
Да, я видел подобные методы. Если посмотреть готовые бибилиотеки cmsis, то можно найти быструю математику c числами в формате q31. Я смотрел их перед тем как писать код. Там, как раз таки таблица с линейной интерполяцией.
Собственно, приблизительно так и происходит в моем случае. Только я не интерполирую табличные значения. Эти значения вычисляются каждый раз, используя ряд Тейлора. Т.е. имеется интерполяция sin(x), на интервале 0..PI/4. Имеется интеполяция cos(x), на интервале 0..PI/4. Аргумент приводится к этому интервалу и нужной функции, и вычисляется точное значение для каждого аргумента. Я просчитал. Точность 32-битного числа, составляет 2.32*10^-10. Я просчитал sin(x) и cos(x) приблизительно с этой точностью на указанных интервалах. Это обычный sin(x), cos(x), но просчитанный в фиксированной точке, и смасшатбированный в диапазон 0..2^32 или [-2^32..2^32-1].
Я тоже большой фанат фиксированных вычислений, они предсказуемы и детерминированы на всех платформах. Мне показалось что в моем случай точности pi/512 за глаза поэтому я использовал табличный подход. Для корня я использовал побитовый поиск, быстр и прост. Но потом понял что нативный способ все равно быстрее и оставил так:
#[inline(always)]
fn sqrt(self) -> FixFloat {
let a = (self as f64).sqrt() as i64;
let b = a + 1;
if b * b > self { a } else { b }
}Еще вопрос как вы храните дробную часть? У вас помножено на какую-то константу?
Начальное приближение корня квадратного, в моем случае, вычисляется аналитически. Это требует 9 умножений.
Первое,- это линейная апроксимация на участке [2^n..2^n+1]. Максимальная абсолютная ошибка результата линейной аппроксимации, порядка 500 (если правильно помню). Линейная аппроксимация требует 1 умножение (обычный kx+b). Коэффициенты при этом являются константами в формате n*sqrt(2)^m или n*2^m, в зависимости от того четное или нечетное число бит аргумента поступило на вход.
Остаток от линейной аппроксимации имеет квадратичную форму. Это порабола Однако, парабола смещенная. Смещение параболы квадратичной ошибки тоже имеет квадрвтичный характер по отношению к аргументу. Корректировка аргумента, требует двух умножений. Первый это вычисление пораболы, второй,- коэффициент масштабирования корректировки. Это разница между центральным значением интервала [2^n..2^n+1], и максимумом аппроксимации. Затем еще два умножения. Это вычисление смещенной пораболы аппроксимации квадратичной ошибки, и некий коэффициент,- максимум квадратичной ошибки. Подобный подход позволил добиться абсолютной точности результата, порядка 6, для 32-битного аргумента.
Затем, аппроксимировал остаток ошибки, уже на глаз. Впринципе, можно посчитать производные, найти максимумы и посчитать точнее. Это еще увеличило точность результата до 0.15. Таким образом начальное значение квадратного корня уложилось в абсолютную погрешность 0.15.
А затем, уже использовал бинарный поиск. Вообще, результат функции sqrt от 32 битного аргумента, в случае получения целочисленного результата, имеет только 16-бит. Я вычисляю 32-бита результата. 16-бит,- целая часть. И 16-бит,- дробная. Дробная часть при начальной аппроксимации вычисляется с точностью 2-9 бит. Т.е. остается бинарным поиском найти 13-7бит.
Константы хранятся непосредственно в коде, зачастую, как пары чисел, сдвинутые на степень 2-ки.
У меня на работе лежит - "Digital Signal Processing Applicaions using the ADSP-2100 Family" от Analog Devices. Там масса информации на которой я учился писать подобные функции. Сама фраза про ряд Тейлора как бы немного уже напрягает. В книге почти как бы то же самое, но используются коэффициенты, которые ищутся на минимизации средне-квадратической ошибки на всем интервале аппроксимации.
У меня, небыло такой книжки. Интересно почитать. Вообще, занимался разработкой кода для микроконтроллеров, а не DSP и не FPGA. FPGA,- это больше уровень хобби. Впрочем понимаю и FSM и синхронную логику. А FixedPoint как-то не требовался раньше. Что касалось "работы на работе", оно обходило меня стороной, и я просто выполнял работу инженера, вот и не лез так глубоко, считая что везде все хорошо. А оказалось,- не всегда и не везде.
Вопрос по теме.
Какая степень аппроксимирующего полинома получилась у Вас при аппроксимации синусойды в 32-битах с точностью, не превышающей 1LSB на всем периоде? У меня максимум 12. Но коэффициенты идут через 1, по этому умножений на степень X меньше (домножение идет на X^2).
Т.е. сколько операций умножения требуется для реализации функции в этом случае?
По правде говоря, я рассматривал аппроксимацию в других точках, т.е. (X-Xn). Но, как оказалось, минимальное количество умножений при разложении в нуле.
Идей там немного и они достаточно просты.
Надо привыкнуть представлять все во fractional format (начальные главы книги..).
Диапазон - [-pi, +pi] - это в int16_t [-(2^15)... (2^15 -1)] - (популярный вариант)
И в большинстве мест придется считать в Q15 (положение разделительной точки..)
Далее мы составляем в матлабе
s = (0: 1/32: 1.0)';
A = [ ones(32,1) s s^2 s^3... ];
Тут "32" и степень аппроксимации - пока от фонаря.
b = sin( (pi/2) *s);
k = A \ b - получаем коэффициенты аппроксимации данного куска синуса от разных этих степеней
Насколько классно получилось это все -
ba = A* k
Тут можно взять разницу исходного синуса и аппроксимированного, потусовать степени полинома, прикинуть желаемую точность
В моей библиотеке есть вариации разной точности - 12 бит (табличный), 16 бит - полином 5 степени, 18 бит - в кодеке AMBE потребовалось...
Саму это книгу вроде бы не сложно было найти в pfd. Я ее распечатал для себя 20 лет назад.
Если достаточно 2% (то есть мы выпиливаем какой-нибудь кружок с точностью 20%, в котором 2% просто утонут) — алгоритм Бхаскары I рулит.
Геометния, это геометрия. В общем и целом, используя лишь апроксимацию я получил точность порядка 0.003%. И на это ушло только 9 умножений, сложения, вычитания и сдвиги. Это о корне квадратном.
Синус и косинус, считаются как полином 12 степени с отбрасыванием младших значащих бит фиксированной точки, после умножения. Я сейчас не дома, и не могу посмотреть. Но, помоему, количество произведений там 6 для вычисления степнеи x + столько же для домножения на коэффициенты. Куда еще проще? Никаких итераций и точность порядка 1младшего разряда.
Добрый день, я не из этого раздела инженерии и не понял, подскажите пожалуйста, в чем проблема если использовать табличные значения этих функций?
Добрый день!
Спасибо за вопрос.
Если говорить о тригонометрии, то она вычисляется для любого X. Поскольку аргумент функции 32-битный, то окружность разбита на 2^32 степени интервала. Это очень большое число. Держать 4 * 2^32 байт памяти на микроконтроллерах непредставляется возможным. Там память ограничена единицами-сотнями тысяч байт.
Есть альтернативные варианты,- это держать табличку на некоторое количество значений, а между ними использовать апроксимацию. Обычно для ускорения вычислений, величина таблицы выбирается достаточно большой. А апроксимацию выбирают линейную. Это быстро, но, как правило приводит к некоторой погрешности. В зависимости от размера таблицы. В этом случае табличка тоже занимает определенный обьес памяти. Так ревлизована математика в библиотеках cmsis от Arm. Ничего плохого в них нет.
Представьте, что Вы не можете опереться ни на что недоказанное при вычислении sin/cos. Как его вычислить? Все просто. Производные и ряд. Кроме того, Вы точно знаете точность (остаток полинома), который используете. Наверное, по этому.
Если посмотреть на корень квадратный, то он частично табличный. Первой версией, была версия, основанная только на битовых сдвигах, и линейной аппроксимации. Там вообще небыло каких-либо констант, а аппроксимация производилась между четными степенями двойки. Естественно, подобная точность меня бы никогда не устроила. В первой итерации, я пришел к апроксимации по степеням двойки, отсюда появилась костанта степени корень из двух.
Затем, посмотрел на результат и аппроксимировал ошибку, используя полином второй степени. Впринципе, на тот момент, рассматривал аппроксимацию разлмчными рядами сразу. Однако, динамический диапазон вычичлений получается широким. Для операций с плавающей точкой, это приведет лишь к увеличению ошибки. А при фиксированной точке, к увеличению разрядности (если требуется сохранить точность).
Я пробовал аппроксимировать ошибку рядами и получал более высокую точность сразу. Однако, это привело к трате ресурсов. По сравнению с итератианым подходом не давало никакого выыигрыша. Аппроксимация сдвига пораболы ошибки увеличила скорость вычислений.
Суть понял, интересный кейс оказался, люблю с таким поковыряться. Попробую подумать как еще можно
Сделал вот заготовку чтобы не пытаться плодить члены Тейлора, а обойтись для корня интерполяцией полиномом 4-й степени и короткой корректировочной таблицей. Я вот это имел ввиду говоря про таблицы. Для 256 точек пока. Если интересно продолжить, напишите в личку.
Примерно понимаю как такой же фокус провернуть для чисел 16-бит и более, попробую придумать.
Для большого количества точек еще надо будет придумать как автоматически границы диапазонов определять.
Еще можно небольшую библиотеку сделать начиная с умножения, и деления, а потом ее на комплексные числа развить. Предполагаю, что вы что то подобное делаете.
Добрый день!
В общем и целом, интересно. Единственное, мне придется освежить в памяти знания по аналитической геометрии. Все-так 20-лет уже этим разделом мат-анализа не пользовался.
Вижу, сначала, Вы используете матрицу с числами i^j. Коэффициенты вычисляете, используя транспонированную и обратную матрицы. Произведение на 10^9, полагаю это округление до 2^32 приблизительно.
Почитаю метод, вспомню, и отвечу Вам. Спасибо.
Там все несколько проще, чем учебник ангема) sin=y/r, cos=x/r, r=(x^2+y^2)^(1/2). А так обращайтесь за материалом, я на хабре как раз в основном про трехмерную геометрию пишу в кватернионах. Формула с матрицей это из линейной алгебры - метод наименьших квадратов заданный через псевдообратную матрицу.
То, что в прямоугольном треугольнике длинна гипотенузы равна сумме квадратов катетов, знают даже дети. Собственно r*sin(x),- есть вериикальна сторона треугольника, r*cos(x) горизонтальная, а r,- это длиина гипотенузы известно было еще в школе.
Я, имел ввиду матричные вычисления. У нас в университете этот предмет назывался "Аналитическая геометрия". Определители матриц, обратные матрицы, численные методы решения систем уравнений, и так далее. Последний раз я работал с этим, когда учился в университете, году в 2003. Сейчас 2025, т.е. 22 года. Если работать, то следует вспомнить методы.
У меня, задачи немного другие на данный момент.
Расскажите, что за задача, что требуются именно функции от углов?
Если говорить о математике, то всего-лишь для станка CNC, без использования FPU, т.е. сопроцессора с плавающей точкой. Просто хочу быть уверен что все работает достаточно быстро и именно так, как было задумано. Возможно, потом переложу на ПЛИС, потому и не использую деления.
В общем и целом, был в моей жизни период, когда меня очень серьезно атаковали. В тот момент, очень хотелось эти атаки прекратить. Вектором атаки могло быть что угодно. И ОС, и BIOS и даже Intel ME или прошивка оборудования.
Тогда, я очень злился и думал как бы сделать так, чтобы в систему на которой работаю никто не вмешивался. Я отключал оборудование от сети, фильтровал питающие кабеля, глушилку использовал внутри системного блока. С программной точки зрения, сначала сменил ОС на Linux, потом фильтровал трафик, затем вообще отказался от сети.
А атаки продолжались и продолжались. Причем атаки слишком умные какие-то. На уровне замены компилируемой мной прошивки после ее сборки, смены библиотек и подобных вещей.
Я ничего не знал о симуляторе и плохо разбирался в людях тогда. А о "работе на работе",- именно в кавычках, не знал вообще ничего. Так вот, в тот период, мне банально мешали сделать хоть что-то. И я понял, что если хочу получить что-то надежное, то придется все делать с основ.
Была даже мысль x86-подобную систему реализовать на FPGA. Использовать минимальный Linux (LFS), и работая в консоли собирать то, что я хотел под целевую систему. Другими словами, шансов у одного против команды / людей у меня небыло.
Вот, с тех пор, вынашиваю идею по созданию защищенного процессора с нуля. Я имею ввиду не защиту памяти на уровне MMU и исключений, а вообще, аппаратно близкую к обьектно ориентированному подходу систему.
Однажды на youtube, обозначил некоторые пункты, которые хотел бы реавлизовать. Основной смысл в том, что каждый вызов функции это вызов объекта за пределы которого функция выйти не может. Отсутствие возможности прочитать свой код, и уж тем более модифицировать его. Вообще без таких команд. Все константы и таблицы в отдельном сегменте только для чтения. Все регистры составляющие контекст функции/задачи сохраняются в свой стек (одинаковой глубины) командами call/ret, и прямой доступ к указателю этого стека возможен лишь из специального режима планировщика задач. При этом нет разницы задача это или функция, и скорость вызова функци равна скорости переключения задач (не совсем так но почти).
Меня раскритиковали, говоря что сейчас самое главное,- производительность, а то о чем я говорю дорого с точки зрения производительности.
Вернемся к теме Авиации, где я хотел бы работать. Так вот, это для домашних проектов можно тяп-ляп, или в коммерции можно тяп-ляп и в продакшн. Для авионики существуют определенные стандарты. Начиная с Misra-C и подобных. Я понимаю, что авиационное ПО, проходит множество проверок и покрывается тестами перед тем как пойдет в продакшн.
Однако, то с чем столкнулся я это не вина компилятора, и не вина того кто зашивает программатором прошивку в изделие. Если бы подобная атака произошла на технологическое приспособление в авиации, то в ряде случаев управление самолетом могли бы в худшем случае перехватить, а влучшем он бы упал. Понятно, что я о критически важных узлах говорю. Отсюда интерес в том чем пытаюсь заниматься. Пытаюсь,- поскольку это вопрос финансирования.
Однако, вопрос финансирования решать не пытаюсь, поскольку привлечь прямого инвестора нет возможности. А все непрямые сейчас смываются и перекладывают вину на тебя. Иногда, даже используя твой счет без тебя чтобы вывести собранные средства. В ход идут и поддельные proxy-карты и многое другое. Так что ечли работать то только с финансистом напрямую.
Ситуация приблизительно такая.
Вот аппроксимация квадратного корня полиномом 5й степени с точностью ~10^(-5). Получил через узлы Чебышева, в диапазоне [1,2], c пяти-шестизначными коэффициентами. То есть вам надо 18 бит, чтобы закодировать числа, и раз вас обработчик процессора не устраивает, то нужно описать преобразования прямо в вашем коде, включая преобразования экспоненциальной части чисел.
Не знаю чем вам не нравилась операция деления, ее тоже сделал для этого же диапазона на всякий случай, полином 6й степени получился.
Посмотрите на qfplib-m3. Может подойдет по ваши цели и задачи
Я посмотрел на библиотеку. Возможно, для многих сечас это не так актуально, но она использует как минимум sdiv. Т.е. деление. У меня в коде, присутсивует умножение 32x32, но деления нет вообще.
Если быть до конца честным, то и умножение хотелось бы свести к меньшей разрядности. Например 18x18. Подобные умножители есть у GoWin в их FPGA.
Вы спросите, а зачем мне это вообще? Я стремлюсь к основам. Т.е. к фундаменту на который можно опереться и не боятся того что оно потом развалится.
Я не могу изготовить чип, поскольку не обладаю доступом к производству. Да и между кодом на Verilog и топологией кристала огромная пропасть.
Но я могу сделать модель чипа с тестами на Verilog, оттестировать ее тоже могу. Кроме того, сейчас доступны недорогие FPGA, на которых эту модель можно запустить.
Синус, косинус, квадратный корень FixedPoint