Арифметика с фиксированной точкой

Формат Q

Число в формате Qm.n — это знаковое целое из 1 + m + n бит, которое интерпретируется как целое, делённое на 2^n. Отсюда:

Перевод: целое = round(значение · 2^n), значение = целое / 2^n. Например, 0,3 в Q15 — round(0,3 · 32 768) = 9830, а обратно 9830 / 32 768 = 0,29999 — погрешность меньше половины шага. Обозначения 1.7, 1.15 и 1.31 использует CMSIS-DSP (в заголовке arm_math_types.h q7_t, q15_t и q31_t описаны как 8-, 16- и 32-битные дробные типы в этих форматах); технический отчёт TR 18037 записывает те же форматы как s.7, s.15 и s.31, где s — знаковый бит.

Формат выбирают по диапазону величины: коэффициенты цифрового фильтра, лежащие между −1 и 1, удобно держать в Q15, а температуру от −40 до 125 °C с шагом около 0,004 °C — в Q8.8.

Сложение и вычитание

Числа одного формата складываются как обычные целые: (a + b) / 2^n = a / 2^n + b / 2^n. Числа разных форматов сначала приводят к одному сдвигом. Результат может выйти за диапазон: 0,75 + 0,5 = 1,25 в Q15 не помещается, и при сложении в int16_t получится отрицательное число. Есть два выхода: считать в более широком типе и ограничить результат (насыщение) или заранее масштабировать данные так, чтобы переполнение было невозможно.

#include <stdint.h>

/* Q15 addition with saturation */
static int16_t q15_add_sat(int16_t a, int16_t b)
{
    int32_t s = (int32_t)a + b;
    if (s > INT16_MAX) return INT16_MAX;    /* 0x7FFF, +0.99997 */
    if (s < INT16_MIN) return INT16_MIN;    /* 0x8000, -1.0 */
    return (int16_t)s;
}

Умножение

Произведение двух чисел Qm.n имеет 2n дробных бит: (a / 2^n) · (b / 2^n) = (a · b) / 2^2n. Чтобы вернуться к формату Qm.n, произведение сдвигают вправо на n бит. Перемножать нужно в удвоенной разрядности — для Q15 в 32 битах, — иначе старшие биты потеряются:

/* Q15 multiplication with rounding and saturation */
static int16_t q15_mul(int16_t a, int16_t b)
{
    int32_t p = (int32_t)a * b;      /* Q30 */
    p += 1L << 14;                   /* round to nearest */
    p >>= 15;                        /* back to Q15 */
    if (p > INT16_MAX) p = INT16_MAX;   /* only -1.0 * -1.0 overflows */
    return (int16_t)p;
}

Приведение (int32_t)a обязательно на 8-битном AVR: там int 16-битный, и a * b без него вычислялось бы в 16 битах. Сдвиг вправо отрицательного произведения стандарт C оставляет на усмотрение реализации; GCC делает его арифметическим, с сохранением знака. Без прибавления 1 << 14 сдвиг округляет вниз, и при многократных умножениях ошибка накапливается в одну сторону.

Единственный случай переполнения при умножении Q15 — (−1) · (−1) = +1, которое не представимо: −32 768 · −32 768 = 2^30, после сдвига 32 768, а максимум int16_t — 32 767. Поэтому в конце стоит ограничение.

Деление

Деление обратно умножению: делимое сначала сдвигают влево на n бит в удвоенной разрядности, затем делят на делитель: (int16_t)(((int32_t)a << 15) / b). Результат помещается в Q15, только если |a| меньше |b|, а деление на ноль нужно исключить заранее. Деление — самая медленная операция на МК без аппаратного делителя (у AVR и Cortex-M0 его нет), поэтому в циклах деление на константу заменяют умножением на заранее вычисленную обратную величину.

Встроенные типы: TR 18037 и GCC

Технический отчёт ISO/IEC TR 18037 «Embedded C» вводит в язык типы _Fract (только дробная часть, значения от −1 до 1) и _Accum (с целой частью), их беззнаковые и насыщающие варианты _Sat. Минимальные форматы — s.7, s.15 и s.23 для short _Fract, _Fract и long _Fract и s4.7, s4.15, s4.23 для соответствующих _Accum (раздел 4.1). Насыщение гарантировано только для типов _Sat или при включённых прагмах FX_FRACT_OVERFLOW и FX_ACCUM_OVERFLOW, а без них переполнение не определено (раздел 4.1.3).

GCC поддерживает эти типы как расширение не на всех платформах. В avr-gcc они есть начиная с версии 4.8, поддержка неполная, а форматы такие:

/* avr-gcc: _Fract is the same Q15 format, handled by the compiler */
_Fract gain = 0.75r;
_Fract y = gain * x;         /* x is a _Fract: a Q15 multiply done by the compiler */

Беззнаковые варианты получают префикс u в суффиксе: 1.0uk — unsigned _Accum. Такой код короче, но привязан к компилятору и платформе, поэтому в переносимых библиотеках, например в CMSIS-DSP для ARM, фиксированная точка реализована на целых типах q7_t, q15_t и q31_t.

Типичные ошибки

• Умножение без расширения разрядности — a * b двух int16_t на AVR считается в 16 битах и теряет старшую часть.
• Сдвиг до умножения вместо сдвига после: (a >> 8) * (b >> 7) сохраняет диапазон, но выбрасывает младшие биты обоих сомножителей.
• Смешение форматов — сложение Q15 с Q8.8 без приведения даёт бессмыслицу без всякой ошибки компиляции.
• Переполнение промежуточной суммы в цикле фильтра — аккумулятор берут шире: сумму произведений Q15 копят в 32 или 64 битах и только в конце возвращают в Q15.
• Округление только вниз — систематическая ошибка, которая накапливается в рекурсивных фильтрах и интеграторах ПИД-регулятора.

Таблицы синусов и других функций в формате с фиксированной точкой удобно генерировать калькулятором таблиц, а особенности целочисленной арифметики C на микроконтроллерах разобраны в статье «Язык C для микроконтроллеров».