Библиотека математических операций с фиксированной точкой

[menzoda]

Вот мой профиль на гитхабе, там две библиотечки, 16-битная и 32-битная. Для каждой добавил несколько Issues чтобы было понятно что я собираюсь делать в ближайшее время. Если вкратце, то я наконец разобрался с делением. Для 32-битной еще нужно доделать квадратный корень. Для 16-битной - синус, косинус и тоже квадратный корень. Потом думаю перенести 16-битную на Cortex-M0 для которого она изначально и предполагалась, просто начал с M3 потому что так легче переносить. Затем, можно обе перенести на Cortex-M4. В принципе, версия для M3 соберется и для M4, но будет не совсем оптимальна.

[menzoda]

Значит дела обстоят следующим образом: 32-битную версию для Cortex-M3 и Cortex-M4 я "релизнул", 16-битную версию для для Cortex-M3 и Cortex-M4 решил отложить и перейти сразу к Cortex-M0 на котором она нужнее. Думал сейчас быстренько все сделаю и тогда уж скажу. Однако, когда я покопался в наборе инструкций Cortex-M0, то оказалось что он настолько урезанный (читай убогий), что просто так взять и быстренько сделать не получится. Он чуть ли не как PIC, просто с регистрами 32-битными. Такой засады я конечно не ожидал от него.

[Andrew Martin]

Такой засады я конечно не ожидал от него.

Ато!
Тут то как раз и необходим асм.

Работаю над БПФ для Cortex-М0 на ассемблере. За основу взято вот это:

fix_fft.c

(11.27 KiB) Скачиваний: 467

только будет пока только прямое БПФ, и размер его задаётся константой.

Код уже работает на тестовом массиве, результат вполне адекватный.
Одна печаль - много мусора. Там, где должны быть 0, периодически попадается то +1, то -1. Так же и амплитуды гармоник отличаются от истинных значений на +/- 1 единицу.
Вероятно это вылазит ошибка округления при арифметическом сдвиге вправо. Как её побороть пока не знаю.

[Andrew Martin]

Дошли руки посмотреть ваши наработки по М0. Функция умножения

Код:

            muls    r0, r1, r0         ;1
            movs   r1, #0             ;1
            asrs     r0, r0, r2         ;1
            adcs    r0, r0, r1         ;1
            bx       lr                    ;3


неоптимальна уже тем, что на её вызов+возврат тратится 4+3=7 тактов, а собственно на умножение с нормализацией - 4 такта.
Кроме того, при вызове функции несколько раз можно выделить "нулевой" регистр один раз, и всё время использовать его для округления. Тогда N умножения с нормализацией займут 3*N+1 такта, а не 11*N тактов, как сейчас.
И это мы ещё не учли расходы на запись аргументов функции в нужные регистры... Здесь особенно бросается в глаза передача числа битов для сдвига, которую вообще можно забить константой для выбранного формата чисел.

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

[menzoda]

Для М0 наработок нет. Я сделал функцию умножения, потом перешел к делению, увидел что у него нет инструкции CLZ, что только лишь её реализация уже будет занимать 1/3 от всей функции деления, и забросил это дело до лучших времен.

на её вызов+возврат тратится 4+3=7 тактов, а собственно на умножение с нормализацией - 4 такта.

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

Кроме того, при вызове функции несколько раз можно выделить "нулевой" регистр один раз, и всё время использовать его для округления.

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

И это мы ещё не учли расходы на запись аргументов функции в нужные регистры...

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

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

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

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

Я думаю как раз наоборот. У M3 есть интересные инструкции и возможности, которые компилятор никогда не догадается использовать, особенно в комплексе. У М0 же все сделано очень топорно, ручная оптимизация практически не даст эффекта.

[Andrew Martin]

Выбор конечно за вами, просто в М0 и так полноценных регистров всего 8, не разгонишься особо... А локальные переменные для больших функций тоже нужно будет базировать на чём-то, итого остаётся 7 регистров.

[menzoda]

Посмотрел-подумал и решил объединить репозитории для отдельных архитектур в один общий. Кроме того добавил поддержку GCC. Еще прикинул как реализовать экспоненту и тангенс, вышло не очень хорошо. Для тангенса наверное нужно будет две lookup таблицы, а с экспонентой... тоже не все гладко.

[Andrew Martin]

Что, некрасивое разложение в ряд Тейлора получается?
Да и ну его, этот тангенс вместе с экспонентой. Гораздо полезнее арктангенс отношения и десятичный логарифм.

[menzoda]

