Разные и самые быстрые способы вычисления синусов и косинусов в Arduino

как правило, справочная таблица > приближение > вычисление. оперативная память > флэш. целое число > фиксированная точка > плавающая точка. предварительный расчет > расчет в реальном времени. зеркальное отображение (синус в косинус или косинус в синус) по сравнению с фактическим вычислением/поиском....

у каждого есть свои плюсы и минусы.

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

изменить: я быстро проверил. при использовании 8-битного целочисленного вывода вычисление 1024 значений sin с помощью справочной таблицы занимает 0,6 мс, а с float — 133 мс, или в 200 раз медленнее.

Двумя основными методами являются математические вычисления (с полиномами) и поисковые таблицы.

Математическая библиотека Arduino (libm, часть avr-libc) использует первую. Он оптимизирован для AVR в том смысле, что он написан на 100% языке ассемблера, и поэтому почти невозможно проследить, что он делает (также нет комментариев). Будьте уверены, что это будет самая оптимизированная реализация чистого числа с плавающей запятой, которую могут придумать умы, намного превосходящие наши.

Однако ключ здесь float. Все, что в Arduino связано с плавающей запятой, будет тяжеловесным по сравнению с чистыми целыми числами, а поскольку вы запрашиваете только целые числа в диапазоне от 0 до 90 градусов, простая таблица поиска — безусловно, самый простой и эффективный метод.

Таблица из 91 значения даст вам все от 0 до 90 включительно. Однако, если вы сделаете эту таблицу значений с плавающей запятой между 0,0 и 1,0, вы все равно столкнетесь с неэффективностью работы с числами с плавающей запятой (предоставленной не так неэффективно, как вычисление sin с числами с плавающей запятой), поэтому сохранение значения с фиксированной точкой вместо этого было бы намного эффективнее.

Это может быть так же просто, как сохранить значение, умноженное на 1000, чтобы у вас было от 0 до 1000, а не от 0,0 до 1,0 (например, sin(30) будет храниться как 500 вместо 0,5). Более эффективным было бы хранить значения как, например, значение Q16, где каждое значение (бит) представляет собой 1/65536 от 1,0. Эти значения Q16 (и связанные с ними Q15, Q1.15 и т. д.) более эффективны для работы, поскольку у вас есть степени двойки, с которыми любят работать компьютеры, а не степени десятых, с которыми они ненавидят работать.

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

Однако возможна комбинация этих двух сценариев. Линейная интерполяция позволит вам получить приблизительное значение угла с плавающей запятой между двумя целыми числами. Это так же просто, как выяснить, как далеко вы находитесь между двумя точками в таблице поиска, и создать средневзвешенное значение на основе этого расстояния двух значений. Например, если у вас 23,6 градуса, вы берете (sintable[23] * (1-0,6)) + (sintable[24] * 0,6). В основном ваша синусоида становится серией дискретных точек, соединенных прямыми линиями. Вы обмениваете точность на скорость.

Вы можете создать пару функций, использующих линейную аппроксимацию для определения sin() и cos() определенного угла.

Я думаю примерно так:

Для каждого я разбил графическое представление функций sin() и cos() на 3 части и сделал линейную аппроксимацию этой части.

В идеале ваша функция должна сначала проверить, находится ли диапазон ангела в диапазоне от 0 до 90.
Затем он будет использовать оператор ifelse, чтобы определить, к какому из трех разделов он принадлежит, а затем выполнить соответствующий линейный расчет (т.е. output = mX + c)

Я искал других людей, которые аппроксимировали cos() и sin(), и наткнулся на этот ответ:

ответ dtb на "Fast Sin/Cos с использованием предварительно вычисленного массива перевода"

По сути, он вычислил, что функция math.sin() из математической библиотеки работает быстрее, чем использование справочной таблицы значений. Но насколько я могу судить, это было вычислено на ПК.

В Arduino включена математическая библиотека, которая может вычислять sin() и cos().

Самым быстрым способом поиска синусов будет таблица поиска. И если вам удобно работать с числами с фиксированной запятой (целые числа, у которых двоичная точка находится не справа от нулевого бита), ваши дальнейшие вычисления с синусом также будут намного быстрее. Затем эта таблица может быть таблицей слов, возможно, во Flash для экономии места в оперативной памяти. Обратите внимание, что в математике вам может понадобиться использовать длинные значения для больших промежуточных результатов.

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

Вот грубый пример на C. По сути, он реализует функцию atan2() из библиотек C, используя CORDIC (т.е. находит угол по двум ортогональным компонентам). Он использует числа с плавающей запятой, но может быть адаптирован для использования с фиксированными точечная арифметика.

