Різні способи (і найшвидші) для обчислення синусів (і косинусів) в Arduino

Я використовую плату Arduino Uno для обчислення кутів моєї системи (робототехніка). Кути фактично 10 бітових значень (від 0 до 1023) з АЦП, використовуючи весь діапазон АЦП. Я буду працювати тільки в 1-му квадранті (від 0 до 90 градусів), де і синуси, і косинуси є позитивними, тому немає жодних проблем з негативними числами. Мої сумніви можна висловити в трьох питаннях:

  1. Які різні способи обчислення цих тригонометричних функцій на Arduino?

  2. Який найшвидший спосіб зробити так само?

  3. У IDE Arduino є функції sin() та cos (), але як насправді їх обчислює Arduino (як у використанні таблиць пошуку або наближень тощо)? Вони виглядають очевидним рішенням, але я хотів би дізнатися про їхню реальну реалізацію, перш ніж спробувати їх.

PS: Я відкритий як для стандартного кодування на IDU Arduino, так і для кодування збірок, а також будь-які інші не згадані опції. Також у мене немає проблем з помилками і наближеннями, які неминучі для цифрової системи; однак, якщо це можливо, було б добре згадати ступінь можливих помилок

8
Я припускаю, що ви хочете працювати в градусах. Чи хочете ви ввести цілі або десяткові числа для кута?
додано Автор Yoni Baciu, джерело
Ви б добре з приблизними значеннями?
додано Автор Yoni Baciu, джерело
Не могли б Ви визначити кількість вимог щодо точності? Наближення cos (π/2x) ≈ 1 − x² має максимальну похибку 5,6e-2. А (1-x²) (1−0.224x²), що коштує 3 множення, добре дорівнює 9.20e-4.
додано Автор Sprogz, джерело
Для всього 90 (цілих) градусів найдосконаліша і найефективніша буде таблиця 90-вхідного пошуку. Насправді для повних 360 градусів можна використовувати 90-вхідну таблицю пошуку. Просто прочитайте його назад за 90-179 і переверніть його на 180-269. Зробіть обидва для 270-359.
додано Автор Majenko, джерело
@EdgarBonet Вибачте за пізню відповідь. Я не маю ніякої поточної кількісної фіксованої точності. Я просто хочу знати всі можливі варіанти
додано Автор Ken Arnold, джерело
Ступені так. Я думаю, що було б легше писати код і тестувати, якщо ми використовуємо цілі числа, так що я пішов би з цим. Я буду надавати більш чітку інформацію про зміни
додано Автор Ken Arnold, джерело
Так, насправді, але я хотів би знати, наскільки поширені помилки різних методів. Це не точний продукт, а мій проект. Насправді наближення неминучі для майже будь-якої (якщо не будь-якої) цифрової системи, що реалізує математичну функцію
додано Автор Ken Arnold, джерело

8 Відповіді

Двома основними методами є математичний розрахунок (з поліномами) і таблиці пошуку.

Бібліотека математики 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/65536th з 1.0. Ці значення Q16 (і відповідні Q15, Q1.15 і т.д.) є більш ефективними для роботи, оскільки у вас є потужність двох, з якою комп'ютери люблять працювати, а не з десятьма силами, з якими вони ненавидять працювати.

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

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

8
додано
Я написав бібліотеку деякий час назад, що використовував поліном Тейлора для sin/cos, який був швидше, ніж бібліотека. Враховуючи, я використовував плаваючою точкою радіани як вхід для обох.
додано Автор tuskiomi, джерело

Є кілька хороших відповідей тут, але я хотів додати метод, який ще не був згаданий, один дуже добре підходить для обчислення тригонометричних функцій на вбудованих системах, і це техніка CORDIC Запис у вікі тут

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

/*
 * Simple example of using the CORDIC algorithm.
 */

#include 
#include 

#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; //Arbitrary radius

        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 - вони можуть бути досить швидкими.

5
додано
Я взяв аналогічний підхід у минулому, на stm8. вона виконує два кроки: 1) обчислює sin (x) і cos (x) від sin (2x), а потім 2) обчислює sin (x +/- x/2) від sin (x), sin (x/2) , cos (x) і cos (x/2) -> через ітерацію можна наблизитися до своєї мети. у моєму випадку, я почав з 45 градусів (0,707), і працював мій шлях до мети. вона значно повільніше, ніж стандартна функція IAR sin ().
додано Автор dannyf, джерело