Ага, некрасивое. Точнее, в нуле оно у всех нормальное, а вот в окрестности произвольной точки у синуса/косинуса оно хорошее, а у тангенса не очень. Я имею ввиду что производные синуса и косинуса это те же самые синусы и косинусы, и их можно получить из таблицы. А вот у тангенса первая и вторая производные разные, их из одной таблицы никак не получишь, нужно две. С экспонентой в этом плане нормально, любая производная будет той же экспонентой, проблема в другом. В реализациях других функций входной аргумент сначала нормализуется, потом вычисляется значение функции, и в конце это значение обратно денормализуется. Это чтобы вычислять с максимальной точностью. В делении, синусе и косинусе, денормализация - это просто битовый сдвиг. В квадратном корне - сдвиг и умножение на константу. А вот в экспоненте все сложнее, тоже наверное надо будет вторую таблицу вводить. Логарифм и арктангенс по прикидкам действительно будут полегче реализовать. К тому же, если они действительно полезнее...

[Andrew Martin]

Да понятно что разложение в окрестности точки, интерполяция то рядом Тейлора между точками таблицы. Я по мотивам ваших функций для себя сделал синус и косинус для М0 и только на 16 бит. Получилось отнюдь неплохо.

Уточню: арктангенс отношения и десятичный логарифм отношения. Первый - для вычисления фазы комплексного аргумента, второй - для перевода в децибелы. Многие наверное согласятся, что это будет полезнее даже синуса/косинуса, а тем паче тангенса.

А квадратный корень можно сделать и без таблиц, но тогда нужен цикл, который намотает тем больше тактов, чем больше бит. За то ошибка вычислений будет минимальной. На М0, где нет многого, в том числе и длинных умножений, это может быть неплохим вариантом.

[menzoda]

Уточню: арктангенс отношения и десятичный логарифм отношения. Первый - для вычисления фазы комплексного аргумента, второй - для перевода в децибелы. Многие наверное согласятся, что это будет полезнее даже синуса/косинуса, а тем паче тангенса.

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

[bad2cat]

Вам надо поглядеть:
1. библиотеку математики для 8080.
Я не шучу. использование 8-ми битных функций для вычислений 16-ти битной и 32-х битной арифметики. Это очень поучительно. Есть такая переводная книга "8085 subroutines".
2. подпрограммы БПФ на С, с алгоритмами Уэлча (могу ошибиться скорее всего). Главное, какой подход к оптимизации предпринят - во-первых вычисление индексов элементов входного массива как двоичный сдвиг. во-вторых использование только целочисленной арифметики, без фиксированной или плавающей точки (программы библиотек для БПФ от МИТ - давненько было - могу спутать - fft.c с разными алгоритмами перемешивания элементов входного массива, прореживание по времени или частоте), в-третьих, умножение не на синусы, а на минус 1 или плюс 1, но входное окно подбирается таким образом, чтобы сретушировать дефекты скоростного вычисления, фактически это одна серия умножений, вместо 10-11 (для 1000 точек).

[Andrew Martin]

А чем целочисленная арифметика отличается от арифметики с фиксированной точкой?

[menzoda]

Ага, я тоже споткнулся об эту фразу, но спрашивать было как-то лень

[Andrew Martin]

Вот неплохой на первый взгляд алгоритм, там, где не подходит ни Ньютон, ни Тейлор.

[menzoda]

Вроде такого способа вычисления до этого не встречал. Посмотрю на досуге.

[Alex_Megavolt_79]

Добрый вечер. А как библиотеку подключить? Я по ходу делаю чего-то не так, 476 ошибок вылезло.Вот часть их

Спойлер

fix32-master\src\cm3\sin.s(2): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(3): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(4): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(5): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(7): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(8): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(11): error: A1163E: Unknown opcode Returns , expecting opcode or Macro
fix32-master\src\cm3\sin.s(12): error: A1163E: Unknown opcode a , expecting opcode or Macro
fix32-master\src\cm3\sin.s(13): error: A1163E: Unknown opcode in , expecting opcode or Macro
fix32-master\src\cm3\sin.s(14): error: A1163E: Unknown opcode in , expecting opcode or Macro
fix32-master\src\cm3\sin.s(15): error: A1167E: Invalid line start
fix32-master\src\cm3\sin.s(16): error: A1163E: Unknown opcode Execution , expecting opcode or Macro
fix32-master\src\cm3\sin.s(17): error: A1163E: Unknown opcode Absolute , expecting opcode or Macro
fix32-master\src\cm3\sin.s(19): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(20): error: A1137E: Unexpected characters at end of line
fix32-master\src\cm3\sin.s(22): error: A1167E: Invalid line start
fix32-master\src\cm3\sin.s(23): error: A1105E: Area directive missing
fix32-master\src\cm3\sin.s(23): warning: A1088W: Faking declaration of area AREA |$$$$$$$|

[bad2cat]

Andrew Martin писал(а):

