atan_mistake_Chebyshev_4_polynomial
Weimich

Functions Implementation: Square Root, Logarithm, Sine, Cosine, Arctangent. C Code and Octave Script

1. Abbreviation

DSP – Digital Signal Processing
RAM – Random Access Memory
ROM – Read Only Memory
DFT – Discrete Fourier Transform
FFT – Fast Fourier Transform

 

2. Introduction

To implement a DSP algorithm we often need quickly to calculate the values of some functions using the minimal RAM, ROM and time of the microprocessor. In this article I want to describe the C Implementation of these functions using the Octave Tool. What functions are most often required for DSP? These are the functions:
  • Square root. It is required for example to calculate the absolute value of a complex number:
    Functions_formula1
  • The decimal logarithm. Required for calculating the magnitude response in decibel:
    Functions_formula2
  • Sine and cosine: sin(x), cos (x). Required for FFT calculation
  • Arctangent: arctan(x) or atan(x). Required for the phase response calculating:
    Functions_formula3
All these functions are used for calculating and analyzing the Discrete Fourier Transform (DFT). See the “Fast Fourier Transform (FFT) and C Implementation Using the Octave GNU Tool“ and “Goertzel Algorithm and C Implementation Using the Octave GNU Tool“.

3. Square Root

To calculate the square root will be used two different algorithms:
  • Calculation using several terms of the Maclaurin series and the Chebyshev polynomial approximation
  • Newton’s iterative algorithm

3.1. Square Root. Taylor/Maclaurin series. Chebyshev Approximation

Taylor series:

Functions_formula4

Using the case:

Functions_formula5

Then (4) will corresponds to the form of Maclaurin series:

Functions_formula6

Let’s choose f(x) as:

Functions_formula7

Using (6) and (7):

Functions_formula8

Let’s take only the first 4 terms of the series:

Functions_formula9

Note: Obviously if the x is smaller then (9) has more accurate result.
Any positive number: a > 0 can be represented as:
Functions_formula10_11
Using the Octave script FunctionsRealization.m (see the array sqrt_relative_mistake2) is calculated the relative error: Δsqrt(1+x)/sqrt (1+x) by the formula (9) in the range x = -0.3…0.4. See the Figure 3-1. In this case the maximum relative error reaches:

Functions_formula12

sqrt_mistake_Maclaurin_4_polynomial

Figure 3-1: Relative mistake: Δsqrt(1+x)/sqrt(1+x) of the square root calculation with 4 elements of the Maclaurin Series

This accuracy can be increased if the coefficients of the polynomial (9) are recalculated using approximation formulas. In the case it is possible to reduce the error caused by discarding the terms of the infinite Maclaurin series.
Let’s rewrite (9) in the form:

Functions_formula13

To find the coefficients, the three equations shall be made. We will use the Chebyshev nodes to reduce the Runge phenomenon (fluctuations are minimized):

Functions_formula14

In our case:

Functions_formula15

Using (13), (14), (15), we make up 3 equations for finding the coefficients:

Functions_formula16

Solving (16):

Functions_formula17

Note: See the Octave script FunctionsRealization.m, the variable Coef_sqrt2_1
Then:

Functions_formula18

Using Octave script FunctionsRealization.m (see the sqrt_relative_mistake2_1) is calculated the relative error: Δsqrt(1+x)/sqrt (1+x) by the (18) in the range x = -0.3…0.4. See the Figure 3-2. In this case the maximum relative error reaches:
Functions_formula19
Comparing (12) and (19), the error has decreased by 3.6 times!
To reduce the number of multiplications in (18), another entry is used:
Functions_formula20

sqrt_mistake_Chebyshev_4_polynomial

Figure 3-2: Relative mistake: Δsqrt(1+x)/sqrt(1+x) of the square root calculation with 4 elements of the Chebyshev Approximation

Note: If the obtained accuracy (19) is insufficient for your applications, then you need to take more terms of the series.

3.2. Square Root. Newton’s Method.

The x shall be found in (21):

Functions_formula21

We can use for it the iterative Newton formula:

Functions_formula22

The formula is easy to prove. Let there is an estimate of the square root:
Functions_formula22_1
Suppose the error of this estimate is Δx. Then:

Functions_formula23

Opening brackets:

Functions_formula24

If Δx is small, then we can neglect the Δx^2. Then the Δx value:

