Different ways (and the fastest) to compute sines (and cosines) in Arduino

Question

I am using an Arduino Uno board to compute the angles of my system (robotic arm). The angles are actually 10 bit values (0 to 1023) from the ADC, using the full range of the ADC. I am only going to be operating in the 1st quadrant (0 to 90 deg), where both sines and cosines are positive, so there is no problem with negative numbers. My doubts can be expressed in 3 questions:

1. What are the different ways to compute these trigonometric functions on Arduino?

2. What is the fastest way to do the same?

3. There are the sin() and cos() functions in the Arduino IDE, but how does the Arduino actually calculate them (as in do they use look-up tables, or approximations etc.)? They seem like an obvious solution, but I would like to know their actual implementation before I try them out.

PS: I am open to both standard coding on the Arduino IDE and assembly coding, as well as any other options not mentioned. Also I have no problems with errors and approximations, which are inevitable for a digital system; however if possible it would be good to mention the extent of possible errors

Answer 1

The two basic methods are mathematical calculation (with polynomials) and lookup tables.

The Arduino's math library (libm, part of avr-libc) uses the former. It is optimised for the AVR in that it is written with 100% assembly language, and as such is almost impossible to follow what it's doing (there's zero comments as well). Rest assured though it will be the most optimised pure-float implementation brains far superior to ours could come up with.

However the key there is float. Anything on the Arduino involving floating point is going to be heavyweight by comparison to pure integer, and since you are only requesting integers between 0 and 90 degrees a simple lookup table is by far the simplest and most efficient method.

A table of 91 values will give you everything from 0 to 90 inclusive. However if you make that a table of floating point values between 0.0 and 1.0 you still then have the inefficiency of working with floats (granted not as inefficient as calculating sin with floats), so storing a fixed point value instead would be far more efficient.

That may be as simple as storing the value multiplied by 1000, so you have between 0 and 1000 instead of between 0.0 and 1.0 (for instance sin(30) would be stored as 500 instead of 0.5). More efficient would be to store the values as, for instance, a Q16 value where each value (bit) represents 1/65536th of 1.0. These Q16 values (and the related Q15, Q1.15, etc) are more efficient to work with since you have powers-of-two which computers love to work with instead of powers-of-ten which they hate working with.

Don't forget as well that the sin() function expects radians, so you first have to convert your integer degrees into a floating point radians value, making the use of sin() even more inefficient compared to a lookup table that can work directly with the integer degrees value.

A combination of the two scenarios, though, is possible. Linear interpolation will allow you to get an approximation of a floating point angle between two integers. It's as simple as working out how far between two points in the lookup table you are and creating a weighted average based on that distance of the two values. For instance if you are at 23.6 degrees you take (sintable[23] * (1-0.6)) + (sintable[24] * 0.6). Basically your sine wave becomes a series of discrete points joined together by straight lines. You trade accuracy for speed.

Answer 2

There are some good answers here but I wanted to add a method which hasn't been mentioned yet, one very well suited to computing trigonometric functions on embedded systems, and that's the CORDIC technique Wiki Entry Here It can compute trig functions using only shifts and adds and a small look-up table.

Here's a crude example in C. In effect, it implements the C libraries' atan2() function using CORDIC (i.e. find an angle given two orthogonal components.) It uses floating point, but can be adapted for use with fixed-point arithmetic.

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

#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;  // 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;
}

But try the native Arduino trig functions first - they might be fast enough anyway.

Answer 3

I have been playing a bit with computing sines and cosines on the Arduino using fixed-point polynomial approximations. Here are my measurements of average execution time and worst case error, compared with the standard cos() and sin() from 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.

If you think this can suit your bill, here is the code: Fixed-point trigonometry. Sorry, I still did not translate this page, which is in French, but you can understand the equations, and the code (variable names, comments...) is in English.

---

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

I wanted to write angles in binary fixed-point, in units of quadrants (or, equivalently, in turns). And I also wanted to use an even polynomial, as these are more efficient to compute than arbitrary polynomials. In other words, I wanted a polynomial P() such that

