MATLAB Examples

# 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

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)]```

or

```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 ```
```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 ```
``` 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...