Fixed-Point ATAN2 using CORDIC

This demo shows how a fixed-point cordic algorithm to calculate a phase from polar coordinates (arctan) can be implemented in MATLAB.

It also shows how:

C Code generation is used to speed-up the MATLAB function
Embedded MATLAB is used to incorporate this function into Simulink diagrams

CORDIC
Floating-point
Fixed-point
Fixed-point versus Floating point
Accelerate using EMLMEX

CORDIC

CORDIC (digit-by-digit method, Volder's algorithm) (for COordinate Rotation DIgital Computer) is a simple and efficient algorithm to calculate hyperbolic and trigonometric functions. It is commonly used when no hardware multiplier is available (e.g., simple microcontrollers and FPGAs) as the only operations it requires are addition, subtraction, bitshift and table lookup. (Source: http:\\wikipedia.org\wiki\cordic)

CORDIC can be used to calculate a number of different functions. This explanation shows how to use CORDIC in rotation mode to calculate the angle of a vector (x,y), which equals MATLAB's built-in function ATAN2.

help atan2

 ATAN2  Four quadrant inverse tangent.
    ATAN2(Y,X) is the four quadrant arctangent of the real parts of the
    elements of X and Y.  -pi <= ATAN2(Y,X) <= pi.
 
    See also ATAN.

    Overloaded methods:
       distributed/atan2

    doc atan2

We will first assume that the angle of the vector is between 0 and 45 degrees or [0 pi/4]. The inverse algorithm (which works for arctan, arccos and arcsin) will try to rotate the vector towards (1,0). So if the y coordinate of the vector (or the imaginary part) is positive CORDIC will rotate with a negative angle, otherwise it will rotate with a positve angle. By keeping track of how many degrees we have rotated the original vector we know the angle of that vector. That's all!

The power of the CORDIC algorithm is that we pick the angles for rotation in such a way that the corresponding rotation matrix only contains powers of 2, which in fixed-point can be easily implemented using bitshifts (no multiplications needed!)

The rotation matrix R for an angle phi equals:

R = [cos(phi) -sin(phi);
     sin(phi) cos(phi)]

R = K*[1        -tan(phi);
       tan(phi) 1         ]

with

K = 1/sqrt(1 + tan(phi)^2)

Because we are only interested in the angle of the vector and not in the magnitude, we can ignore the factor K for now. We have to deal with it later when we determine the appropriate fixed-point data types.

Now we let tan(phi) equal powers of 2:

phi = 2.^-(1:5)

phi =

  Columns 1 through 3

    0.5000    0.2500    0.1250

  Columns 4 through 5

    0.0625    0.0313

This corresponds to the following angles:

atan(phi)

ans =

  Columns 1 through 3

    0.4636    0.2450    0.1244

  Columns 4 through 5

    0.0624    0.0312

or in degrees:

atan(phi)*180/pi

ans =

  Columns 1 through 3

   26.5651   14.0362    7.1250

  Columns 4 through 5

    3.5763    1.7899

It can be proven that you can rotate to any angle between 0 and 45 degrees by rotating clock-wise or counter clock-wise with these angles. The number of iterations (in this case only 5) determines the precicion of the algorithm.

Floating-point

In floating-point the algorithm is astonishing simple. Let's try with an angle of 0.25 (rad).

v = [cos(0.25), sin(0.25)]'; % Starting vector (input)
a = 0;
for n = 1:5;
        a = a + sign(v(2))*phi(n); % Keep track of rotation
        R = [1,                 sign(v(2))*2^(-n);
            -sign(v(2))*2^(-n), 1                   ]; % Rotation matrix
        v = R*v; % Rotate
end

The result after 5 iterations equals:

a =

    0.2813

Fixed-point

The input should accomodate the unit circle, making necessary a datatype of 16 bits with a 14 bit fraction length. The range of this fixed-point data type is [-2 2]. This is usefull, because the gain of the CORDIC is max 1.1644. So overflows are not possible in the rotation of the input vector.

Tcar = numerictype(1,16,14)
Fcar = fimath('SumWordLength', 16, 'SumFractionLength', 14, 'SumMode', 'Specify');

 
Tcar =
 

          DataTypeMode: Fixed-point: binary point scaling
                Signed: true
            WordLength: 16
        FractionLength: 14

Because we want to work with bitshifts we have to abandon the matrix notation and work with the seperate real and imaginary parts of the complex vector (or x and y from a the vector (x,y), whatever you prefer).

im = fi(sin(0.25), 'NumericType', Tcar,'FiMath', Fcar);
re = fi(cos(0.25), 'NumericType', Tcar,'FiMath', Fcar);

With regards to the output: the temporary result of the cordic algorithm is always between +pi/4 and -pi/4 which is within the range of 16 bit fixed-point number with a 15 bit fraction length. However, the output of the algorithm is always between -pi and pi, so a fraction length of 13 bits makes more sense, as more precision is lost in the output.

Trad = numerictype(1,16,13);
Frad = fimath('SumWordLength', 16, 'SumFractionLength', 13, 'SumMode', 'Specify');

Now with regards to the number of rotations/iterations: the trick is to find trade-off between fixed-point precision and algorithmic precision. The algorithmic precision is approximately atan(2.^-N). The fixed-point precision is 2^-13. The trade-off is around N = 12.

N = 12;
angles = fi(atan(2.^-(1:N)),'NumericType', Trad, 'FiMath', Frad);

The fixed-point CORDIC including bitshifts can be written as:

reg = fi(0, 'NumericType',Tcar, 'FiMath', Fcar); % Temporary register

a = fi(0,'NumericType', Trad, 'FiMath', Frad); % Starting angle (0 radians)

for i = 1:N,

    if bitget(im,16), % negative angle
        a = a - angles(i);

        reg = re - bitshift(im,-i);
        im = bitshift(re,-i) + im;
        re = reg;

    else
        a = a + angles(i);

        reg = re + bitshift(im,-i);
        im = -bitshift(re,-i) + im;
        re = reg;

    end

end

The fixed-point result after 5 iterations equals:

 
a =
 
    0.2501

          DataTypeMode: Fixed-point: binary point scaling
                Signed: true
            WordLength: 16
        FractionLength: 13

             RoundMode: nearest
          OverflowMode: saturate
           ProductMode: FullPrecision
  MaxProductWordLength: 128
               SumMode: SpecifyPrecision
         SumWordLength: 16
     SumFractionLength: 13
         CastBeforeSum: true

For angles between pi/4 and pi/2 we use the following trigonometric equality:

tan(pi/2 - a) == 1/tan(a);

This means that we can simply swap the imaginary and real part of the input vector and then substract the result from pi/2.

For starting vectors in one of the other three quadrants (or angles other than [0 pi/2] we always start in the first quadrant by taking the absolute values of both the imaginary and real part. We then have to substract or add the result of the CORDIC algorithm from/to pi/2 or 0 radians. See the implementation in atan2_fixpt.m.

Fixed-point versus Floating point

Let's test how accurate the fixed-point implementatoin is in comparison with floating-point. Because we are interested in the accuracy of the algorithm we look at how well the CORDIC perform in calculating the angle of a fixed-point input vector:

phis = 2*pi*(rand(1,1e3) -0.5); % Some random angles
phis_fixpt = fi(phis,'NumericType', Trad, 'FiMath', Frad); % Rounded towards fixed-point representation
phis_cordic = phis_fixpt; % Pre-allocate fixed-point output vector.

Floating point:

phis_atan2 = atan2(sin(double(phis_fixpt)),cos(double(phis_fixpt)));

Fixed-point

sinphis = fi(sin(phis), 'NumericType', Tcar, 'FiMath', Fcar);
cosphis = fi(cos(phis), 'NumericType', Tcar, 'FiMath', Fcar);

phis_cordic = atan2_fixpt(sinphis, cosphis);

Now look at the root mean squared error of the algorithm and the histogram:

err_cordic = sqrt(mean((double(phis_cordic)-phis_atan2).^2))

err_cordic =

  1.7811e-004

Accelerate using EMLMEX

The fixed-point toolbox provides the possibility of accelerating fixed-point m-code using EMLMEX, which generates C-MEX code from m-code:

emlmex -eg {sinphis,cosphis} -o atan2_fixpt_mex atan2_fixpt

Now let's compare the speed of the original m-file

tic
   phis_cordic = atan2_fixpt(sinphis, cosphis);
toc

Elapsed time is 8.839059 seconds.

with the MEX file:

tic
    phis_cordic = atan2_fixpt_mex(sinphis, cosphis);
toc

Elapsed time is 0.014633 seconds.

The speed-up is more then 100 times...

Fixed-Point ATAN2 using CORDIC

Contents

CORDIC

Floating-point

Fixed-point

Fixed-point versus Floating point

Accelerate using EMLMEX