/*
 * Простой пример использования алгоритма CORDIC.
 */

#include <stdio.h>
#include <math.h>

#define CORDIC_TABLE_SIZE  16

double cordic_table[CORDIC_TABLE_SIZE];

void init_table(void);
double angle(double I, double Q);

/*
 * Given a sine and cosine component of an
 * angle, compute the angle using the CORIDC
 * algoritm.
 */
double angle(double I, double Q)
{
    int L;
    double K = 1;
    double angle_acc = 0;
    double tmp_I;

    if (I < 0) {
        /* rotate by an initial +/- 90 degrees */
        tmp_I = I;
        if (Q > 0.0) {
            I = Q;           /* subtract 90 degrees */
            Q = -tmp_I;
            angle_acc = -90;
        } else {
            I = -Q;          /* add 90 degrees */
            Q = tmp_I;
            angle_acc = 90;
        }
    } else {
        angle_acc = 0;
    }

    /* rotate using "1 + jK" factors */
    for (L = 0, K = 1; L <= CORDIC_TABLE_SIZE; L++) {
        tmp_I = I;
        if (Q >= 0.0) {
            /* angle is positive: do negative roation */
            I += Q * K;
            Q -= tmp_I * K;
            angle_acc -= cordic_table[L];
        } else {
            /* angle is negative: do positive rotation */
            I -= Q * K;
            Q += tmp_I * K;
            angle_acc += cordic_table[L];
        }
        K /= 2.0;
    }
    return -angle_acc;
}

void init_table(void)
{
    int i;
    double K = 1;

    for (i = 0; i < CORDIC_TABLE_SIZE; i++) {
        cordic_table[i] = 180 * atan(K) / M_PI;
        K /= 2.0;
    }
}
int main(int argc, char **argv)
{
    double I, Q, A, Ar, R, Ac;

    init_table();

    printf("# Angle,    CORDIC Angle,  Error\n");
    for (A = 0; A < 90.0; A += 0.5) {

        Ar = A * M_PI / 180; /* convert to radians for C's sin & cos fn's */

        R = 5;  // Произвольный радиус

        I = R * cos(Ar);
        Q = R * sin(Ar);

        Ac = angle(I, Q);
        printf("%9f, %9f,   %12.4e\n", A, Ac, Ac-A);
    }
    return 0;
}

Но сначала попробуйте родные триггерные функции Arduino — они в любом случае могут быть достаточно быстрыми.

Я немного поигрался с вычислением синусов и косинусов на Arduino использует полиномиальные приближения с фиксированной точкой. Вот мои измерения среднего времени выполнения и наихудшей ошибки по сравнению со стандартными cos() и sin() из avr-libc:

function    max error   cycles   time
-----------------------------------------
cos_fix()   9.53e-5     108.25    6.77 µs
sin_fix()   9.53e-5     110.25    6.89 µs
cos()       2.98e-8     1720.8   107.5 µs
sin()       2.98e-8     1725.1   107.8 µs

Она основана на полиноме 6-й степени, рассчитанном только с 4 умножения. Сами умножения выполняются в сборке, так как я обнаружил, что gcc реализовал их неэффективно. Углы выражаются как uint16_t в единицах 1/65536 оборота, что делает арифметику углов естественной по модулю одного оборота.

Если вы считаете, что это может удовлетворить ваши потребности, вот код: Фиксированная точка тригонометрия. Извините, я так и не перевел эту страницу, которая на французском языке, но вы может понимать уравнения и код (имена переменных, комментарии...) на английском языке.

Изменить: Поскольку сервер, похоже, исчез, вот некоторая информация о приближения, которые я нашел.

Я хотел записать углы в двоичной фиксированной точке, в квадрантах (или, что то же самое, по очереди). И я также хотел использовать даже полиномиальные, так как они более эффективны для вычисления, чем произвольные многочлены. Другими словами, мне нужен полином P() такой, что

cos(π/2 x) ≈ P(x²) для x ∈ [0,1]

Я также требовал, чтобы аппроксимация была точной на обоих концах интервал, чтобы убедиться, что cos(0) = 1 и cos(π/2) = 0. Эти ограничения привели к форме

P(u) = (1 − u)(1 + uQ(u))

где Q() — произвольный полином.

Затем я искал наилучшее решение в зависимости от степени Q() и нашел это:

        Q(u)             │ degree of P(x²) │ max error
─────────────────────────┼─────────────────┼──────────
          0              │         2       │  5.60e-2
       −0.224            │         4       │  9.20e-4
−0.2335216 + 0.0190963 u │         6       │  9.20e-6

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