Functions_formula25

Then:

Functions_formula26

Newton Method. Calculation example:

Functions_formula26_1

Let’s take the x initial value:

Functions_formula26_2

Iteration 1:

Functions_formula26_3

Iteration 2:

Functions_formula26_4

Iteration 3:

Functions_formula26_5

 

Notes:
  • Despite the fact that this iterative formula converges quite quickly, it is impossible to say exactly how many iterations will be required to obtain a result with sufficient accuracy. It depends on the choice of the initial x value. Usually x=a is taken as the initial value
  • This algorithm works well, for example, in the case of a slowly changing parameter measurements. Then you can take the already calculated value of the square root from the previous dimension as the initial value

4. Logarithm

To calculate the logarithm, we use Maclaurin series (6):

Functions_formula27

Reasoning in the same way as in paragraph 3.1, we will leave the first 3 terms of the series:

Functions_formula28

Note: Obviously if the x is smaller then (28) has more accurate result.
Any positive number: a > 0 can be represented as:

Functions_formula29_30

Using the Octave script FunctionsRealization.m (see the ln_absolute_mistake) is calculated the absolute error: Δln (1+x) by (28) in the range x = -0.32…0.36. See the Figure 4-1. In this case the maximum absolute error reaches:

Functions_formula31

ln_mistake_Maclaurin_3_polynomial

Figure 4-1: Absolute mistake: Δln(1+x) of the natural logarithm calculation with 3 elements of the Maclaurin Series

 

As in paragraph 3.1 to increase the accuracy of the calculations we will use the approximation Chebyshev nodes:

Functions_formula31_1

from (14). Let’s rewrite (28) in the form:

Functions_formula32

Using (32) and (14), we make up 3 equations for finding the coefficients:

Functions_formula33

Solving (33):
Functions_formula34
Note: See the Octave script FunctionsRealization.m, the variable Coef_ln
Then:

Functions_formula35

Using the Octave script FunctionsRealization.m (see the ln_absolute_mistake1) is calculated the absolute error: Δln(1+x) by (35) in the range x = -0.32…0.36. See the Figure 4-2. In this case, the maximum absolute error reaches:

Functions_formula36

Comparing (31) and (36), the error has decreased by 4 times!
To reduce the number of multiplications in (35), another entry is used:
Functions_formula37

ln_mistake_Chebyshev_3_polynomial

Figure 4-2: Absolute mistake: Δln(1+x) of the natural logarithm calculation with 3 elements of the Chebyshev Approximation

 

Notes:
  • If the obtained accuracy (36) is insufficient for your applications, then you need to take more terms of the series
  • In the future we will need a decimal logarithm. To move from the natural logarithm is used the basic logarithmic formula:

     

    Functions_formula38_39

5. Sine, Cosine

5.1 Sine

To calculate the sine, we first use the Maclaurin series (6):

Functions_formula40

Reasoning similarly, as in paragraph 3.1, we will leave the first 4 terms of the series:

Functions_formula41

Using the Octave script FunctionsRealization.m (see the sin_absolute_mistake) is calculated the absolute error: Δsin(x) by (41) in the range x = -π…π. See the Figure 5-1. In this case the maximum absolute error reaches:

Functions_formula42

sin_mistake_Maclaurin_4_polynomial

Figure 5-1: Absolute mistake: Δsin(x) of the sine calculation with 4 elements of the Maclaurin Series

 

As in paragraph 3.1 to increase the accuracy of calculations we will use the Chebyshev approximation. Rewrite the (41) in the form:

Functions_formula43

Taking advantage of the odd sine: sin (- x) = – sin(x) is chosen 9 = 4*2+1 Chebyshev nodes and using only the first four nodes to determine the coefficients a1, …, a4 :

Functions_formula44

Using (43) and (44), we make up 4 equations for finding the coefficients of a:

Functions_formula45

Solving (45):

Functions_formula46

Note: See the Octave script FunctionsRealization.m, the variable Coef_sin
Then:

Functions_formula47_1

Using Octave script FunctionsRealization.m (see the sin_absolute_mistake1) is calculated the absolute error: Δsin(x) by (47) in the range x = -π…π. See the Figure 5-2. In this case the maximum absolute error reaches:

Functions_formula48

Comparing (42) and (48) the error has decreased by more than 200 times!
To reduce the number of multiplications in (47), another entry is used:

