How can a microcontroller run like flying? An efficient way to implement mathematical functions!

Hello everyone, I am Xiaomai. Today I will share with you how to implement the trigonometric function algorithm on a microcontroller with tight resources and low computing power.

I have previously published an IQMath article about , which is ti a mathematical operation library on the company's platform, which encapsulates many efficient mathematical operation methods.

For example, when using fixed-point arithmetic on a fixed-point processor that does not have a floating-point arithmetic unit, I have written an article in Q format before, briefly introducing this knowledge.

So here comes the question. There is a reader friend whose hardware platform cannot be used IQMath , but he wants to perform some trigonometric function calculations. So how can he implement it himself?

Let’s briefly introduce the overall idea, because the resources of the hardware platform are relatively tight;

• RAM is relatively small;
• ROM is relatively small;
• CPU processing speed is slower;

Therefore, the more common method here is to exchange space for time and store the value of in the array in advance sin . cos When needed, the specific data can be obtained by accessing the array. This is the look-up table method we often mention .

Let’s introduce it in detail below.

sine table

The expression of this sine function is like this,

The details are shown in the figure below;

First, let’s briefly analyze this waveform:

• Within the blue box is a full cycle waveform;
• Within the red box is a quarter-cycle waveform;

In fact, it is not difficult to find that as long as we express the data of this quarter of the waveform, the remaining waveforms can be expressed through conversion.

This greatly saves the space required by the lookup table method.

Below we will introduce how to implement it;

First of all, we have to figure out one point, which is the dimension, and do it in the form of normalization.

• The range of y is [-1, 1] ;
• The range of x is [0, 2π] , of course, the range of x [-π, π] is no problem, and we will continue to introduce it below;

In actual programs, we cannot do this. We expect to access these values ​​​​through integer types, so we have to do a few things:

• Try to avoid using floating point operations;
• Try to avoid division;
• Try to avoid multiplication;

Therefore, it is necessary to understand the Q format first and use left shift and right shift to replace multiplication and division to improve operation efficiency;

For the X-axis data, it can be [0, 2π] subdivided into 128, 256, 512 or 1024, etc.;

Here we first subdivide it into 1024 equal parts. As mentioned before, we only need to select the content of the first quarter period;

    #define POINT_NUM  256
 #define PI          3.141592f
 for (int i = 0; i < POINT_NUM; i++) {

        printf( "[%03d:%1.4f]\t", i , (sin(i*PI/2 / POINT_NUM)));
        if((i+1) % 8 == 0){
            printf("\r\n");
        }
    }

The printed output is as follows:

Here we can simply take a few special points to verify and find that the overall performance is acceptable;

The next step is to convert the floating point data y into Q1.15 format.

#define POINT_NUM  256
#define PI          3.141592f
 printf("sin=============================================\r\n");
    for (int i = 0; i < POINT_NUM; i++) {

        printf("[ %d:\t0x%04X ]", i, (int16_t)(sin(i*PI/2 / POINT_NUM) * 32768));
        if((i+1) % 8 == 0){
            printf("\r\n");
        }
    }

The final output is as follows;

Source code part

The following part of the code is an example of reference ST . mcsdk Below we will analyze the role of each part in turn. The overall code is as follows;