У меня был вопрос, похожий на ОП. Я хотел создать таблицу LUT для вычисления первого квадранта синусоидальной функции как беззнаковых 16-битных целых чисел, начиная с 0x8000 до 0xffff. И я закончил тем, что написал это для удовольствия и прибыли. Примечание. Это работало бы более эффективно, если бы я использовал операторы «если». Также это не очень точно, но будет достаточно точно для синусоиды в звуковом синтезаторе

void sin_lut_ctor(){

//Составить справочную таблицу для 511 условий функции синуса.
// Подключаем некоторые полиномы, чтобы творить чудеса
// и вы получите приближение для синусов до π/2.
//

//Все синусы, кроме π/2, можно вывести с помощью математики

const uint16_t uLut_d = 0x0200; //максимальная глубина LUT для π/2 термов.
uint16_t uLut_0[uLut_d];        //Сам LUT.
//Поместите 2 выше перед вашим void setup() в качестве глобальных переменных.
//Эти коэффициенты будут работать только для uLut_d = 511.

uint16_t arna_poly_0 = 0x000a; // 11
uint16_t arna_poly_1 = 0x0001; // 1
uint16_t arna_poly_2 = 0x0007; // 7
uint16_t arna_poly_3 = 0x0001; // 1 предварительно вычисленный полином
uint16_t arna_poly_4 = 0x0001; // 1
uint16_t arna_poly_5 = 0x0007; // 7
uint16_t arna_poly_6 = 0x0002; // 2
uint16_t arna_poly_7 = 0x0001; // 1

uint16_t Imm_UI_0 = 0x0001;              // Итератор
uint16_t Imm_UI_1 = 0x007c;              // Инкремент, уменьшающийся во времени

uint16_t Imm_UI_2 = 0x0000;              //
uint16_t Imm_UI_3 = 0x0000;              //
uint16_t Imm_UI_4 = 0x0000;              //
uint16_t Imm_UI_5 = 0x0000;              //
uint16_t Imm_UI_6 = 0x0000;              // Временные переменные
uint16_t Imm_UI_7 = 0x0000;              //
uint16_t Imm_UI_8 = 0x0000;              //
uint16_t Imm_UI_9 = 0x0000;              //
uint16_t Imm_UI_A = 0x0000;
uint16_t Imm_UI_B = 0x0000;

uint16_t Imm_UI_A = uLut_d - 0x0001;     // 510

uLut_0[0x0000] = 0x8000;        // Предположим, что средняя точка равна 32768 (0x8000 hex)
while (Imm_UI_0 < Imm_UI_A) //Построить четверть таблицы синусов
  {
Imm_UI_2++;                                   //Увеличиваем временную переменную на 1

Imm_UI_B = Imm_UI_2 / arna_coeff_0;           // Делим с первым коэффициентом (примечание: целочисленное деление)
Imm_UI_3 += Imm_UI_B;                         //Увеличиваем следующее временное значение, если первое увеличилось до 1-го коэффициента
Imm_UI_1 -= Imm_UI_B;                         // Уменьшаем инкремент, если это так
Imm_UI_2 *= 0x001 - Imm_UI_B;                 // Устанавливаем первую временную переменную обратно в 0

Imm_UI_B = Imm_UI_3 / arna_poly_1;           // Делаем то же самое, что и раньше, со следующим набором временных переменных
Imm_UI_4 += Imm_UI_B;
Imm_UI_1 -= Imm_UI_B;
Imm_UI_3 *= 0x0001 - Imm_UI_B;

Imm_UI_B = Imm_UI_4 / arna_poly_2;           // И снова... и снова... ну вы поняли.
Imm_UI_5 += Imm_UI_B;
Imm_UI_1 -= Imm_UI_B;
Imm_UI_4 *= 0x0001 - Imm_UI_B;

Imm_UI_B = Imm_UI_5 / arna_poly_3;
Imm_UI_6 += Imm_UI_B;
Imm_UI_1 -= Imm_UI_B;
Imm_UI_5 *= 0x0001 - Imm_UI_B;

Imm_UI_B = Imm_UI_6 / arna_poly_4;
Imm_UI_7 += Imm_UI_B;
Imm_UI_1 -= Imm_UI_B;
Imm_UI_6 *= 0x0001 - Imm_UI_B;

Imm_UI_B = Imm_UI_7 / arna_poly_5;
Imm_UI_8 += Imm_UI_B;
Imm_UI_1 -= Imm_UI_B;
Imm_UI_7 *= 0x0001 - Imm_UI_B;

Imm_UI_B = Imm_UI_8 / arna_poly_6;
Imm_UI_9 += Imm_UI_B;
Imm_UI_1 -= Imm_UI_B;
Imm_UI_8 *= 0x0001 - Imm_UI_B;

Imm_UI_B = Imm_UI_9 / arna_poly_7          //в последнем наборе не нужно будет увеличивать следующую переменную, поэтому пропустите шаг, на котором вы хотите ее увеличить.
Imm_UI_1 -= Imm_UI_B;
Imm_UI_9 *= 1 - Imm_UI_B;

uLut_0[Imm_UI_0] = (uLut_0[Imm_UI_0 - 0x0001] + Imm_UI_1); //Устанавливаем текущее значение как предыдущее, увеличенное нашим инкрементом
Imm_UI_0++;              //Увеличиваем итератор
  }   
  uLut_0[Imm_UI_A] = 0xffff; //Наконец, устанавливаем последнее значение на 0xffff

  //И вот оно. Таблица синусов только с одним оператором if (цикл while)
}