А чем целочисленная арифметика отличается от арифметики с фиксированной точкой?

1.Степень не обрабатываете
2.после каждого раунда (серии бабочек) результат мантиссы сдвигаете вправо (делите на два) - операция элементарная и очень быстрая, но не даёт переполниться сумматору.
3.умножение только сдвигами, т.е. умножение на 2,4, 8 (где-то есть такой вариант БПФ, специально для целочисленной арифметики), поищите.

Я перестал заниматься БПФ, после обрушения диска в 2001.
А там были мои наработки - во-первых АОН из звуковой карты на PC и на 8085, во-вторых попытка проработки анализатора спектра (АС) на компе.

Заготовка для АС на Duron 900 МГц делала 800 кадров по 1024 точек в секунду, а на Intel Celeron 600 МГц она просчитывала по 700 кадров. Но это без вывода на экран - в память.

Andrew Martin писал(а):

Да понятно что разложение в окрестности точки, интерполяция то рядом Тейлора между точками таблицы. Я по мотивам ваших функций для себя сделал синус и косинус для М0 и только на 16 бит. Получилось отнюдь неплохо.

Это лишнее.
Вы слишком увлекаетесь программированием и перестаёте замечать простые вещи. Которые к тому же работают быстрее чем ваши. Это "опьянение" от хорошей, быстрой, современной техники, прощающей ошибки программирования. Поучитесь у старичков, которые раньше на примитивных процах работали.
В вашем случае достаточно линейной интерполяции и не надо никаких рядов Тейлора.

menzoda писал(а):

Во-вторых, я ради примера посмотрел на их реализацию sin (PU).В-третьих, они используют таблицу на 512 значений. 512 значений, Карл! Каждое по четыре байта, всего 2 Кб памяти. .

Ага, вот сейчас не поленился и провёл исследование по точности аппроксимации.
Между таблицей, её ЛИНЕЙНОЙ аппроксимацией и настоящим синусом (т.е. рядами Тейлора) в экселе.

Предположим надо точность 10^-6, чего на практике хватает за глаза для любых случаев - будь то Фурье (там вообще иногда однобитное преобразование делают , будь то построение метровой окружности для распечатки с большим dpi.
Я просто беру таблицу 90 градусов через 0,1 градус, 900 слов. Что в 32-х битной системе - 900 х4 = 3,6 кБ.

Нужно сравнить точность в ключевых точках, там где производная функции (а значит отклонение) максимально. Это точка 45 градусов (может быть, или начало и конец, т.е. 0 или 90 градусов).

Так вот максимальное отклонение для такой таблицы всего 0,38 ppm в точке 90 градусов (или 0 градусов для косинуса). А минимальное значение 0,00033 ppm в точке 0 градусов. Сравнение производится в середине интервала, разумеется.
А для таблицы с шагом 1 градус размер таблицы будет 360 байт, а наихудшее отклонение от идеального синуса (в точке 90 градусов) - 38 ppm. Внизу в файликах - оба случая для трёх точек аппроксимации.

Стоит ли городить сложные и медлительные аппроксимации по рядам Тейлора? Вряд ли.

Вложения

snag-2017-07-010-a.JPG

(60.1 KiB) Скачиваний: 631

snag-2017-07-009-a.JPG

(106.85 KiB) Скачиваний: 607

Последний раз редактировалось bad2cat Чт июл 13, 2017 09:51:19, всего редактировалось 1 раз.

[jcxz]

bad2cat писал(а):

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

На 32-битном МК для большинства случаев достаточно полиномиального вычисления синуса - всего 7 MAC-операций чтобы получить ~13 бит точности.
И не надо никаких таблиц и аппроксимаций.
А если уж хочется ещё быстрее и нужна именно генерация синуса (вычисление значений с неким угловым шагом), то можно считать одно значение полиномом, а потом N значений -
рекуррентно тратя всего одну(!) MAC-операцию на одно значение синуса (достаточно вспомнить школьный курс тригонометрии).

[bad2cat]

jcxz писал(а):

для большинства случаев достаточно полиномиального вычисления синуса - всего 7 MAC-операций чтобы получить ~13 бит точности.

Вот надо всегда оговаривать какой случай вас интересует - быстродействие или компактность кода. А то вы смешиваете две задачи. А на практике всегда приходится выбирать компромиссный путь.

jcxz писал(а):

И не надо никаких таблиц и аппроксимаций. А если уж хочется ещё быстрее и нужна именно генерация синуса (.

Надо таблицы, надо. Речь как раз идёт о скорости вычислений. И не для 13 бит (т.е. 4,8 десятичных цифры), а о разработке библиотечной функции с точностью 9 десятичных знаков, т.е. 32 бит. В этом случае таблицы быстрее работают.