Я грав трохи з обчислювальними синусами і косинусами на 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

It's based on a 6th degree polynomial computed with only 4 multiplications. The multiplications themselves are done in assembly, as I found that gcc implemented them inefficiently. The angles are expressed as uint16_t in units of 1/65536 of a revolution, which makes the arithmetic of angles naturally work modulo one revolution.

Якщо ви вважаєте, що це може відповідати вашому рахунку, ось код: Фіксована точка тригонометрія . На жаль, я ще не переклав цю сторінку, яка є французькою, але ви може зрозуміти рівняння, а код (імена змінних, коментарі ...) англійською мовою.


Edit: Since the server seems to have vanished, here is some info on the approximations I found.

Я хотів записати кути в двійковій фіксованій точці, в одиницях квадрантів (або, що еквівалентно, по черзі). І я також хотів використати парне поліном, оскільки вони є більш ефективними для обчислення, ніж довільні поліноми. Іншими словами, я хотів поліном P() такий, що

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

I also required the approximation to be exact at both ends of the interval, to ensure that cos(0) = 1 and cos(π/2) = 0. These constraints led to the form

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

Вибір серед вищезазначених рішень є компромісом швидкості/точності. The третє рішення дає більшу точність, ніж досяжну при 16-бітовому, і це той, який я вибрав для 16-бітової реалізації.

5
додано
@TLW: додав, що до відповіді.
додано Автор Sprogz, джерело
@TLW: Я вимагав, щоб у нього були деякі "приємні" властивості (наприклад, cos (0) = 1), які обмежені формою (1 − x²) (1 + x²Q (x²)), де Q (u) довільний поліном (це пояснюється на сторінці). Я взяв першу ступінь Q (тільки 2 коефіцієнти), знайшов приблизні коефіцієнти по придатності, потім вручну налаштував оптимізацію методом проб і помилок.
додано Автор Sprogz, джерело
Це дивовижно, @Edgar.
додано Автор SDsolar, джерело
Що ви зробили, щоб знайти поліном?
додано Автор ThomasX, джерело
@EdgarBonet - цікавий. Зауважте, що ця сторінка не завантажується для мене, хоча кешована. Не могли б ви додати поліном, який використовується для цієї відповіді?
додано Автор ThomasX, джерело

Таблиця пошуку буде найшвидшим способом пошуку синусів. І якщо вам зручно обчислюватися з номерами з фіксованою точкою (цілі числа, двійкова точка яких знаходиться десь інше, ніж праворуч від біта-0), ваші подальші розрахунки з синусами будуть також набагато швидше. Ця таблиця може бути таблицею слів, можливо, у Flash для збереження простору оперативної пам'яті. Зауважте, що у вашій математиці може знадобитися використовувати довгі проміжні результати.

2
додано

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

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

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

2
додано
Не обов'язково. Ви могли б мати його, так що вихід масштабується між 0-100 замість 0-1. Таким чином ви маєте справу з цілими числами, а не з плаваючою точкою. Примітка: 100 було довільним. Немає жодних причин, щоб ви не могли масштабувати вихід між 0-128 або 0-512 або 0-1000 або 0-1024. Використовуючи кратне 2, вам потрібно лише зробити правильні зміни, щоб масштабувати результат назад.
додано Автор Yoni Baciu, джерело
Досить розумний, @sa_leinad. Проголосувати. Я пам'ятаю, що це робилося при роботі з зміщенням транзисторів.
додано Автор SDsolar, джерело
Чи не буде це пов'язано з множенням з плаваючою точкою?
додано Автор Ken Arnold, джерело

Я шукав інших людей, які наблизили cos() і sin (), і я натрапив на цю відповідь:

Відповідь dtb на "Швидкий гріх/Cos з використанням попередньо обчисленого масиву перекладу"

В основному він обчислив, що функція math.sin() з бібліотеки математики була швидшою, ніж використання таблиці пошуку значень. Але з того, що я можу сказати, це було обчислено на ПК.

У Arduino є бібліотека математики , яка може обчислювати sin() та cos ().

2
додано
Комп'ютери мають вбудовані FPU, які роблять його швидким. Arduino не, і це робить його повільним.
додано Автор Majenko, джерело
Відповідь також для C#, який робить такі речі, як перевірка меж масиву.
додано Автор Longdaysjourneyintocode, джерело

