CORDIC-atn.bpe

Given x and y, computes the angle
and the length of the vector from the origin to point (x,y). (revised 12/05)

CORDIC-sincos.bpe

Given an angle in degrees, 0-360, calculates the sine and cosine to x.xxx.

CORDIC is an acronym for **COordinate Rotation Digital
Computer**. The notions behind this computing machinery were
motivated by the need to calculate the trig functions and
inverse trig functions in real time navigation systems. The
CORDIC algorithms require only shifts, adds and table lookups, simple
integer math. So, it is possible to implement a specialized
CORDIC machine, small and fast enough for real time calculations,
dedicated to that one purpose. The algorithm can be coded in
firmware on even small microcontrollers. With slight
modifications in initial conditions and data tables,
the core algorithm can multiply, divide, modulate,
and calculate square roots, hyperbolic functions, exponentials and
logs. It
is this versatility and simplicity, that make CORDIC the preferred
implementation of math functions
on small hand
calculators.

The purpose of this note is to show how to code the algorithms BASIC Stamp. First off, it is an educational exercise, to see how these fascinating algorithms work. Beyond that, sometimes you really need one of these math functions. The alternative is often a table lookup with interpolation, an external coprocessor, a different algorithm (which may turn out to be quite beyond the capability of the Stamp), or a rethinking of what is really required. The BASIC Stamp of course has no floating point unit, although you can add one as an external chip. The same can be said for most small embedded processors, like the PIC or the SX series, or the MSP430 series. The CORDIC method is attractive because it is an integer algorithm close to the hardware. Of course the Stamp will not be fast in its implementation, but in some applications it is fast enough.

You will find very little about CORDIC in math books or other
literature. For example, Jack Crenshaw does not mention it in his
excellent book, *Math Toolkit for
Real Time Programming*, nor is there a single mention of it in
Donald Knuth's seminal tome, *The
Art of Computer Programming*.
I don't know why. Maybe the algorithm is too unsophisticated
or so specialized as to be relegated to hand calculators,
coprocessors and navigational engines. The computational methods
of choice, explained in detail in college math courses and texts,
are based on power series expansions of the functions, or solvers
such as Newton's method. The
series expansion for each mathematical function can present
unique difficulties in terms of convergence, so each function ends
up with a
different piece of computer code, different from the next function in
the library.
They run best on a full general purpose
computer with floating point math capabilities built into the hardware,
and under such conditions they can be optimized
computational speed and accuracy. But for generality and for compatibility with simple integer hardware, the CORDIC
method has no peer.

Here are some links to URLs relating to the CORDIC method (also try Google, I know some of these are out of date):

- Ray Andraka, Consulting
- gurus' guru. By far the most incisive discussion of the subject I know of. Download the pdf file, "Survey of CORDIC Algorithms for FPGAs".
- Grant R. Griffin
- dspGURU. An informative FAQ, with additional links to the online and printed literature.
- Ingo Cyliax
- CORDIC swiss army knife for computing math functions.. An article from Circuit Cellar INK. It even includes a program for calculating sine and cosine on the BASIC Stamp II. A different approach from the one taken here.
- Pitts Jarvis
- "IMPLEMENTING CORDIC ALGORITHMS A single compact routine for computing transcendental functions" An tutorial from Dr. Dobbs Journal, October 1990. code
- Helmut Knaust
- "How Do Calculators Calculate? A tutorial from University of Texas.
- Volder, J.E., 1959, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, V. EC-8, No. 3, pp. 330-334
- seminal paper
- Vladimir Baykov
- Interesting newsgroup posting by Vladimir BAYKOV
- Norbert Lindlbauer, Berkeley Center for New Music
- Discussion of CORDIC architectures with a view toward music synthesizer

Givens transform, the underpinning of CORDIC |
---|

The CORDIC method depends on the Givens rotation transform, which expresses the coordinates of a point P in two different coordinate systems, rotated by a certain angle about their common origin. In the first system point P has coordinates (x,y) and in the second system the same point P has coordinates (x'y'). The second system is rotated by an angle H with relation to the first. (H=phi in the diagram to the right-sorry about the notation). The relation between the coordinates is:

