STM32. About sine
Instead of an introduction
Trigonometric functions, characterized by high consumption of processor time, can negatively affect the choice of budget microcontrollers (without FPU module) for tasks where counting speed is important, for example, attitude control.
Double-precision library trigonometric functions are twice as slow as single-precision trigonometric functions. This is an unfortunate fact.
In the illustration, the operating time of “sin()” and “sinf()”, measured in CPU cycles, on the range [-420° .. +420°] in 5° increments, plus low noise.
The performance of the library “sin()” and “sinf()” is high near the point “0” and decreases by a factor of two in proportion to the modulus of the argument.
There is a need to calculate trigonometry somewhat faster, without FPU, but allowing for worse accuracy:
argument – within the sensitivity of available electronic position sensors, for example, LSM6DSO;
the result is within the accuracy of available stepper motors and servos, for example, HS-5082MG.
Algorithm Requirements
Source Code Requirements:
a) Source code without third party intellectual ballast;
b) Self-contained compiled module based on standard libraries;
c) Cross-platform – no assembler inserts.Requirements for the executable code:
a) Limited amount of binary code – the less, the better;
b) Limited amount of static and dynamic data – the less the better.Interface requirements (usual):
a) float sint( float x );
b) float cost( float x ).Argument requirements:
a) Entity: geometric flat corner;
b) Unit of measure: Angular degree Decimal: 1.0°;
c) Measurement accuracy: +/- 0.05°;
d) Operating range: +/- 5400°;
e) Data type: float, according to IEEE 754.Result requirement:
a) Entity: trigonometric functions “sin” and “cos”;
b) Counting accuracy: +/- 0.005;
c) Data type: float, according to IEEE 754.
Applicable methodology
Comparative analysis of methods for approximating trigonometric functions is beyond the scope of this publication.
The applicable method considers the sine of a flat angle in the operating range from 0 to ¼π.
The sines of angles outside the working range are found through axial symmetry in relation to the result from the working range.
The calculation of the sine of the angle in the working range is replaced by the search for a point on a piece of a one-dimensional spline:
Under certain conditions, spline and sine exhibit similar behavior.
The approximating spline is represented by a Bezier shape. The spline point is calculated geometrically by the de Castelljou reduced method.
The algorithm is optimized to match the sine function.
The cosine is considered as a sine shifted in phase by ¾ π.
Achieved result
Characteristics of the binary code on the Cortex M0 platform:
binary code size – 168 bytes;
stack size – 40 bytes (local variables – 24 bytes);
static and dynamic memory – unnecessarily.
Result with three exact decimal places and a guaranteed error in the fourth decimal place — | 0.0003 |.
The speed of the newly created function on the Cortex M0 platform exceeds the library “sin” and “sinf” by 10 or more times.
The algorithm has options selected via conditional compilation:
slow counting – guaranteed error |0.0003| and less;
fast counting – gives an advantage of ~20 cycles on the Cortex M0, worsening the error to |0.0004|.
The advantage of 15-25 cycles ( ~1us ) is a significant performance boost for low-end controllers at low frequencies.
On the Cortex M4 platform, the fast counting method gives a dubious advantage of 3-5 cycles, commensurate with the measurement error and a deterioration in accuracy to |0.0004|.
In general, testing the algorithm on the Cortex M4 showed an acceptable result with a slight superiority in speed of the “sinf” library function, optimized for the existing FPU.
There is no “sin” function on the chart [ double sin ( double x ) ]showing the worst time in decimal order.
On the Cortex M4 platform, each “sin” call [ double sin ( double x ) ] takes about 2500 cycles.
The final algorithm accepts arguments in degrees of angle and in radians.
The choice of the argument unit is done through conditional compilation.
Instead of a conclusion
The algorithm meets the technical requirements, showing a satisfactory speed and the expected accuracy of the calculation.
There is some potential to improve the performance of the algorithm. From a practical point of view, it is advisable to postpone this potential for the future.
Application Library
File header
/******************************************************************************
* File: cosint.h Created on 31 мар. 2022 г.
*
* Copyright (c) 2022, "Nikolay E. Garbuz" <nik_garbuz@list.ru>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 3 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Authored by "Nikolay E. Garbuz" <nik_garbuz@list.ru>
* Modified by
*
* TAB Size .EQ 4
********************************************************************************/
#ifndef _COSINT_H_
#define _COSINT_H_
#ifdef __cplusplus
extern "C"
{
#endif
//**********************************************************
//---------------- module interface ------------------------
/**
* @brief Fast sine function for approximate calculations
* @param Geometric angle between -5400.00 and 5400.00 degrees
* with two correct numbers after the decimal separator
* @retval Sine with three correct numbers after the decimal separator
*/
float sint( float x );
/**
* @brief Fast cosine function for approximate calculations
* @param Geometric angle between -5400.00 and 5400.00 degrees
* with two correct numbers after the decimal separator
* @retval Cosine with three correct numbers after the decimal separator
*/
float cost( float x );
//**********************************************************
//---------------- module options --__----------------------
/*
* The [COSINT_RAD] directive sets the angle units.
* By default, module functions take a parameter measured in degrees.
*
* Remove or comment out the [#undef COSINT_RAD] line
* to switch the interface [cosint] from degrees to radiant.
*/
#define COSINT_RAD
#undef COSINT_RAD
/*
* The [COSINT_FAST] directive defines the linear proportion method.
* The default is the slow method with a maximum error modulus is |0.0003|.
*
* Remove or comment out the [#undef COSINT_FAST] line to switch module [cosint] to fast mode.
* About ~1.6% speedup, maximum error modulus is 0.0004
*/
#define COSINT_FAST
#undef COSINT_FAST
#ifdef __cplusplus
}
#endif
#endif /* _COSINT_H_ */Source
/******************************************************************************
* File: cosint.c Created on 31 мар. 2022 г.
*
* Copyright (c) 2022, "Nikolay E. Garbuz" <nik_garbuz@list.ru>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 3 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Authored by "Nikolay E. Garbuz" <nik_garbuz@list.ru>
* Modified by
*
* TAB Size .EQ 4
********************************************************************************/
#pragma GCC push_options
#pragma GCC optimize ("O3")
#include "cosint.h"
#include <stdint.h>
//**********************************************************
//---------------- module implementation -------------------
typedef enum
{
coname = -1, siname
} EFName;
inline float cosint( float x, EFName fn );
//---------------- exported module functions ---------------
float sint( float x )
{
return cosint( x, siname );
}
float cost( float x )
{
return cosint( x, coname );
}
//**********************************************************
//------------- helper data and function -------------------
#ifndef M_PI
# define M_PI 3.14159265358979323846 /* pi */
#endif
#ifdef COSINT_RAD
# define R_UNIT M_PI
#else
# define R_UNIT 180.
#endif
#define R_TU ( R_UNIT / 2. ) // The real measure unit in use
#define M_PRC 0x0A // The precision of calculations
#define M_MAX 0x1F // The bit range of calculations
#define M_DAT ( M_MAX >> 1 ) // The bit range of a data
#define M_SGN 0x02 // The sign mask for a result of calculations
#define M_QTR 0x01 // The mask of a circle quarter
#define M_ASQ ( M_SGN | M_QTR ) // The mask of common characteristics for the angular data
#define M_ADV ( 1 << M_DAT ) // The allowed data value
#define M_ADM ( M_ADV - 1 ) // The mask of allowed data value
#define M_MTR (int32_t)( ( M_ADV << M_PRC ) / R_TU )
// Slow and precise, maximum error modulus is |0.0003|
#define DO_SLOW( t, pw, A, B ) ( ( ( t * ( B - A ) + ( ( 1 << pw ) - 1 ) ) >> pw ) + A )
// Fast and rough, about ~1.6% speedup, maximum error modulus is |0.0004|
#define DO_FAST( t, pw, A, B ) ( ( ( t * ( B - A ) ) >> pw ) + A )
#ifdef COSINT_FAST
# define LN_ DO_FAST
#else
# define LN_ DO_SLOW
#endif
// Control data for the spline
#define Ax 0x0000
#define Bx 0x320B
#define Cx 0x655D
#define Dx 0x8000
// helper function
float cosint( float x, EFName fn )
{
int32_t sq, tm;
int32_t A, B, C, D;
// Converting from real to integer
tm = M_MTR;
tm *= x;
tm = tm >> M_PRC;
// Choosing a function via a phase shift by it name
tm += fn & ( M_ADM );
// Extracting a quarter of the angle and a sign of the result
sq = ( tm >> M_DAT ) & M_ASQ;
// Removing periods and reverse from the angle value
// Transforming the angle value to the first quarter
if ( sq & M_QTR )
tm = ( tm & M_ADM ) ^ M_ADM;
else
tm &= M_ADM;
// Calculating a result thru a line proportion
A = Dx;
B = LN_( tm, M_DAT, Cx, Dx );
C = LN_( tm, M_DAT, Bx, Cx );
D = LN_( tm, M_DAT, Ax, Bx );
A = LN_( tm, M_DAT, B, A );
B = LN_( tm, M_DAT, C, B );
C = LN_( tm, M_DAT, D, C );
A = LN_( tm, M_DAT, B, A );
B = LN_( tm, M_DAT, C, B );
A = LN_( tm, M_DAT, B, A );
// Correcting a numeric sign of the result
x = ( sq & M_SGN ) ? -A : A;
// Converting from integer to real and return
x /= M_ADV;
return x;
}
#pragma GCC pop_options