Functions_formula49

sin_mistake_Chebyshev_4_polynomial

Figure 5-2: Absolute mistake: Δsin(x) of the sine calculation with 4 elements of the Chebyshev Approximation

 

Note:
If the obtained accuracy (36) is insufficient for your applications, then you need to take more terms of the series

 

5.2 Cosine

To calculate the cosine is used the formula:

Functions_formula50

The input x value is -π ≤ x ≤ π, then the range for the sine argument in (50) is -π/2 ≤ x + π/2 ≤ 3π/2. If the argument x + π/2 > π that corresponds to the range π… 3π/2, then it is necessary to subtract 2π from the x: π-2π… 3π/2-2π = -π… -π/2

 

6. Arctangent

To calculate the arctangent we first use the Maclaurin series (6):
Functions_formula51
This series converges only for x: -1 … 1, while the arctan function corresponds to the values α = -π/4 … π/4. Let’s define a series for the values of the argument: x > 1 for the function values: α = π/4 … π/2. The angle α can be written in the form:
Functions_formula52
In this case β = 0 … π/4, which corresponds to (51). Then the tan(α) can be calculated:

Functions_formula53

Rewrite (53) in the form:

Functions_formula54

Using (52), (54) and (51) we can write for x > 1:

Functions_formula55

Reasoning similarly for x < -1:

Functions_formula56

Notes:
  • The series (51), (55) and (56) completely define arctan (x) for any x
  • Formula (55) is also true for x=1, and formula (56) is also true for x = -1
Reasoning in the same way as in the previous paragraphs, we will leave the first 4 terms of the series (51):

Functions_formula57

Using Octave script FunctionsRealization. m ( see the sin_absolute_mistake) is calculated the absolute error: Δarctan(x) by (57) in the range x = -1…1. See the Figure 6-1. In this case the maximum absolute error reaches:

Functions_formula58

atan_mistake_Maclaurin_4_polynomial

Figure 6-1: Absolute mistake: Δatan(x) of the sine calculation with 4 elements of the Maclaurin Series

 

To increase the accuracy of calculations we will use the Chebyshev approximation. Rewrite (57) in the form:

Functions_formula59

Taking advantage of the odd arctan: arctan (- x) = – arctan(x) is chosen 9 = 4*2+1 Chebyshev nodes and using only the first four nodes to determine the coefficients a1, …, a4 :

Functions_formula60

Using (59) and (60) we make up 4 equations for finding the coefficients:

Functions_formula61

Solving (61):

Functions_formula62

Note: See the Octave script FunctionsRealization.m, the variable Coef_atan
Then:

Functions_formula63

Using Octave script FunctionsRealization.m (see the atan_absolute_mistake1) is calculated the absolute error: Δatan(x) by (63) in the range x = -1…1. See the Figure 6-2. In this case, the maximum absolute error reaches:

Functions_formula64

Comparing (58) and (64) the error has decreased by more than 590 times!

atan_mistake_Chebyshev_4_polynomial

Figure 6-2: Absolute mistake: Δatan(x) of the sine calculation with 4 elements of the Chebyshev Approximation

 

To reduce the number of multiplications in (63), another entry is used:

Functions_formula65

Note
If the obtained accuracy (36) is insufficient for your applications, then you need to take more terms of the series

 

7. Octave GNU file FunctionsRealization.m

FunctionsRealization.m Script:
  • Calculates the coefficients of the Chebyshev approximation for the functions square root, logarithm, sine, arctangent
  • Calculates the approximated functions and compares them with the Octave Tool samples and estimate an approximation error
  • Prepares all the graphs used in the article
  • Generates test data arrays for the C code

8. C Language Files

8.1 square_root.c/h

square_root.c/h is supported square root functions

/*
Square Root
Float Arithmetic
Polynomial sqrt(1+x) = 1 + x (0.50016487 + x (-0.12822807 + 0.059227649 x))
0.7 < 1 + x <= 1.4
-0.3 < x <= 0.4
Input
float input – input value
Return
float – sqrt(input) Note: if input number is 0.0 or negative then return value is 0.0
*/
float sgrt_polynomial_float(float input);