x' = x cos H - y sin H

y' = y cos H + x sin H

Look at a derivation of the Givens transform which illustrates the above formula with a geometrical proof. It also shows the close relation between coordinate rotation and vector rotation.

CORDIC is a successive approximation algorithm. A sequence of successively smaller rotations based on binary decisions hone in on the value we want to find. For example, we start with initial values of x and y, and end up with what we want, the angle whose tangent is y/x and also with the length of the vector SQR(x^2,+y^2).Engineers are familiar with successive approximation, where, say, an unknown voltage is presented at one input, and is first compared with 1/2 of a reference voltage. If the unknown is greater than 1/2 the reference, then the most significant bit of the answer is 1 and the next value to compare with is 3/4 of the reference, otherwise, the most significant bit of the answer is 0 and the next value to compare with is 1/4 of the reference. The sequence of decisions up=1 and down=0 at each step generates the digital binary output for as many steps as the system precision will allow.

In CORDIC, the comparison takes place between the one of the variables and a target value for that variable.

For example, suppose you start with a given angle, and what you want to find is the sine and cosine of that angle. Assume there is an initial unit vector on the x axis. The first decision is whether the given angle is greater than or less than zero. The first rotation is accordingly then either (tan H= +1) or (tan H = -1), that is 45 degrees positive or negative. The angle accumulator now contains either +45 degrees or -45 degrees. The original vector is rotated to that angle with new coordinates (x'y'). The next rotation is (tan H = +/- 1/2) in the direction that rotates toward the given angle. All the computer has to decide is the direction, based on comparing the value in the angle accumulator with the starting angle. The angle accumulator at each step is updated by lookup of the corresponding atn(H) in the table, added to or subtracted from the current value in the angle accumulator. Angles add up. The net rotation is the simple sum of all the small rotations. At each step the coordinates x',y' are updated by means of the Givens transform. Through a sequence of smaller and smaller steps the algorithm ends up with original unit vector rotated so that it coincides with the given starting angle, and the final coordinates x' and y' equal to the sine and cosine of that angle.

Now, suppose you start with a value of y/x, and you want to find the angle that has y/x as its TANgent. That is, you want to calculate the arctangent. The coordinates are first rotated by an angle of 45 degrees (tan H = +/- 1) the direction being counterclockwise if y>0 and clockwise if y<0. That determines the most significant bit. Then the next step is a smaller rotation, (tan H = +/- 1/2), also based on the sign of y. The action taken in sequence is to rotate the coordinates in progressively smaller steps until the transformed x axis lines up with the point P, that is, the y coordinate finally is reduced to zero. At each step the computer has to make only the binary decision of clockwise or counter clockwise depending on the current value of the y component. The binary number generated by the sequence of up/down decisions is the unique angle subtended by (x,y) in the original coordinate system. While the y value is reduced to zero in the final coordinate system, the x value is, in the end, and after adjustment by a constant scale factor, equal to the length of the original vector, SQR(x^2+y^2).

In the calculation of the arctangent and length of the vector, it is the coordinate y that is minimized, while in the calculation of the sine and cosine, it is the difference between the given angle and the rotated angle that is minimized. Both are CORDIC algorithms, one called "rotation" and the other called "vectoring", and there are many more twists of the same kind for calculation of other kinds of quantities by means of successive rotations. The difference has to do with the starting condition and which variable is minimized. But all involve a successive approximation and a binary decision at each step. The variables are linked through the geometry. Very clever.

Recall that tan H = sin H/ cos H, so the above formula can be rewritten as,

x' = cos H (x - y tan H)

y' = cos H (y + x tan H)

One neat trick of CORDIC is to allow only angular rotations of +/- tan H = 1, 1/2, 1/4, 1/8 and so on out to the number of bits to be acquired. In that way, the math computation requires only division by powers of two, shifts, and addition or subtraction. The value of each bit in angular base is arctan 1/2, arctan 1/4 and so on. For example, the translation from the binary value theta = 1011 into standard angular measure is simply theta = arctan(1/2) + -arctan(1/4)+ arctan (1/8) + arctan(1/16). This computation is usually done at the step by step as the x' and y' are calculated. The angles that correspond to ATN 1, ATN 1/2, ATN 1/4,... are stored in a lookup table in the computer memory, and the angle is accumulated at each step by adding or subtracting the successively smaller values. (Note, the Andraka article points out situations where the funny binary base is itself sufficient, without translation to standard measure.)