Теперь, чтобы вернуть значения, используйте эту функцию. Она принимает значение от 0x0000 до 0x0800 и возвращает соответствующее значение из LUT

uint16_t lu_sin(uint16_t lu_val0)
{
  //Получить значение от 0x0000 до 0x0800. Вернуть соответствующий грех (значение)
  Imm_UI_0 = lu_val0/0x0200; //определение квадранта
  Imm_UI_1 = lu_val0%0x0200; //Получить какое значение
  if (Imm_UI_0 == 0x0000)
  {
    return uLut_0[Imm_UI_1];
  }
  if (Imm_UI_0 == 0x0001)
  {
    return uLut_0[0x01ff - Imm_UI_1];
  }
  if (Imm_UI_0 == 0x0002)
  {
    return 0xffff - uLut_0[Imm_UI_1];
  }
  if (Imm_UI_0 == 0x0003)
  {
    return 0xffff - uLut_0[0x01ff - Imm_UI_1];
  }
}// Здесь я использую операторы if, но аналогично блоку кода выше,
 // можно обойтись без. только с целочисленными делениями и модулями

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

Просто для удовольствия и чтобы доказать, что это можно сделать, я завершил процедуру сборки AVR для вычисления результатов sin(x) в 24 битах (3 байта) с одним битом ошибки. Входной угол указывается в градусах с одной десятичной цифрой, от 000 до 900 (0~90,0) только для первого квадранта. Он использует менее 210 инструкций AVR и работает в среднем за 212 микросекунд, варьируя от 211 мкс (угол = 001) до 213 мкс (угол = 899).

На все это ушло несколько дней, более 10 дней (свободных часов) только размышления о наилучшем алгоритме расчета, учитывая микроконтроллер AVR, отсутствие плавающей запятой, устранение всех возможных делений. Больше времени потребовалось, чтобы сделать правильные повышающие значения для целых чисел, чтобы иметь хорошую точность, необходимую для повышения значений 1e-8 до двоичных целых чисел 2 ^ 28 или более. После того, как все ошибки, виновные в точности и округлении, были найдены, увеличили разрешение их вычислений на дополнительные 2 ^ 8 или 2 ^ 16, были достигнуты наилучшие результаты. Сначала я смоделировал все вычисления в Excel, позаботившись о том, чтобы все значения были как Int(x) или Round(x,0), чтобы точно представлять обработку ядра AVR.

Например, в алгоритме угол должен быть в радианах, а ввод в градусах для удобства пользователя. Чтобы преобразовать градусы в радианы, тривиальная формула: рад=градусы*PI/180, это кажется красивым и простым, но это не так, PI - бесконечное число - если использовать несколько цифр, это создаст ошибки на выходе, деление на 180 требует Манипуляции с битами AVR, поскольку в нем нет команды деления, и, более того, для результата потребуется число с плавающей запятой, поскольку оно включает числа, намного меньшие целого числа 1. Например, радиан 1 ° (градус) равен 0,017453293. Поскольку числа PI и 180 являются константами, почему бы не обратить это в обратное для простого умножения? PI/180 = 0,017453293, умножьте его на 2^32 и получите константу 74961320 (0x0477D1A8), умножьте это число на ваш угол в градусах, скажем, 900 для 90° и сдвиньте его на 4 бита вправо (÷16), чтобы получить 4216574250 (0xFB53D12A), то есть радианы 90° с расширением 2^28, помещаются в 4 байта без единого деления (кроме 4-битного сдвига вправо). В некотором смысле ошибка, заложенная в такой трюк, меньше, чем 2^-27.