cos(π/2 x) ≈ P(x²) for 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))

where Q() is an arbitrary polynomial.

Next, I searched for the best solution as a function of the degree of Q() and found this:

        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 choice among the solutions above is a speed/accuracy trade-off. The third solution gives more accuracy than achievable with 16-bits, and it's the one I chose for the 16-bit implementation.

Answer 4

You could create a couple of functions that uses linear approximation to determine the sin() and cos() of a particular angle.

I am thinking something like this:

For each I have broken the graphical representation of sin() and cos() into 3 sections and have done a linear approximation of that section.

Your function would ideally first check that the range of the angel is between 0 and 90.
Then it would use an ifelse statement to determine what of the 3 sections it belongs and then does the corresponding linear calculation (ie output = mX + c)

Answer 5

I looked for other people that had approximated cos() and sin() and I came across this answer:

dtb's answer to "Fast Sin/Cos using a pre computed translation array"

Basically he computed that the math.sin() function from the math library was faster than using a look-up table of values. But from what I can tell, this was computed on a PC.

Arduino has a math library included that can calculate sin() and cos().

Answer 6

A lookup table will be the fastest way to find sines. And if you're comfortable computing with fixed-point numbers (integers whose binary-point is somewhere other than to the right of bit-0), your further calculations with the sines will be much faster as well. That table can then be a table of words, possibly in Flash to save RAM space. Note that in your math you may need to use longs for large intermediate results.

Answer 7

Just for the fun of it, and to prove it can be done, I finished an AVR assembly routine to calculate sin(x) results in 24 bits (3 bytes) with one bit of error. The input angle is in degrees with one decimal digit, from 000 to 900 (0~90.0) for the first quadrant only. It uses less than 210 AVR instructions and runs on average of 212 microseconds, varying from 211us (angle=001) to 213us (angle=899).

It took several days to do it all, more than 10 days (free hours) just thinking about the best algorithm for the calculation, considering the AVR microcontroller, no floating point, eliminating all possible divisions. What took more time was to make the right step-up values for integers, to have good precision it needs to step-up values of 1e-8 to binary integers 2^28 or more. Once all the errors culprits of precision and rounding up were found, increased their calculation resolution by extra 2^8 or 2^16, the best results were met. I first simulated all the calculations on Excel taking care of having all values as Int(x) or Round(x,0) to represent exactly the AVR core processing.

For example, in the algorithm the angle must be in Radians, the input is in Degrees to facilitate for the user. To convert Degrees to Radians the trivial formula is rad=degrees*PI/180, it seems nice and easy, but it is not, PI is an infinite number - if using few digits it will create errors at the output, division by 180 requires AVR bit manipulation since it has no division instruction, and more than that, the result would require floating point since involves numbers far below integer 1. For example, Radian of 1° (degree) is 0.017453293. Since PI and 180 are constants, why not reverse this thing for simple multiplication? PI/180 = 0.017453293, multiply it by 2^32 and it results as a constant 74961320 (0x0477D1A8), multiply this number by your angle in degrees, lets say 900 for 90° and shift it 4 bits right (÷16) to obtain 4216574250 (0xFB53D12A), that is the radians of the 90° with 2^28 expansion, fit in 4 bytes, without a single division (except the 4 bit shift right). In a way, the error included in such trick is smaller than 2^-27.

So, all further calculations need to remember it is 2^28 higher and taken care of it. You need to divide the on-the-go results by 16, 256 or even 65536 just to avoid it use unnecessary growing hunger bytes that would not help resolution. That was a painstaking job, just finding the minimum quantity of bits in each calculation results, keeping the results precision around 24 bits. Each one of the several calculations where done in try/error with higher or lower bits count in the Excel flow, watching the overall quantity of error bits at the result in a graph showing 0-90° with a macro running the code 900 times, once per tenth of a degree. That "visual" Excel approach was a tool I created, helped a lot to find the best solution for every single part of the code.