#define SIN_COS_TABLE {\
    0x0000,0x00C9,0x0192,0x025B,0x0324,0x03ED,0x04B6,0x057F,\
    0x0648,0x0711,0x07D9,0x08A2,0x096A,0x0A33,0x0AFB,0x0BC4,\
    0x0C8C,0x0D54,0x0E1C,0x0EE3,0x0FAB,0x1072,0x113A,0x1201,\
    0x12C8,0x138F,0x1455,0x151C,0x15E2,0x16A8,0x176E,0x1833,\
    0x18F9,0x19BE,0x1A82,0x1B47,0x1C0B,0x1CCF,0x1D93,0x1E57,\
    0x1F1A,0x1FDD,0x209F,0x2161,0x2223,0x22E5,0x23A6,0x2467,\
    0x2528,0x25E8,0x26A8,0x2767,0x2826,0x28E5,0x29A3,0x2A61,\
    0x2B1F,0x2BDC,0x2C99,0x2D55,0x2E11,0x2ECC,0x2F87,0x3041,\
    0x30FB,0x31B5,0x326E,0x3326,0x33DF,0x3496,0x354D,0x3604,\
    0x36BA,0x376F,0x3824,0x38D9,0x398C,0x3A40,0x3AF2,0x3BA5,\
    0x3C56,0x3D07,0x3DB8,0x3E68,0x3F17,0x3FC5,0x4073,0x4121,\
    0x41CE,0x427A,0x4325,0x43D0,0x447A,0x4524,0x45CD,0x4675,\
    0x471C,0x47C3,0x4869,0x490F,0x49B4,0x4A58,0x4AFB,0x4B9D,\
    0x4C3F,0x4CE0,0x4D81,0x4E20,0x4EBF,0x4F5D,0x4FFB,0x5097,\
    0x5133,0x51CE,0x5268,0x5302,0x539B,0x5432,0x54C9,0x5560,\
    0x55F5,0x568A,0x571D,0x57B0,0x5842,0x58D3,0x5964,0x59F3,\
    0x5A82,0x5B0F,0x5B9C,0x5C28,0x5CB3,0x5D3E,0x5DC7,0x5E4F,\
    0x5ED7,0x5F5D,0x5FE3,0x6068,0x60EB,0x616E,0x61F0,0x6271,\
    0x62F1,0x6370,0x63EE,0x646C,0x64E8,0x6563,0x65DD,0x6656,\
    0x66CF,0x6746,0x67BC,0x6832,0x68A6,0x6919,0x698B,0x69FD,\
    0x6A6D,0x6ADC,0x6B4A,0x6BB7,0x6C23,0x6C8E,0x6CF8,0x6D61,\
    0x6DC9,0x6E30,0x6E96,0x6EFB,0x6F5E,0x6FC1,0x7022,0x7083,\
    0x70E2,0x7140,0x719D,0x71F9,0x7254,0x72AE,0x7307,0x735E,\
    0x73B5,0x740A,0x745F,0x74B2,0x7504,0x7555,0x75A5,0x75F3,\
    0x7641,0x768D,0x76D8,0x7722,0x776B,0x77B3,0x77FA,0x783F,\
    0x7884,0x78C7,0x7909,0x794A,0x7989,0x79C8,0x7A05,0x7A41,\
    0x7A7C,0x7AB6,0x7AEE,0x7B26,0x7B5C,0x7B91,0x7BC5,0x7BF8,\
    0x7C29,0x7C59,0x7C88,0x7CB6,0x7CE3,0x7D0E,0x7D39,0x7D62,\
    0x7D89,0x7DB0,0x7DD5,0x7DFA,0x7E1D,0x7E3E,0x7E5F,0x7E7E,\
    0x7E9C,0x7EB9,0x7ED5,0x7EEF,0x7F09,0x7F21,0x7F37,0x7F4D,\
    0x7F61,0x7F74,0x7F86,0x7F97,0x7FA6,0x7FB4,0x7FC1,0x7FCD,\
    0x7FD8,0x7FE1,0x7FE9,0x7FF0,0x7FF5,0x7FF9,0x7FFD,0x7FFE}

const int16_t hSin_Cos_Table[256] = SIN_COS_TABLE;

typedef struct
{
  int16_t hCos;
  int16_t hSin;
} Trig_Components;

/**
  * @brief  This function returns cosine and sine functions of the angle fed in
  *         input
  * @param  hAngle: angle in q1.15 format (-1~0.9999)
  * @retval Trig_Components Cos(angle) and Sin(angle) in Trig_Components format
  */
Trig_Components trig_functions( int16_t hAngle )
{
 int32_t shindex;
 uint16_t uhindex;

 Trig_Components Local_Components;

 /* 10 bit index computation  */
 shindex = ( ( int32_t )32768 + ( int32_t )hAngle );
 uhindex = ( uint16_t )shindex;
 //uhindex /= ( uint16_t )64;
    uhindex = uhindex >> 6;
 /**
  | hAngle    | angle  | std   |
  | (0,16384]   | U0_90  | (0,0.5] |
  | (16384,32767]  | U90_180  | (0.5,0.99]|
  | (-16384,-1]   | U270_360  | (0,-0.5]  |
  | (-16384,-32768] | U180_270  | (-0.5,-1) |
 */
//SIN_MASK        0x0300u
 switch ( ( uint16_t )( uhindex ) & SIN_MASK )
 {          
//0x0200u
   case U0_90:
  Local_Components.hSin = 
            hSin_Cos_Table[( uint8_t )( uhindex )];
  Local_Components.hCos = 
            hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  break;
//0x0300u
   case U90_180:
  Local_Components.hSin = 
            hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  Local_Components.hCos = 
            -hSin_Cos_Table[( uint8_t )( uhindex )];
  break;
//0x0000u
   case U180_270:
  Local_Components.hSin = 
            -hSin_Cos_Table[( uint8_t )( uhindex )];
  Local_Components.hCos = 
            -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  break;
//0x0100u
   case U270_360:
  Local_Components.hSin =  
            -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  Local_Components.hCos =  
            hSin_Cos_Table[( uint8_t )( uhindex )];
  break;
   default:
  break;
 }
 return ( Local_Components );
}

Since the input hAngle is in Q1.15 the format, you can simply draw a picture here; the following is a schematic diagram of the angle hAngle from 0x0000~0xFFFF , as shown below;

