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

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

Бібліотека математики 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) . В основному ваша синусова хвиля стає серією дискретних точок, з'єднаних прямими лініями. Ви торгуєте точністю для швидкості.

Є кілька хороших відповідей тут, але я хотів додати метод, який ще не був згаданий, один дуже добре підходить для обчислення тригонометричних функцій на вбудованих системах, і це техніка 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 з використанням поліноміальних наближень з фіксованою точкою. Ось мої вимірювання середнього часу виконання і найгіршого випадку помилки, порівняння зі стандартним 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 ² ) для 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-бітової реалізації.

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

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

Я маю на увазі щось подібне:

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

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

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

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

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

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

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 повільніше.

Я мав скрутне питання до ОП. Я хотів зробити таблицю 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

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