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

Наконец-то переходим к самой важной части — как это вот всё использовать. Здесь я расскажу, как выбрать подходящую комбинацию длины таблицы и степени полинома, как создать эти таблицы и приведу примеры, как посчитать быстро (ну или относительно быстро).
Картинки не будет, не нашёл подходящей…

Как это использовать на практике

Для использования аппроксимации в своём проекте для начала нужно определиться, с какой точностью нужно вычислять синус. Потом нужно оценить ресурсы (объём памяти и быстродействие вычислителя), которые у вас есть.

К примеру, вы хотите, используя 32-битный вычислитель, сделать генератор с 24-битным ЦАП. Старший отведён на знак, поэтому точность аппроксимации должна быть не меньше 23 бит. На диаграмме видно, что заданная точность достигается при:
таблице на 8192 строки при использовании аппроксимации 1-й степени,

1. 512 строк при 2-й степени,
2. 64 строки при 3-й степени,
3. 32 при 4-й степени,
4. 16 при 5-й степени, и
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 — Точность вычислений