/*
Square Root
Fixed Point Arithmetic
Polynomial sqrt(1+x) = 1 + x (0.50016487 + x (-0.12822807 + 0.059227649 x))
0.7 < 1 + x <= 1.4
-0.3 < x <= 0.4
Input
uint16 input – input value: 0…65535
Return
uint16 – sqrt(input) in format: b15*2^7 + b14*2^6 + … + b8*2^0 + b7*2^-1 + b6*2^-2 + … + b0*2^-8
Note: Result = ReturnValue/OUT_SQRT_FACTOR, where OUT_SQRT_FACTOR = 256
*/
uint16 sgrt_polynomial_int(uint16 input);

/*
Square Root
Float Arithmetic
Newton Method x(n+1) = (x(n)+ a/x(n))/2
Input
float input – input value
Return
float – sqrt(input) Note: if input number is 0.0 or negative then return value is 0.0
*/
float sgrt_Newton_float(float input);

/*
Square Root
Fixed Point Arithmetic
Newton Method x(n+1) = (x(n)+ a/x(n))/2
Input
uint16 input – input value: 0…65535
Return
uint16 – sqrt(input) in format: b15*2^7 + b14*2^6 + … + b8*2^0 + b7*2^-1 + b6*2^-2 + … + b0*2^-8
Note: Result = ReturnValue/OUT_SQRT_FACTOR, where OUT_SQRT_FACTOR = 256
*/
uint16 sgrt_Newton_int(uint16 input);

The square root functions were tested in main() using the test arrays: input_array_float[ ] and input_array_int[ ]. The test results are shown in Figure 8-1, Figure 8-2, Figure 8-3, Figure 8-4

sqrt_polyn_float_1_100_relative_error

Figure 8-1: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_polynomial_float() for x=1…100. Max relative mistake is 0.000163

 

sqrt_polyn_int_ffff_1_100_relative_error

Figure 8-2: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_polynomial_int() for x= 65535, 1, 2, … 100. Max relative mistake is 0.001302

 

sqrt_Newton_float_1_100_relative_error

Figure 8-3: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_Newton_float() for x=1…100. Max relative mistake is 0.0000000946

 

sqrt_Newton_int_ffff_1_100_relative_error

Figure 8-4: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_Newton_int() for x= 65535, 1, 2, … 100. Max relative mistake is 0.000913

 

8.2 logarithm.c/h

logarithm.c/h is supported logarithm functions

/*
Square Root
Float Arithmetic
Polynomial log10(1+x) = log10(e) * ln(1+x) = 0.43429449 ( x (1.0004359 + x (-0.52194273 + 0.33578411 x)))
0.68 < 1 + x <= 1.36
-0.32 < x <= 0.36
Input
float input – input value
Return
float – log10(input) Note: if input number is 0.0 or negative then return value is 0.0
*/
float logarithm10_float(float input);

/*
Decimal Logarithm
Fixed Point Arithmetic
Polynomial log10(1+x) = log10(e) * ln(1+x) = 0.43429449 ( x (1.0004359 + x (- 0.52194273 + 0.33578411 x)))
0.68 < 1 + x <= 1.36
-0.32 < x <= 0.36
Input
uint16 input – input value: 1…65535 => log10(1)=0 … log10(65535)=4.8164735
Return
sint16 – log10(input) in format: b14*2^2 + b13*2^1 + b12*2^0 + b11*2^-1 + … + b0*2^-12
Note: Result = ReturnValue/OUT_LOGARITHM_FACTOR, where OUT_LOGARITHM_FACTOR = 4096
*/
sint16 logarithm10_int(uint16 input);

The logarithm10 functions were tested in main() using the test arrays: input_array_float[ ] and input_array_int[ ]. The test results are shown in Figure 8-5, Figure 8-6

log10_polyn_float_0_1_1_100_absolute_error

Figure 8-5: Absolute mistake: Δlogarithm of the logarithm10_float() for x= 0.1, 1, 2, …, 99, 100. Max absolute mistake is 0.000365

 

log10_polyn_int_FFFF_1_100_absolute_error

Figure 8-6: Absolute mistake: Δlogarithm of the logarithm10_int() for x= 65535, 1, 2, …, 99, 100. Max absolute mistake is 0.000539

 

8.3 sine.c/h

Sine.c/h is supported sine and cosine functions

/*
Sine
Float Arithmetic
Polynomial sin(x) = x (0.9992497 + x^2 (-0.16564582 + x^2 (7.9537760e-3 – 1.4482903e-4 x^2)))
-pi <= x <= pi
Input
float input – input value
Return
float – sin(input) Note: if input > pi or input < -pi then return value is 0.0
*/
float sin_float(float input);