Итак, во всех дальнейших расчетах нужно помнить, что это на 2 ^ 28 больше, и позаботиться об этом. Вам нужно разделить результаты на ходу на 16, 256 или даже 65536, чтобы избежать использования ненужных растущих голодных байтов, которые не помогут разрешению. Это была кропотливая работа, просто найти минимальное количество битов в каждом результате вычислений, сохраняя точность результатов около 24 бит. Каждое из нескольких вычислений выполнялось в режиме «попытка/ошибка» со старшими или младшими битами в потоке Excel, наблюдая за общим количеством ошибочных битов в результате на графике, показывающем 0-90° с макросом, запускающим код 900 раз, один раз на десятую градуса. Этот «визуальный» подход к Excel был инструментом, который я создал, и он очень помог найти лучшее решение для каждой отдельной части кода.

Например, округление этого конкретного результата расчета 13248737,51 до 13248738 или просто потеря десятичных знаков "0,51", насколько это повлияет на точность конечного результата для всех 900 входных углов (00,1 ~ 90,0) тестов?

Мне удавалось удерживать животное в пределах 32 бит (4 байта) при каждом вычислении, и в итоге чудом удалось получить точность в пределах 23 бит результата. При проверке всех 3 байтов результата ошибка составляет ±1 младший бит, ожидание.

Пользователь может получить один, два или три байта из результата для собственных требований к точности. Конечно, если достаточно одного байта, я бы порекомендовал использовать одну таблицу sin размером 256 байт и использовать инструкцию AVR 'LPM' для ее захвата.

После того, как последовательность действий в Excel стала гладкой и аккуратной, окончательный перевод из Excel в сборку AVR занял меньше двух часов, как обычно, сначала нужно больше думать, а потом меньше работать.

В то время я смог еще больше сократить использование регистров. Фактический (не окончательный) код использует около 205 инструкций (~ 410 байт), выполняет вычисление sin(x) в среднем за 212 мкс, тактовая частота 16 МГц. На такой скорости он может вычислить 4700+ sin(x) в секунду. Неважно, но он может воспроизводить точную синусоиду до 4700 Гц с точностью и разрешением 23 бита без каких-либо таблиц поиска.

Базовый алгоритм основан на ряде Тейлора для sin(x), но сильно изменен, чтобы соответствовать моим намерениям с учетом микроконтроллера AVR и точности.

Даже несмотря на то, что использование таблицы размером 2700 байт (900 записей * 3 байта) было бы привлекательным с точки зрения скорости, что в этом интересного или обучающего опыта? Конечно, подход CORDIC также рассматривался, может быть позже, здесь речь идет о том, чтобы втиснуть Тейлора в ядро AVR и взять воду из сухой породы.

Интересно, может ли Arduino «sin(78.9°)» выполнять обработку (C++) с точностью 23 бита менее чем за 212 мкс, а необходимый код — менее чем за 205 инструкций. Может быть, если С++ использует CORDIC. Скетчи Arduino могут импортировать ассемблерный код.

Нет смысла публиковать код здесь, позже я отредактирую это сообщение, включив в него веб-ссылку, возможно, в моем блоге по адресу этот URL. Блог в основном на португальском языке.

Это интересное предприятие, не требующее денег, раздвинуло пределы движка AVR почти до 16MIPS на частоте 16 МГц, без инструкции деления, умножение только в 8x8 бит. Это позволяет вычислить sin(x), cos(x) [=sin(900-x)] и tan(x) [=sin(x)/sin(900-x)].

Прежде всего, это помогло моему 63-летнему мозгу отполироваться и смазаться. Когда подростки говорят, что «старики» ничего не смыслят в технологиях, я отвечаю: «Подумайте еще раз, кто, по-вашему, создал основу для всего, чем вы сегодня наслаждаетесь?».

Ура!

Как уже упоминалось, таблицы поиска — это то, что вам нужно, если вам нужна скорость. Недавно я исследовал вычисление триггерных функций на ATtiny85 для использования быстрых векторных средних (в моем случае ветра). Всегда есть компромисс... для меня мне нужно было только угловое разрешение в 1 градус, поэтому таблица поиска из 360 целых чисел (масштабирование от -32767 до 32767, работа только с целыми числами) была лучшим способом. Чтобы получить синус, достаточно ввести индекс 0-359... очень быстро! Некоторые цифры из моих тестов:

Время поиска во FLASH (мкс): 0,99 (таблица сохранена с помощью PROGMEM)

Время поиска в ОЗУ (мкс): 0,69 (таблица в ОЗУ)

Время библиотеки (мкс): 122,31 (с использованием Arduino Lib)

Обратите внимание, что это средние значения по выборке из 360 точек для каждого из них. Тестирование проводилось на нано.