Compare the error between CORDIC and Polynomial approximation for trigonometic functions

Intro

Computing sin(), and cos() function are heavily needed in signal processing systems, FFT, DCT. Usually, we relly on the math.h package in C for the floating-point computation.

However, the full-precesion calculation for such functions (in math.h) could be very time consuming. Two methods are compared here, mainly for the complexity vs. error.

CORDIC and Polynomical

  • CORDIC is believed to be superior when there’s no hardware multipliers available, because it only need international shift and add.
  • Polynomial Approximation uses the Taylor series to calculate the sin() and cos().

A quick note on CORDIC can be found in this Wikipeadia page on CORDIC. It is an iterative approach that each iteration reduces the error.

Using Taylor series (that expand at 0), the sin() and cos() can be represented as:
$$cos(x)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+ …$$
$$
sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+…$$
The result becomes more accurate along with the increase of orders. For the example, cos(x) use 6 orders, while sin(x) uses 7 orders.

For sin() and cos(), the error can be quantified as by the number of binary bits (nb), which represent the dynamical range:

$$nb = log_2\frac{range}{resolution}=log_2\frac{max-min}{2*error}=log_2\frac{1-(-1)}{2*error}$$

  • The “CORRECT” number is given by the matlab function “sin()”, (possiblility using the .math.h C function) which might not be so accurate in iteself. But it is believed to be accurate enough to calculate the error of our scheme (since it’s much more accurate than these 2 approaches mentioned before).

  • When calculating the value of sin(), we are using the format of (64-bit) double floating point in matlab. This might be very different if you are implemented that in a FPGA or ASIC, when you might be restricted to use a fixed-point ALU. Since normally the calculation precision is less than in this test, your error will be larger in reality.

  • the “order” in the polynomial approximation means the maximum number of order in the Taylor series. A n-order polynomial results only adds up n/2 items($\frac{x^n}{n!} $).

result

The dynamical ranges of those two approaches are listed below, where the input ranges form $-\pi/2$ to $\pi/2$.

Dynamical range of CORDIC

Dynamical range of Taylor approximation

For CORDIC, the error is more uniform across inputs range. However, for Taylor series, the more distant form the expanding root (0), the larger the error. To reduce the error around $\pi/2$, larger inputs are flipped.

Rule of thumb results is that: 30 CORDIC iterations results to 30-bit resolutions, which is about 10 taylor iterations.

The complexity for 30-bits is compared below. Each CORDIC iteration need 2 additions and a comparison with 0. The terms for Taylor series ($\frac{1}{n!}$) are pre-stored, so for each order, we need 2 multiplications, one for $x^n$ (generated form $x^(n-2)$), one for $\frac{1}{n!}$.

add. compare mul.
CORDIC 60 (30x2) 30 0
Taylor polynominal 5 (10/2) 0 10 (10/2x2)

Sum

The error plot can be used to understand the usage of CORDIC and Taylor approximation.
So generally:

  • CORDIC is prefered if no FPU is presented.
  • Taylor can be faster if FPU is presented, since for FPU, the latency for mul. and add. are more or less equal.
  • For low-power, CORDIC is prefered, since add, especially fixed-point add is much more cheaper than FPU. But it takes longer latencies, and if the control power dominates (no execution power), then Taylor should be used.

Code

Compare error script:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
clear;close all;%be careful with it. You might lose your work state.
iterations=5:5:35; %number of iterations for CORDIC
orders=2:2:14; %number of orders for taylor polynimial approximation
x=(-pi/2:pi/4/200:pi/2); % x axis. x is the 'x' (input) in 'cos(x)';
sinx = sin(x); % exact result
colors = {'-b', '-g', '-r', '-c', '-m', '-y', 'k', 'w'};
%% plot cordic error
figure(1);
legendCordic = {}; %plot legend
for index = 1:length(iterations),
iteration = iterations(index);
% Cordic computation
sinxCordic = zeros(size(x));
for num = 1:length(x),
tempOut = cordic(x(num), iteration); % compute cos, sin
sinxCordic(num) = tempOut(2);
end;
% Assume that the LSB denote the twice as the error.
% Thus, Number_of_bit = log2( (max-value - min-value) / (LSB distance) ).
% This is the dynamical range.
nbCordic = log2 ( (1-(-1)) ./ (2.* abs(sinxCordic-sinx)) ) ;
legendCordic{index} = sprintf('%d iterations', iteration); % update legend
plot(x, nbCordic, colors{index});hold on;
end
legend(legendCordic);
xlabel('input values (-pi/2 to pi/2)');
ylabel('dynamical range of error (number of bits)');
title('For sin(x), dynamical range increases with iterations of CORDIC');
axis([-1.6 1.6 0 60]);
set(gca,'fontsize',16);
%% plot taylor series, polynominal error
figure(2);
legendTaylor = {}; %plot legend
for index = 1:length(orders),
order = orders(index);
% taylor computation
sinxTaylor = zeros(size(x));
for num = 1:length(x),
tempOut = taylor(x(num), order); % compute cos, sin
sinxTaylor(num) = tempOut(2);
end;
% Assume that the LSB denote the twice as the error.
% Thus, Number_of_bit = log2( (max-value - min-value) / (LSB distance) ).
% This is the dynamical range.
nbCordic = log2( (1-(-1)) ./ (2.* abs(sinxTaylor-sinx)) );
legendTaylor{index} = sprintf('%d orders', order); % update legend
plot(x, nbCordic, colors{index});hold on;
end
legend(legendTaylor);
xlabel('input values (-pi/2 to pi/2)');
ylabel('dynamical range of error (number of bits)');
title('For sin(x), dynamical range increases with taylor series orders');
set(gca,'fontsize',16);
axis([-1.6 1.6 0 60]);