Here it is in terms that relate x and y at step i+1, to x and y at step i:

The iteration starts with i=0. The decision at each step is the factor d(i), which will take on a value of either +1 or - 1, depending on the sign of the decision variable. The particular variable or condition that determines the d(i), that is the main thing that distinguishes one CORDIC algorithm from another. The core calculation is exactly what is written there. Simply by changing the initial conditions and decision criteria, the core algorithm can be made to perform very different tasks. That is what gives it its versatility. Same machinery, many functions.x(i+1) = cos(atn(1/2^i)) * {x(i) - y(i) * d(i) * (1/2^i)}

y(i+1) = cos(atn(1/2^i)) * {y(i) + x(i) * d(i) * (1/2^i)}

theta = theta + d(i) * arctan (1/2^i)

A second trick, or insight, of CORDIC is that the value of cos H does not depend on the direction of rotation. That is a property of cosine, that cos H = cos -H. That fact is very convenient, because the cosine terms are a product that does not depend on d(i),

The value of cos(atn(H)) approaches unity as H approaches zero, and the limit of the product of terms is a constant number, 0.607252935... It is that many significant digits even after 16 steps. The exact value of the product is easily computed in advance and stored as a constant. This factor is called the CORDIC gain. Some computations (like the length of the vector, or the values of cosine and sine, have to be compensated at some point in the calculation by this gain. That can be done either in the initial conditions or at the end result. In Stamp math, multiplying by 0.607252935 is done as **39797, because 39797/65536 ~= 0.607252935. {Side math note, from simple geometry, cos(atn(1/2^i)) = 1/sqrt(1 + 2^-2i)}cos(atn(1)) * cos(atn(1/2)) * cos(atn(1/4)) * cos(atn(1/8)) * cos(atn(1/16)) * ...

= 0.707107 * 0.894427 * 0.970143 * 0.992278 * 0.998052 * 0.999512 * 0.999878 * ...

scaling for the microcontroller, BASIC Stamp |
---|

The BASIC Stamp and microcontrollers in general work with integers. If we want to represent a fraction, we have to assign some integer to be our unity, and then smaller numbers represent fractional parts of that unity. For example, if number 1 in this system has the value 10000, then the smaller numbers represent fractions from 0/10000, 1/10000, ... , 9999/10000. In some of the other math notes on this site, I have had reason to choose other values for the unity. For example, with a unity of 65536, the fractions are represented in steps of 1/65536, and that is quite a good approximation and is well suited to fractional multiplication using the ** operator.

To use CORDIC on the BASIC Stamp, a scale has to be chosen for both angles and coordinates.

The choice of units is quite arbitrary, and it can be made to work in any units you want. However, the CORDIC algorithm works only with angles in quadrants 1 and 4 (no matter what units are chosen). Angles outside of that range have to be handled with a wrapper that rotates or reflects the problem into that range and compensates at the end, using standard trigonometric identities.

The BASIC Stamp has built in operators for SIN, COS, ATN and HYP. Those functions represent the circle in units of brads, with 256 brads equal to 360 degrees or 2 pi radians. The angular resolution is thus a little less than one degree. Similarly, SIN and COS take on values from -128 (in twos complement) to +128. Thus those values scale the standard range of SIN and COS from -1 to +1, and the full range is thus capable of a resolution of about 1/2%.

However, CORDIC can do better with 16 bits to work with. In the programs below, I have chosen to use the number 10000 for the unity for the coordinates, x and y, and the number 16384 to represent the angle 45 degrees. The angles that can be represented on the Stamp thus range from -32768 (-90 degrees) to +32767 (+90 degrees). The example programs show how to provide the wrapper to translate to standard units and how to cover all four quadrants, or a continuous input variable.

The programs depend on a table of data that translates from the CORDIC binary representation of angle to our standard representation, where 45 degrees is 16384. The words in the table are calculated in advance on a calculator or in Excel as

16384 * arctan (1), 16384 * arctan (1/2), ... , 16384 * arctan (1/32768)

In this calculating these value for the table, it will be necessary to round off. This step leads to a loss of precision in the final answer, and for that reason it pays to use as large a number as possible to represent the range of the 1st and 4th quadrant. For example, instead of 16384 to represent 45 degrees, you could use the number 4500 to have the answer come out directly in degrees * 100, or 7854 to have the answer come out in radians * 10000, or 25000 to have the answer come out as a fraction of pi * 100000. But the round off would be more significant when using these smaller numbers. It is better to start with the highest number possible within the 16 bit framework, to optimize the precision, and then rescale at the end. And of course, going to double precision would be best of all, and the principles described here are the same. But the BASIC Stamp is slow and clunky with those calculations, so I'll start simple and leave the double precision for later.

Granted, it is issues of scaling that will take the most consideration in most applications. The CORDIC machinery will seem trivial in comparison, sitting at the core of wrappers upon wrappers.

Specifics of the BASIC Stamp implementation |
---|

The demo programs are found at the links at the top of this article, there is one program for calculating the angle and vector length, given a point (x,y), and another program for calculating the sine and cosine of a given angle. The programs make use of some program tricks unique to the BASIC Stamp.

CORDIC makes frequent use of the sign of a number. The decisions on whether to add or subtract are made based on the sign of one particular variable. Negative numbers in the BASIC Stamp are represented in twos complement, where bit 15 is the sign. In these programs each variable defined is accompanied by an alias for its sign bit, as

z VAR Word

zs VAR z.Bit15

When it is necessary to calculate a number conditional on its own sign, or on the sign of another number, you will see this kind of sequence.

y = y + ( -zs ^ x + zs )

The goal of this calculation is to either add or subtract the value of a variable x from another variable y, and the add or subtract depends on the sign of a third variable z. The formula in () follows directly from the definition of twos complement. The negative of a number is its ones complement plus 1. If z is negative, then zs in the above equals 1, and -zs= -1, which in twos complement is all bits equal to 1. The exclusive OR of that with the number x forms the ones complement of x and adding zs=1 completes the twos complement negation. On the other hand, if z is positive, then zs=0, then -zs also equals zero, and the XOR with x and adding zero leaves plain old positive x. Overall, the variable x takes on the sign of z. Often the program uses this in the context of a calculation. For example, division of x by a power of two, 2^I, can be done naively like this.

y = y + (x >> i) ' divide x by 2^i

Nevertheless that does not work if x is negative, as it often will be during the CORDIC procedure. Instead, this is done this way:

y = y + (-xs ^ (ABS x >> i) + xs)

That is, the division by 2^i, shift right i, is applied to the absolute value, and the sign is restored. It would be possible to do this with IF-THEN statements, but there are so many of these kind of statements that I find that the IF-THEN logic quickly becomes unwieldy. It is also possible to do these divisions by two simply by shifting the value right, and extending the sign. (This works for powers of two, not for division in general) The BASIC Stamp does not provide an easy way to do this, but I think that would be the way to go in machine language implementations. In PBASIC, this comes out as

y = y + (x>>i | (-xs << (16-i))

That is, the value of x is shifted right i bits, to divide by 2^i, and then the bits on the left are filled in with the original sign, sign extended.

The attached programs show several ways to do the core rotation. One based on IF-THEN decision, others based on applying the sign logic. The different versions use the conditional compilation directive of the latest Stamp IDE, and the selection directive, #SELECT method = 1 (enter the number of the method you want to try).

There is a DEBUG line commented out within the core CORDIC loop. If you put that line back in the program, you can better visualize the internal operation of the CORDIC machinery as it rotates through smaller and smaller angles.

more or less precision, further developments |
---|

Things will become only a little more difficult when we move to double precision. But not too much so. It is still just shifts, no multiply or divide. We need to be able to compare magnitudes, and add and subtract. And it requires a double precision table of the arctangents and a new larger scaling for unity and 360 degrees.