Считаем синус быстро и точно: ч.3 - Практика

Картинки не будет, не нашёл подходящей…
Как это использовать на практике
Для использования аппроксимации в своём проекте для начала нужно определиться, с какой точностью нужно вычислять синус. Потом нужно оценить ресурсы (объём памяти и быстродействие вычислителя), которые у вас есть.К примеру, вы хотите, используя 32-битный вычислитель, сделать генератор с 24-битным ЦАП. Старший отведён на знак, поэтому точность аппроксимации должна быть не меньше 23 бит. На диаграмме видно, что заданная точность достигается при:
таблице на 8192 строки при использовании аппроксимации 1-й степени,
- 512 строк при 2-й степени,
- 64 строки при 3-й степени,
- 32 при 4-й степени,
- 16 при 5-й степени, и
- 8 при 6-й степени.
Теперь следует определиться с ресурсами и приоритетами их использования. Если вы хотите получить максимальное быстродействие, нужно использовать аппроксимацию 1-й степени, таблица на 8192 строки займёт 64к. Если вы используете какой-нибудь STM32 на Cortex-M3, то размещение таблицы в ОЗУ уменьшит время вычисления в отличии от размещения её в памяти программ. Возможно, это работает и для других МК, я проверял только на STM.
Если у вас очень мало памяти и есть запас по быстродействию (FPGA например), можно использовать аппроксимацию с высокими степенями. Таблицы 5-й и 6-й степеней займут, соответственно, 16 * (5+1) * 4 = 384 и 8 * (6+1) * 4 = 224 байт. Уменьшить объём ещё в 4 раза можно, используя симметричность функции синуса, правда точность при этом может упасть, т.к. при аргументе функции 0x40000000, что равно 90°, будет использован 0x3fffffff. Для использоания этого фокуса, при сборке проекта в опциях компилятора нужно задать дефайн Q1_ONLY — это уменьшит размер таблицы и добавит лишние вычисления в функции get_sine_int32() и get_sine_float().
Таблицы с коэффициентами
Их можно получить программой make_table. По ссылке можно найти файл проекта на code::blocks, исходные тексты и скрипт для самостоятельной сборки с помощью gcc, а так же готовый .exe.
В качестве аргументов программа принимает:
- размер таблицы, равный 2^n;
- степень полинома;
- множитель для коэффициентов для приведения их к целочисленным значениям;
- параметр AC_SHIFT, показывающий, на сколько бит множитель при коэффициенте, например, 3-й степени будет «весомее» множителя предыдущей 2-й.
Результатом работы программы является исходный файл на языке си. Программа выводит его непосредственно на экран, его можно перенаправить в файл. Например:
.\make_table.exe 64 3 0x40000000LL 3 > ..\tables\sine_approx_64_3_3.c
создаст файл sine_approx_64_3_3.c. Этот файл содержит таблицу коэффициентов для аппроксимации полномами 3-й степени по 64 отрезкам, чего достаточно для получения 25-битной точности. Множитель при коэффициенте A0 = 0x40000000, а множители при коэффициентах более высоких степеней будут отличаться от соседних в 2³ = 8 раз.
В начале файла есть набор директив препроцессора, с помошью которых можно настраивать некоторые параметры таблицы:
- TABLE_TYPE — тип данных таблицы, по умолчанию int32_t. Его можно «перебить», указав компилятору другой.
- Q1_ONLY — указание использовать только первый квадрант, позволяет уменьшить объём таблицы в 4 раза.
- DC — «добавка» к ко всем коэффициентам перед превращением их в целочисленный вид, позволяет делать арифметическое округление. По умолчанию равно 0.5 (для int32_t), для чисел с плавающей точной следует задать его равным 0.
- AC0 — множитель для коэффициента 0-й степени, по умолчанию 0x40000000LL. Может быть задан другой, для таблицы типа float можно задать 1.
- AC_SHIFT — на сколько бит будет сдвинут влево множитель AC0 для каждой последующей степени, механизм этого понятен при формировании множителей AC1, AC2 и так далее. Для таблицы длиной 64 строки опттимальное значение 3, для таблиц другой длины оптимальные значения лежат в районе от 1 до 5.
Каждая строка полученной таблицы выглядит примерно так:
-0.00015749713825 * AC3 + DC, -0.000000170942 * AC2 + DC, +0.098174807817 * AC1 + DC, -0.000000001187 * AC0 + DC, // 0
где
- -0.00015749713825096520473455514 * AC3 — коэффициент для x³,
- -0.00000017094269773828251637917 * AC2 — для x²,
- и так далее.
В конце файла, после таблицы, сформирован отчёт, какой «запас по битам» у каждого коэффициента. В процессе рассчёта множителей находятся таковые с максимальными абсолютными значениями. Для множителя 0 порядка это 0.999999969786949, что почти рано 1.0, и запаса по битам нет. А максимальный коэффициент 3-го порядка равен 0.00015749713825, и его можно безболезненно умножить на 2^12, и результат будет всё ещё меньше 1.0. С помощью запаса по битам при желании можно скорректировать множители AC1… ACn.
В папке проекта есть скрипт на lua, который сформирует оставшиеся файлы, необходимые для включения в проект, и сформирует cmd-файл, запустив который получатся 90 файлов с таблицами 4..65536 строк и степенями от 1 до 6.
Примеры использования
Есть 2 готовых примера — для ПК и для популярного STM32F103C8.
Они реализуют цепочку вычислений умножение — сложение — умножение — сложение и т.д., формула
y = A2·x² + A1·x + A0
преобразована к виду
y = ((0 + A2)·x + A1)·x + A0
что бы избавиться от лишних умножений для возведения в степень.
ПК
Там есть 2 пары файлов — для целочислееных вычислений get_sine_int32 (.h/.c) и для вычисления с плавающей точкой get_sine_float (.h/.c), кторые нужно включить в свой проект. Так же в проект нужно включить нужную таблицу с коэффициентами.
Использование таблиц и функций с плавающей точкой оправдано, т.к. это позволяет ускорить вычисления.
Cortex-M3
Проект на embitz для STM32F103C8.
Вычисления синуса для полиномов от 1-й до 6-й степеней оформлены в виде ассемблерных вставок в код на с. Перед каждой вставкой в комментариях есть примерно измеренное кол-во тактов на 1 цикл вычислений. Я не считаю себя знатоком ARMовского ассемблера, возможно кому-то удасться вычисления ускорить.
Размешение таблиц в ОЗУ, как уже было сказано, позволяет немного поднять быстродействие.
Здесь надо рассказать о том, почему использование коэффициента AC_SHIFT ускоряет вычисление. У Cortex-M3 есть инструкции умножения двух 32-битных чисел с получением 64-битного произведения: UMULL и SMULL. Первая умножает два беззнаковых числа, и нам не подходит. Вторая, которую мы используем, умножает знаковые числа: одно из них — это коэффициент из таблицы, знаковое число, второе — это смещение угла внутри интервала, число беззнаковое. Если второе число превысит 0x7fffffff, то оно будет интерпретировано как отрицательное и результат будет не тот, который мы ожидаем. Инструкции, которая перемножала бы знаковое и беззнаковое число, в наборе команд нет, и компилятор вынужден в таких случаях городить огород.
Избежать этого можно, если оставить второе число знаковым, но не позволить ему превысить 0x7fffffff, например сделать его в 2 раза меньше. При этом надо будет после каждого умножения корректировать результат. Сделать это можно двумя способами:
- каждый раз удваивать результат, на что уйдут машинные такты;
- сделать в два раза больше то число, которое на него умножаем.
Более экономичным является второй способ. Если x уменьшен в 2 раза, то из формулы
y = (((0 + A3)·x + A2)·x + A1)·x + A0
понятно, что в 2 раза должен быть увеличен кэффициент A1, в 4 раза A2, и в 8 раз A3.
Если же x уменьшен в 4 раза, то множители становятся соответственно 4, 16 и 64.
Именно за это при компиляции таблиц отвечает параметр AC_SHIFT. Кроме того, его использование убивает второго зайца — коэффициенты, содержащие много нулей после запятой — а они относятся к высоким степеням полиномов, не так сильно теряют при округлении.
За сим позвольте откланяться, и спасибо за внимание!
Часть 1 — Теория
Часть 2 — Точность вычислений
- +6
- 08 октября 2021, 14:49
- vix
Комментарии (1)
RSS свернуть / развернуть