For example, rounding up this particular calculation result 13248737.51 to 13248738 or just lose the "0.51" decimals, how much it will affect the final result precision for all the 900 input angles (00.1 ~ 90.0) tests?

I was able to keep the animal contained within 32 bits (4 bytes) on every calculation, and ended up with the magic to obtain precision within 23 bits of the result. When checking the whole 3 bytes of the result, the error is ±1 LSB, outstanding.

The user may grab one, two or three bytes from the result for its own requirements of precision. Of course, if just one byte is enough I would recommend to use a single 256 bytes sin table and use AVR 'LPM' instruction to grab it.

Once I had the Excel sequence running smooth and neat, the final translation from Excel to AVR assembly took less than 2 hours, as usual you should think more first, work less later.

At that time I was able to squeeze yet more and reduce registers usage. The actual (not final) code uses around 205 instructions (~410 bytes), runs a sin(x) calculation in average of 212us, clock at 16MHz. At that speed it can calculate 4700+ sin(x) per second. Not being important, but it can run a precise sinewave up to 4700Hz with 23 bits of precision and resolution, without any lookup tables.

The base algorithm is based on Taylor series for sin(x), but modified a lot to fit my intentions with the AVR microcontroller and precision in mind.

Even that using a 2700 bytes table (900 entries * 3 bytes) would be speed attractive, what is the fun or learning experience on that? Of course, CORDIC approach was also considered, maybe later, the point here is to squeeze Taylor into the AVR core and take water from a dry rock.

I wonder if Arduino "sin(78.9°)" can run Processing (C++) with 23 bits of precision in less than 212us and the necessary code smaller than 205 instructions. May be if C++ uses CORDIC. Arduino sketches can import assembly code.

Makes no sense to post the code here, later I will edit this post to include a weblink to it, possibly on my blog at this url. The blog is mostly in portuguese.

This hobby-no-money venture was interesting, pushing the limits of the AVR engine of almost 16MIPS at 16MHz, without division instruction, multiplication only in 8x8 bits. It allows to calculate sin(x), cos(x) [=sin(900-x)] and tan(x) [=sin(x)/sin(900-x)].

Above it all, this helped to keep my 63 years old brain polished and oiled. When teenagers say the 'old people' know sh*nothing about technology, I answer "think again, who do you think created the bases for everything you enjoy today?".

Cheers

Answer 8

I had a simillar question to OP. I wanted to make a LUT table to calculate the first quadrant of the sine function as unsigned 16 bit integers starting from 0x8000 to 0xffff. And i ended up writing this for fun and profit. Note: This would work more efficiently if i used 'if' statements. Also it's not very accurate, but would be accurate enough for a sine wave in a sound synthesiser

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)
}

Now to get back the values, use this function.It accepts a value from 0x0000 to 0x0800 and returns the appropriate value from the 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

Remember, this is not the most efficient approach to this task, i just couldn't figure out how to make taylor series to give out results in the appropriate range.

Answer 9

As others have mentioned lookup tables are the way to go if you want speed. I've recently been investigating computation of trig functions on an ATtiny85 for use of fast vector averages (wind in my case). There's always trade off's...for me I only needed 1 deg angular resolution so a lookup table of 360 int's (scaling -32767 to 32767, only working with int's) was the best way to go. Retrieving the sine is just a matter of supplying an index 0-359...so very fast! Some numbers from my tests:

FLASH lookup time (us): 0.99 (table stored using PROGMEM)

RAM lookup time (us): 0.69 (table in RAM)

Lib time (us): 122.31 (Using Arduino Lib)

Note these are averages across a 360 point sample for each. Testing was done on a nano.

Answer 10

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....

each has its pluses and minuses.

you can make all sorts of combinations to see which works best for your application.

edit: I did a quick checking. using 8-bit integer output, calculating 1024 sin values with look-up table takes 0.6ms, and 133ms with floaters, or 200x slower.