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:

The decimal logarithm. Required for calculating the magnitude response in decibel:

Sine and cosine: sin(x), cos (x). Required for FFT calculation

Arctangent: arctan(x) or atan(x). Required for the phase response calculating:
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:
Using the case:
Then (4) will corresponds to the form of Maclaurin series:
Let’s choose f(x) as:
Using (6) and (7):
Let’s take only the first 4 terms of the series:
Note: Obviously if the x is smaller then (9) has more accurate result.
Any positive number: a > 0 can be represented as:
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 31. In this case the maximum relative error reaches:
Figure 31: 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:
To find the coefficients, the three equations shall be made. We will use the Chebyshev nodes to reduce the Runge phenomenon (fluctuations are minimized):
In our case:
Using (13), (14), (15), we make up 3 equations for finding the coefficients:
Solving (16):
Note: See the Octave script FunctionsRealization.m, the variable Coef_sqrt2_1
Then:
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 32. In this case the maximum relative error reaches:
Comparing (12) and (19), the error has decreased by 3.6 times!
To reduce the number of multiplications in (18), another entry is used:
Figure 32: 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):
We can use for it the iterative Newton formula:
The formula is easy to prove. Let there is an estimate of the square root:
Suppose the error of this estimate is Δx. Then:
Opening brackets:
If Δx is small, then we can neglect the Δx^2. Then the Δx value:
Then:
Newton Method. Calculation example:
Let’s take the x initial value:
Iteration 1:
Iteration 2:
Iteration 3:
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):
Reasoning in the same way as in paragraph 3.1, we will leave the first 3 terms of the series:
Note: Obviously if the x is smaller then (28) has more accurate result.
Any positive number: a > 0 can be represented as:
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 41. In this case the maximum absolute error reaches:
Figure 41: 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:
from (14). Let’s rewrite (28) in the form:
Using (32) and (14), we make up 3 equations for finding the coefficients:
Solving (33):
Note: See the Octave script FunctionsRealization.m, the variable Coef_ln
Then:
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 42. In this case, the maximum absolute error reaches:
Comparing (31) and (36), the error has decreased by 4 times!
To reduce the number of multiplications in (35), another entry is used:
Figure 42: 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:
5. Sine, Cosine
5.1 Sine
To calculate the sine, we first use the Maclaurin series (6):
Reasoning similarly, as in paragraph 3.1, we will leave the first 4 terms of the series:
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 51. In this case the maximum absolute error reaches:
Figure 51: 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:
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 :
Using (43) and (44), we make up 4 equations for finding the coefficients of a:
Solving (45):
Note: See the Octave script FunctionsRealization.m, the variable Coef_sin
Then:
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 52. In this case the maximum absolute error reaches:
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:
Figure 52: 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:
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π/22π = π… π/2
6. Arctangent
To calculate the arctangent we first use the Maclaurin series (6):
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:
In this case β = 0 … π/4, which corresponds to (51). Then the tan(α) can be calculated:
Rewrite (53) in the form:
Using (52), (54) and (51) we can write for x > 1:
Reasoning similarly for x < 1:
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):
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 61. In this case the maximum absolute error reaches:
Figure 61: 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:
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 :
Using (59) and (60) we make up 4 equations for finding the coefficients:
Solving (61):
Note: See the Octave script FunctionsRealization.m, the variable Coef_atan
Then:
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 62. In this case, the maximum absolute error reaches:
Comparing (58) and (64) the error has decreased by more than 590 times!
Figure 62: 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:
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 81, Figure 82, Figure 83, Figure 84
Figure 81: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_polynomial_float() for x=1…100. Max relative mistake is 0.000163
Figure 82: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_polynomial_int() for x= 65535, 1, 2, … 100. Max relative mistake is 0.001302
Figure 83: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_Newton_float() for x=1…100. Max relative mistake is 0.0000000946
Figure 84: 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 85, Figure 86
Figure 85: Absolute mistake: Δlogarithm of the logarithm10_float() for x= 0.1, 1, 2, …, 99, 100. Max absolute mistake is 0.000365
Figure 86: 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.9537760e3 – 1.4482903e4 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.9537760e3 – 1.4482903e4 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 b15sign
Return
sint16 – sin(input) in format 1Q15: b14*2^1 + b13*2^2 + b12*2^3 + b11*2^4 + … + b0*2^15 and b15sign
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 b15sign
Return
sint16 – cos(input) in format 1Q15: b14*2^1 + b13*2^2 + b12*2^3 + b11*2^4 + … + b0*2^15 and b15sign
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 87, Figure 88
Figure 87: Absolute mistake: Δsine of the sin_float() for x= π … +π. Max absolute mistake is 2.587e04
Figure 88: Absolute mistake: Δsine of the sin_int() for x= π … +π. Max absolute mistake is 2.694e02
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 89, Figure 810, Figure 811
Figure 89: Absolute mistake: Δatan of the atan_float() for x = 100 … 1. Max absolute mistake is 1.036e4
Figure 810: Absolute mistake: Δatan of the atan_float() for x = 1 … 1. Max absolute mistake is 1.032e04
Figure 811: Absolute mistake: Δatan of the atan_float() for x = 1 … 100. Max absolute mistake is 1.036e4