CORDIC function (reference: jacoby, or wikipedia )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
function v = cordic(beta,n)
% This function computes v = [cos(beta), sin(beta)] (beta in radians)
% using n iterations. Increasing n will increase the precision.
beta = mod(beta,2*pi);
sign = 1; %sign of change;
while beta < -pi/2 || beta > pi/2
if beta < 0
beta = beta + pi;
else
beta = beta - pi;
end
sign = - sign; % flip the sign for second or third quadrant
end
% Initialization of tables of constants used by CORDIC
% need a table of arctangents of negative powers of two, in radians:
% angles = atan(2.^-(0:27));
angles = [ ...
0.78539816339745 0.46364760900081 0.24497866312686 0.12435499454676 ...
0.06241880999596 0.03123983343027 0.01562372862048 0.00781234106010 ...
0.00390623013197 0.00195312251648 0.00097656218956 0.00048828121119 ...
0.00024414062015 0.00012207031189 0.00006103515617 0.00003051757812 ...
0.00001525878906 0.00000762939453 0.00000381469727 0.00000190734863 ...
0.00000095367432 0.00000047683716 0.00000023841858 0.00000011920929 ...
0.00000005960464 0.00000002980232 0.00000001490116 0.00000000745058 ];
% and a table of products of reciprocal lengths of vectors [1, 2^-2j]:
Kvalues = [ ...
0.70710678118655 0.63245553203368 0.61357199107790 0.60883391251775 ...
0.60764825625617 0.60735177014130 0.60727764409353 0.60725911229889 ...
0.60725447933256 0.60725332108988 0.60725303152913 0.60725295913894 ...
0.60725294104140 0.60725293651701 0.60725293538591 0.60725293510314 ...
0.60725293503245 0.60725293501477 0.60725293501035 0.60725293500925 ...
0.60725293500897 0.60725293500890 0.60725293500889 0.60725293500888 ];
Kn = Kvalues(min(n, length(Kvalues)));
% Initialize loop variables:
v = [1;0]; % start with 2-vector cosine and sine of zero
poweroftwo = 1;
angle = angles(1);
% Iterations
for j = 0:n-1;
if beta < 0
sigma = -1;
else
sigma = 1;
end
factor = sigma * poweroftwo;
R = [1, -factor; factor, 1];
v = R * v; % 2-by-2 matrix multiply
beta = beta - sigma * angle; % update the remaining angle
poweroftwo = poweroftwo / 2;
% update the angle from table, or eventually by just dividing by two
if j+2 > length(angles)
angle = angle / 2;
else
angle = angles(j+2);
end
end
% Adjust length of output vector to be [cos(beta), sin(beta)]:
v = v * Kn * sign;
return

taylor approximation function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
function v = taylor(beta,n)
% This function computes v = [cos(beta), sin(beta)] (beta in radians)
% using n order approximation. Increasing n will increase the precision.
sign = 1; %sign of change;
beta = mod(beta,2*pi);
% convert to form second or third quadrant
while beta < -pi/2 || beta > pi/2
if beta < 0
beta = beta + pi;
else
beta = beta - pi;
end
sign = - sign;
end
% convert from fourth quadrant
if beta < 0,
beta = -beta;
flip_1 = true;
else
flip_1 = false;
end;
% flip from the second half of first quadrant
% if value is more closed to pi/2 than 0, then approximate on pi/2, reduce error;
if beta > pi/4 - pi/40 , % why it should flip on here, not pi/4? to minimize error
flip_2 = true;
beta = pi/2 - beta;
else
flip_2 = false;
end
sign_sin = 1;
sign_cos = -1;
v = [1, 0];% 0th order
base = 1;% 1/n!
beta_n = 1;% beta^n
for x = 1:n,
base = base/x; % 1/n!
beta_n = beta_n * beta;% beta^n
if mod(x,2)==0;% cos()
v(1) = v(1) + sign_cos * beta_n * base;
sign_cos = - sign_cos;
else % sin ()
v(2) = v(2) + sign_sin * beta_n * base;
sign_sin = - sign_sin;
end
end
% flip back, to the whole first quadrant
if flip_2,
temp = v(1); v(1) = v(2); v(2) = temp;
end
% flip back, to first quadrant and fourth quadrant
if flip_1,
v(2) = -v(2);
end
% convert to quadrant 1,2,3,and 4
v = sign*v;
return