generally, look-up table > approximation -> calculation. ram > flash. integer > fixed point > floating point. pre-calclation > real time calculation. mirroring (sine to cosine or cosine to sine) vs. actual calculation/look-up....

у кожного є свої плюси і мінуси.

Ви можете зробити всілякі комбінації, щоб побачити, що найкраще підходить для вашої програми.

edit: Я зробив швидку перевірку. з використанням 8-бітового цілочисельного виводу, обчислення 1024 значень гріха з таблицею пошуку займає 0.6ms, а 133ms з floaters, або 200x повільніше.

1
додано

Я мав скрутне питання до ОП. Я хотів зробити таблицю LUT для обчислення першого квадранта функції синуса як беззнакових цілих чисел, починаючи з 0x8000 до 0xffff. Та я закінчився писати це для утіхи та прибутку. Примітка: Це буде працювати більш ефективно, якщо я використовую "якщо" заяви. Крім того, він не дуже точний, але був би достатньо точним для синусоїди у звуковому синтезаторі

void sin_lut_ctor(){

//Make a Look Up Table for 511 terms of the sine function.
//Plugin in some polynomials to do some magic
//and you get an aproximation for sines up to π/2.
//

//All sines yonder π/2 can be derived with math

const uint16_t uLut_d = 0x0200; //maximum LUT depth for π/2 terms. 
uint16_t uLut_0[uLut_d];        //The LUT itself.
//Put the 2 above before your void setup() as global variables.
//This coefficients will only work for 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   Precalculated Polynomials
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;             // Itterator
uint16_t Imm_UI_1 = 0x007c;             // An incrementor that decreases in time

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;             // Temporary variables
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;        //Assume that the middle point is 32768 (0x8000 hex)
while (Imm_UI_0 < Imm_UI_A) //Construct a quarter of the sine table
  {
Imm_UI_2++;                                   //Increase temporary variable by 1

Imm_UI_B = Imm_UI_2/arna_coeff_0;           //Divide it with the first coefficient (note: integer division)
Imm_UI_3 += Imm_UI_B;                         //Increase the next temporary value if the first one has increased up to the 1st coefficient
Imm_UI_1 -= Imm_UI_B;                         //Decrease the incrementor if this is the case
Imm_UI_2 *= 0x001 - Imm_UI_B;                 //Set the first temporary variable back to 0

Imm_UI_B = Imm_UI_3/arna_poly_1;           //Do the same thing as before with the next set of temporary variables
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;           //And again... and again... you get the idea.
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          //the last set won't need to increment a next variable so skip the step where you would increase it.
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); //Set the current value as the previous one increased by our incrementor
Imm_UI_0++;              //Increase the itterator
  }   
  uLut_0[Imm_UI_A] = 0xffff; //Lastly, set the last value to 0xffff

  //And there you have it. A sine table with only one if statement (a while loop)
}

Тепер, щоб повернути значення, скористайтеся цією функцією. Приймає значення від 0x0000 до 0x0800 і повертає відповідне значення з LUT

uint16_t lu_sin(uint16_t lu_val0)
{
  //Get a value from 0x0000 to 0x0800. Return an appropriate sin(value)
  Imm_UI_0 = lu_val0/0x0200; //determine quadrant
  Imm_UI_1 = lu_val0%0x0200; //Get which value
  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];
  }
}// I'm using if statements here but similarly to the above code block, 
 //you can do without. just with integer divisions and modulos

Пам'ятайте, що це не найефективніший підхід до цього завдання, я просто не міг зрозуміти, як зробити ряди Тейлора, щоб дати результати у відповідному діапазоні.

1
додано
Ваш код не компілюється: Imm_UI_A оголошується двічі, ; і деякі декларації змінних відсутні, а uLut_0 має бути глобальним. Маючи необхідні виправлення, lu_sin() є швидким (від 27 до 42 циклів ЦП), але дуже неточним (максимальна помилка ≈ 5.04e-2). Я не можу зрозуміти, в чому полягають ці “поліноми Арнадата”: це здається досить важким обчисленням, але результат майже так само поганий, як просте квадратичне наближення. Метод також має величезну вартість пам'яті. Було б краще обчислити таблицю на вашому ПК і покласти її в вихідний код як масив PROGMEM .
додано Автор Sprogz, джерело