/*
Cosine
Float Arithmetic
cos(x) = sin(x + pi/2)
-pi <= x <= pi
Input
float input – input value
Return
float – cos(input) Note: if input > pi or input < -pi then return value is 0.0
*/
float cos_float(float input);

/*
Sine
Fixed Point Arithmetic
Polynomial sin(x) = x (0.9992497 + x^2 (-0.16564582 + x^2 (7.9537760e-3 – 1.4482903e-4 x^2)))
-pi <= x <= pi
Input
sint16 input – input value: -PI_INT_VALUE=-25736 … PI_INT_VALUE=25736
input format 3Q13: b14*2^1 + b13*2^0 + b12*2^-1 + b11*2^-2 + … + b0*2^-13 and b15-sign
Return
sint16 – sin(input) in format 1Q15: b14*2^-1 + b13*2^-2 + b12*2^-3 + b11*2^-4 + … + b0*2^-15 and b15-sign
Note: Result = ReturnValue/OUT_SINE_FACTOR, where OUT_SINE_FACTOR = 32768
*/
sint16 sin_int(sint16 input);

/*
Cosine
Fixed Point Arithmetic
cos(x) = sin(x + pi/2)
-pi <= x <= pi
Input
sint16 input – input value: -PI_INT_VALUE=-25736 … PI_INT_VALUE=25736
input format 3Q13: b14*2^1 + b13*2^0 + b12*2^-1 + b11*2^-2 + … + b0*2^-13 and b15-sign
Return
sint16 – cos(input) in format 1Q15: b14*2^-1 + b13*2^-2 + b12*2^-3 + b11*2^-4 + … + b0*2^-15 and b15-sign
Note: Result = ReturnValue/OUT_SINE_FACTOR, where OUT_SINE_FACTOR = 32768
*/
sint16 cos_int(sint16 input);

The sine / cosine functions were tested in main() using the test arrays: input_array_float[ ] and input_array_int[ ]. The test results are shown in Figure 8-7, Figure 8-8

sin_polyn_float_minus_pi_plus_pi_absolute_error

Figure 8-7: Absolute mistake: Δsine of the sin_float() for x= -π … +π. Max absolute mistake is 2.587e-04

 

sin_polyn_int_minus_pi_plus_pi_absolute_error

Figure 8-8: Absolute mistake: Δsine of the sin_int() for x= -π … +π. Max absolute mistake is 2.694e-02

 

Note: The error of the sin_int() is large (0.02694). This is caused by using 16 bits fixed point arithmetic. To reduce the error, you can use 32 bits fixed point arithmetic instead of 16 bits.

 

8.3 atan.c/h

atan.c/h is supported arctangent function

/*
Arctangent
Float Arithmetic
Polynomial atan(x) = x (0.99904621 + x^2 (-0.3198728 + x^2 (0.14370242 – 0.037537854 x^2)))
-1 <= x <= 1
Polynomial atan(x) = +/- pi/2 – x^-1 (0.99904621 + x^-2 (-0.3198728 + x^-2 (0.14370242 – 0.037537854 x^-2)))
|x| > 1
Input
float input – input value
Return
float – atan(input)
*/
float atan_float(float input);

The atan function was tested in main() using the test arrays: input_array_float 1/2/3[ ]. The test results are shown in Figure 8-9, Figure 8-10, Figure 8-11

atan_polyn_float_minus_1_minus_100_absolute_error

Figure 8-9: Absolute mistake: Δatan of the atan_float() for x = -100 … -1. Max absolute mistake is 1.036e-4

 

atan_polyn_float_minus_1_plus_1_absolute_error

Figure 8-10: Absolute mistake: Δatan of the atan_float() for x = -1 … 1. Max absolute mistake is 1.032e-04

 

atan_polyn_float_plus_1_plus_100_absolute_error

Figure 8-11: Absolute mistake: Δatan of the atan_float() for x = 1 … 100. Max absolute mistake is 1.036e-4

 

9. Download the .m /.c /.h Files

You can download the files:
square_root.c/h
logarithm.c/h
sine.c/h
atan.c/h
FunctionsRealization.m
with the button:
Download
 
 

Views: 1,018