Note here that negative numbers are stored in complement form, and the complement of a positive number is equal to itself;

The complement of a negative number is to negate the other bits except the sign bit, and then add 1;

So you can do the math 0xFFFF to express -1 ;

0x8000 express -32768 ;

Because there are unsigned ranges and signed ranges in the Q format, here makes hAngle full use of this 16 bit data, and is compatible with the incoming parameters that can be signed int16 or unsigned uint16 . This is a bit confusing, first look at the picture below;

In the picture above;

• The left side is signed int16 and the right side is unsigned uint16 ;
• The two circles represent the numerical ranges of int16 and respectively uint16 ;
• The waveforms in the green box on the left correspond to the waveforms in the orange box;

There are a few things we should pay attention to here. Whether they are signed or unsigned, their periods are the same;

• Signed integer int16: -32768 ~ 32765,
• Unsigned integer uint16: 0 ~ 65535,

So both use 65536 numbers to represent one period of the sine, which is 2π .

This is a more critical point, so for the key point 0x800 0 , the values ​​represented by signed and unsigned are different;

• Signed integer int16: 0x8000 is represented as -32768;
• Unsigned integer uint16: 0x8000 is represented as 32768;

Therefore, they are exactly one period apart, 65536, so the sine value y represented is the same, as shown by the blue arrows ① and ② in the figure above .

internal implementation

Since the highest bit of the signed integer int16 is the sign bit, here we first convert it into an unsigned integer;

The type used earlier int32 is to prevent data overflow. Adding it here is equivalent to shifting the sine wave by half a cycle, so the mapping relationship between y and x 32768 below needs to be modified according to the actual situation;

/* 10 bit index computation  */
shindex = ( ( int32_t )32768 + ( int32_t )hAngle );
uhindex = ( uint16_t )shindex;
//uhindex /= ( uint16_t )64;
uhindex = uhindex >> 6;

Because the quarter of the sine table that was improved earlier is 256 data, the entire sine cycle should be 1024 subdivided data, which is 10 times of 2, which requires 10 bits ;

• 10 bit The data range is 0~1023 ;
• 16 bit The data range is 0~65535 ;

In order to obtain effective high 10 bit data, the data is shifted to the right 6 bit , as shown below;

Therefore, we can get the following range of data 0 ~ 1023 , 0 ~ 0x400

Therefore, we introduce four masks into the program as identification bits for which quadrant the sine waveform falls in. This also avoids the use of division operations and improves efficiency, as shown below;

#define SIN_MASK        0x0300u
#define U0_90           0x0200u
#define U90_180         0x0300u                  
#define U180_270        0x0000u
#define U270_360        0x0100u

Among them, U0_90 means 0° ~ 90° , and so on;

So why is this mapping relationship?

Shouldn't 0~90° be from 0x000u~0x100u ? Here we briefly explain;

There is such an operation before, the details are as follows;

shindex = ( ( int32_t )32768 + ( int32_t )hAngle );
uhindex = ( uint16_t )shindex;

hAngle The addition here 32768 is equivalent to adding one π , and the sine waveform moves half a cycle to the left ; therefore, the overall mapping relationship must correspond to the original data, as shown below;

Finally, now that we already know which quadrant the waveform is in, we can get a new waveform based on the relationship between the current quadrant and our sine table. There is central symmetry here, symmetry about the y-axis, and we can get the sine and cosine values ​​by simply doing a transformation value;

//SIN_MASK        0x0300u
 switch ( ( uint16_t )( uhindex ) & SIN_MASK )
 {          
//0x0200u
   case U0_90:
  Local_Components.hSin = 
            hSin_Cos_Table[( uint8_t )( uhindex )];
  Local_Components.hCos = 
            hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  break;
//0x0300u
   case U90_180:
  Local_Components.hSin = 
            hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  Local_Components.hCos = 
            -hSin_Cos_Table[( uint8_t )( uhindex )];
  break;
//0x0000u
   case U180_270:
  Local_Components.hSin = 
            -hSin_Cos_Table[( uint8_t )( uhindex )];
  Local_Components.hCos = 
            -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  break;
//0x0100u
   case U270_360:
  Local_Components.hSin =  
            -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
  Local_Components.hCos =  
            hSin_Cos_Table[( uint8_t )( uhindex )];
  break;
   default:
  break;
 }

Summarize

This article briefly introduces the implementation of the sine and cosine functions, and mcsdk makes a simple analysis with reference to ST's algorithm. In it, you need to understand part of the knowledge of Q format for fixed-point operations. You can refer to my " One article to teach you to understand the Q format of C language " This article may contain errors and omissions, please point them out. If you have a better method, please discuss it below.

Latest articlesabout