A survey of CORDIC algorithms for FPGA based computers
Ray Andraka
Andraka Consulting Group, Inc
16 Arcadia Drive
North Kingstown, RI 02852
401/884-7930 FAX 401/884-7950
email:randraka@ids.net
1. ABSTRACT
The current trend back toward hardware
intensive signal processing has uncovered a
relative lack of understanding of hardware
signal processing architectures. Many
hardware efficient algorithms exist, but these
are generally not well known due to the
dominance of software systems over the past
quarter century. Among these algorithms is a
set of shift-add algorithms collectively known
as CORDIC for computing a wide range of
functions including certain trigonometric,
hyperbolic, linear and logarithmic functions.
While there are numerous articles covering
various aspects of CORDIC algorithms, very
few survey more than one or two, and even
fewer concentrate on implementation in
FPGAs. This paper attempts to survey
commonly used functions that may be
accomplished using a CORDIC architecture,
explain how the algorithms work, and explore
implementation specific to FPGAs.
1.1 Keywords
CORDIC, sine, cosine, vector magnitude, polar conversion
2. INTRODUCTION
The digital signal processing landscape has long been
dominated by microprocessors with enhancements such as
single cycle multiply-accumulate instructions and special
addressing modes. While these processors are low cost
and offer extreme flexiblility, they are often not fast enough
for truly demanding DSP tasks. The advent of
reconfigurable logic computers permits the higher speeds of
dedicated hardware solutions at costs that are competitive
with the traditional software approach. Unfortunately,
algorithms optimized for these microprocessor based
systems do not usually map well into hardware. While
hardware-efficient solutions often exist, the dominance of
the software systems has kept those solutions out of the
spotlight. Among these hardware-efficient algorithms is a
class of iterative solutions for trigonometric and other
transcendental functions that use only shifts and adds to
perform. The trigonometric functions are based on vector
rotations, while other functions such as square root are
implemented using an incremental expression of the desired
function. The trigonometric algorithm is called CORDIC,
an acronym for COordinate R
otation DIgital Computer.
The incremental functions are performed with a very simple
extension to the hardware architecture, and while not
CORDIC in the strict sense, are often included because of
the close similarity. The CORDIC algorithms generally
produce one additional bit of accuracy for each iteration.
The trigonometric CORDIC algorithms were originally
developed as a digital solution for real-time navigation
problems. The original work is credited to Jack Volder
[4,9]. Extensions to the CORDIC theory based on work by
John Walther[1] and others provide solutions to a broader
class of functions. The CORDIC algorithm has found its
way into diverse applications including the 8087 math
coprocessor[7], the HP-35 calculator, radar signal
processors[3] and robotics. CORDIC rotation has also been
proposed for computing Discrete Fourier[4], Discrete
Cosine[4], Discrete Hartley[10] and Chirp-Z [9] transforms,
filtering[4], Singular Value Decomposition[14], and solving
linear systems[1].
This paper attempts to survey the existing CORDIC and
CORDIC-like algorithms with an eye toward
implementation in Field Programmable Gate Arrays
(FPGAs). First a brief description of the theory behind the
algorithm and the derivation of several functions is
presented. Then the theory is extended to the so-called
unified CORDIC algorithms, after which implementation of
FPGA CORDIC processors is discussed.
Permission to make digital or hard copies of part or all of this work for
personal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that
copies bear this notice and the full citation on the first page. Copyrights
for components of this work owned by others than ACM must be
honored. Abstracting with credit is permitted. To copy otherwise, to
republish, to post on servers, or to redistribute to lists, requires prior
specific permission and/or a fee. Request permissions from Publications
Dept, ACM Inc., fax +1 (212) 869-0481, or permissions@acm.org.
FPGA 98 Monterey CA USA
Copyright 1998 ACM 0-89791-978-5/98/01..$5.00
3. CORDIC THEORY: AN ALGORITHM
FOR VECTOR ROTATION
All of the trigonometric functions can be computed or
derived from functions using vector rotations, as will be
discussed in the following sections. Vector rotation can
also be used for polar to rectangular and rectangular to
polar conversions, for vector magnitude, and as a building
block in certain transforms such as the DFT and DCT. The
CORDIC algorithm provides an iterative method of
performing vector rotations by arbitrary angles using only
shifts and adds. The algorithm, credited to Volder[4], is
derived from the general (Givens) rotation transform:
xx y
yy x
’cos sin
’ cos sin
=−
=+
φφ
φφ
which rotates a vector in a Cartesian plane by the angle
φ
.
These can be rearranged so that:
[]
[]
xxy
yyx
’cos tan
’cos tan
=⋅−
=⋅+
φφ
φφ
So far, nothing is simplified. However, if the rotation
angles are restricted so that tan(
φ)=±
2
-i
, the multiplication
by the tangent term is reduced to simple shift operation.
Arbitrary angles of rotation are obtainable by performing a
series of successively smaller elementary rotations. If the
decision at each iteration, i, is which direction to rotate
rather than whether or not to rotate, then the cos(
δ
i
) term
becomes a constant (because cos(
δ
i
) = cos(-
δ
i
)). The
iterative rotation can now be expressed as:
[]
[]
xKxyd
yKyxd
iiiii
i
iiiii
i
+
−
+
−
=−⋅⋅
=+⋅⋅
1
1
2
2
where:
()
K
d
i
ii
i
==+
=±
−− −
cos tan
12
2112
1
Removing the scale constant from the iterative equations
yields a shift-add algorithm for vector rotation. The
product of the K
i
’s can be applied elsewhere in the system
or treated as part of a system processing gain. That product
approaches 0.6073 as the number of iterations goes to
infinity. Therefore, the rotation algorithm has a gain, A
n
,
of approximately 1.647. The exact gain depends on the
number of iterations, and obeys the relation
A
n
i
n
=+
−
∏
12
2
The angle of a composite rotation is uniquely defined by the
sequence of the directions of the elementary rotations. That
sequence can be represented by a decision vector. The set
of all possible decision vectors is an angular measurement
system based on binary arctangents. Conversions between
this angular system and any other can be accomplished
using a look-up. A better conversion method uses an
additional adder-subtractor that accumulates the elementary
rotation angles at each iteration. The elementary angles can
be expressed in any convenient angular unit. Those angular
values are supplied by a small lookup table (one entry per
iteration) or are hardwired, depending on the
implementation. The angle accumulator adds a third
difference equation to the CORDIC algorithm:
()
zzd
iii
i
+
−−
=−⋅
1
1
2tan
Obviously, in cases where the angle is useful in the
arctangent base, this extra element is not needed.
The CORDIC rotator is normally operated in one of two
modes. The first, called rotation by Volder[4], rotates the
input vector by a specified angle (given as an argument).
The second mode, called vectoring, rotates the input vector
to the x axis while recording the angle required to make that
rotation.
In rotation mode, the angle accumulator is initialized with
the desired rotation angle. The rotation decision at each
iteration is made to diminish the magnitude of the residual
angle in the angle accumulator. The decision at each
iteration is therefore based on the sign of the residual angle
after each step. Naturally, if the input angle is already
expressed in the binary arctangent base, the angle
accumulator may be eliminated. For rotation mode, the
CORDIC equations are:
()
xxyd
yyxd
zzd
iiii
i
iiii
i
iii
i
+
−
+
−
+
−−
=−⋅⋅
=+⋅⋅
=−⋅
1
1
1
1
2
2
2tan
where
d
i
= -1 if z
i
< 0, +1 otherwise
which provides the following result:
[]
[]
xAx zy z
yAy zx z
z
A
nn
nn
n
n
i
n
=−
=+
=
=+
−
∏
0000
0000
2
0
12
cos sin
cos sin
In the vectoring mode, the CORDIC rotator rotates the
input vector through whatever angle is necessary to align
the result vector with the x axis. The result of the vectoring
operation is a rotation angle and the scaled magnitude of
the original vector (the x component of the result). The
vectoring function works by seeking to minimize the y
component of the residual vector at each rotation. The sign
of the residual y component is used to determine which
direction to rotate next. If the angle accumulator is
initialized with zero, it will contain the traversed angle at
the end of the iterations. In vectoring mode, the CORDIC
equations are:
()
xxyd
yyxd
zzd
iiii
i
iiii
i
iii
i
+
−
+
−
+
−−
=−⋅⋅
=+⋅⋅
=−⋅
1
1
1
1
2
2
2tan
where
d
i
= +1 if y
i
< 0, -1 otherwise.
Then:
xAxy
y
zz
y
x
A
nn
n
n
n
i
n
=+
=
=+
=+
−
−
∏
0
2
0
2
0
1
0
0
2
0
12
tan
The CORDIC rotation and vectoring algorithms as stated
are limited to rotation angles between -
π
/2 and
π
/2. This
limitation is due to the use of 2
0
for the tangent in the first
iteration. For composite rotation angles larger than
π
/2, an
additional rotation is required. Volder[4] describes an
initial rotation
±π
/2. This gives the correction iteration:
xdy
ydx
zzd
’
’
’
=− ⋅
=⋅
=+⋅
π
2
where d = +1 if y<0, -1 otherwise.
There is no growth for this initial rotation. Alternatively, an
initial rotation of either
π
or 0 can be made, avoiding the
reassignment of the x and y components to the rotator
elements. Again, there is no growth due to the initial
rotation:
xdx
ydy
’
’
=⋅
=⋅
z’ = z if d= 1, or z - π if d= -1
d = -1 if x<0, +1 otherwise.
Both reduction forms assume a modulo 2π representation of
the input angle. The style of first reduction is more
consistent with the succeeding rotations, while the second
reduction may be more convenient when wiring is
restricted, as is often the case with FPGAs.
The CORDIC rotator described is usable to compute
several trigonometric functions directly and others
indirectly. Judicious choice of initial values and modes
permits direct computation of sine, cosine, arctangent,
vector magnitude and transformations between polar and
Cartesian coordinates.
3.1 Sine and Cosine
The rotational mode CORDIC operation can simultaneously
compute the sine and cosine of the input angle. Setting the
y component of the input vector to zero reduces the rotation
mode result to:
xAx z
yAx z
nn
nn
=⋅
=⋅
00
00
cos
sin
By setting x
0
equal to 1/ A
n
, the rotation produces the
unscaled sine and cosine of the angle argument,
z
0
. Very
often, the sine and cosine values modulate a magnitude
value. Using other techniques (e.g., a look up table)
requires a pair of multipliers to obtain the modulation. The
CORDIC technique performs the multiply as part of the
rotation operation, and therefore eliminates the need for a
pair of explicit multipliers. The output of the CORDIC
rotator is scaled by the rotator gain. If the gain is not
acceptable, a single multiply by the reciprocal of the gain
constant placed before the CORDIC rotator will yield
unscaled results. It is worth noting that the hardware
complexity of the CORDIC rotator is approximately
equivalent to that of a single multiplier with the same word
size.
3.2 Polar to Rectangular Transformation
A logical extension to the sine and cosine computer is a
polar to Cartesian coordinate transformer. The
transformation from polar to Cartesian space is defined by:
x = rcos
θ
y = rsin
θ
As pointed out above, the multiplication by the magnitude
comes for free using the CORDIC rotator. The
transformation is accomplished by selecting the rotation
mode with
x
0
=polar magnitude, z
0
=polar phase, and y
0
=0.
The vector result represents the polar input transformed to
Cartesian space. The transform has a gain equal to the
rotator gain, which needs to be accounted for somewhere in
the system. If the gain is unacceptable, the polar magnitude
may be multiplied by the reciprocal of the rotator gain
before it is presented to the CORDIC rotator.
3.3 General vector rotation
The rotation mode CORDIC rotator is also useful for
performing general vector rotations, as are often
encountered in motion correction and control systems. For
general rotation, the 2 dimensional input vector is presented
to the rotator inputs. The rotator rotates the vector through
the desired angle. The output is scaled by the CORDIC
rotator gain, which must be accounted for elsewhere in the
system. If the scaling is unacceptable, a pair of constant
multipliers is required to compensate for the gain.
CORDIC rotators may be cascaded in a tree architecture for
general rotation in n-dimensions. Some optimization of
multidimensional rotation is possible to permit
computational savings over the general n-dimensioned case,
as reported by Hsiao et al. [4]
3.4 Arctangent
The arctangent,
θ
=Atan(y/x), is directly computed using
the vectoring mode CORDIC rotator if the angle
accumulator is initialized with zero. The argument must be
provided as a ratio expressed as a vector (x, y). Presenting
the argument as a ratio has the advantage of being able to
represent infinity (by setting x=0). Since the arctangent
result is taken from the angle accumulator, the CORDIC
rotator growth does not affect the result.
zz
y
x
n
=+
−
0
1
0
0
tan
3.5 Vector Magnitude
The vectoring mode CORDIC rotator produces the
magnitude of the input vector as a byproduct of computing
the arctangent. After the vectoring mode rotation, the
vector is aligned with the x axis. The magnitude of the
vector is therefore the same as the x component of the
rotated vector. This result is apparent in the result
equations for the vector mode rotator:
xAxy
nn
=+
0
2
0
2
The magnitude result is scaled by the processor gain, which
needs to be accounted for elsewhere in the system. This
implementation of vector magnitude has a hardware
complexity of roughly one multiplier of the same width.
The CORDIC implementation represents a significant
hardware savings over an equivalent Pythagorean
processor. The accuracy of the magnitude result improves
by 2 bits for each iteration performed.
3.6 Cartesian to Polar transformation
The Cartesian to Polar transformation consists of finding
the magnitude (r=sqrt(x
2
+y
2
)) and phase angle (
φ
=atan[y/x])
of the input vector, (x, y). The reader will immediately
recognize that both functions are provided simultaneously
by the vectoring mode CORDIC rotator. The magnitude of
the result will be scaled by the CORDIC rotator gain, and
should be accounted for elsewhere in the system. If the
gain is unacceptable, it can be corrected by multiplying the
resulting magnitude by the reciprocal of the gain constant.
3.7 Inverse CORDIC functions
In most cases, if a function can be generated by a CORDIC
style computer, its inverse can also be computed. Unless
the inverse is calculable by changing the mode of the
rotator, its computation normally involves comparing the
output to a target value. The CORDIC inverse is illustrated
by the Arcsine function.
3.8 Arcsine and Arccosine
The Arcsine can be computed by starting with a unit vector
on the positive x axis, then rotating it so that its y
component is equal to the input argument. The arcsine is
then the angle subtended to cause the y component of the
rotated vector to match the argument. The decision
function in this case is the result of a comparison between
the input value and the y component of the rotated vector at
each iteration:
()
xxyd
yyxd
zzd
iiii
i
iiii
i
iii
i
+
−
+
−
+
−−
=−⋅⋅
=+⋅⋅
=−⋅
1
1
1
1
2
2
2tan
where
d
i
= +1 if y
i
< c, -1 otherwise, and
c = input argument.
Rotation produces the following result:
()
xAxc
yc
zz
c
An x
A
nn
n
n
n
i
n
=⋅−
=
=+
⋅
=+
−
∏
0
2
2
0
0
2
12
arcsin
The arcsine function as stated above returns correct angles
for inputs -1 < c/A
n
x
0
< 1, although the accuracy suffers as
the input approaches
±
1 (the error increases rapidly for
inputs larger than about 0.98). This loss of accuracy is due
to the gain of the rotator. For angles near the y axis, the
rotator gain causes the rotated vector to be shorter than the
reference (input), so the decisions are made improperly.
The gain problems can be corrected using a “double
iteration algorithm”[9] at the cost of an increase in
complexity.
The Arccosine computation is similar, except the difference
between the x component and the input is used as the
decision function. Without modification, the arccosine
algorithm works only for inputs less than 1/A
n
, making the
double iteration algorithm a necessity. The Arccosine
could also be computed by using the arcsine function and
subtracting π/2 from the result, followed by an angular
reduction if the result is in the fourth quadrant.
3.9 Extension to Linear functions
A simple modification to the CORDIC equation permits the
computation of linear functions:
()
xx yd x
yyxd
zzd
ii ii
i
i
iiii
i
iii
i
+
−
+
−
+
−
=−⋅⋅⋅ =
=+⋅⋅
=−⋅
1
1
1
02
2
2
For rotation mode (d
i
= -1 if z
i
< 0, +1 otherwise) the linear
rotation produces:
xx
yyxz
z
n
n
n
=
=+
=
0
000
0
This operation is similar to the shift-add implementation of
a multiplier, and as multipliers go is not an optimal
solution. The multiplication is handy in applications where
a CORDIC structure is already available. The vectoring
mode (
d
i
= +1 if y
i
< 0, -1 otherwise) is more interesting, as
it provides a method for evaluating ratios:
xx
y
zzyx
n
n
n
=
=
=−
0
000
0
The rotations in the linear coordinate system have a unity
gain, so no scaling corrections are required.
3.10 Extension to Hyperbolic Functions
The close relationship between the trigonometric and
hyperbolic functions suggests the same architecture can be
used to compute the hyperbolic functions. While, there is
early mention of using the CORDIC structure for
hyperbolic coordinate transforms [4], the first description of
the algorithm is that by Walther [1]. The CORDIC
equations for hyperbolic rotations are derived using the
same manipulations as those used to derive the rotation in
the circular coordinate system. For rotation mode these are:
()
xxyd
yyxd
zzd
iiii
i
iiii
i
iii
i
+
−
+
−
+
−−
=+⋅⋅
=+⋅⋅
=−⋅
1
1
1
1
2
2
2tanh
where
d
i
= -1 if z
i
< 0, +1 otherwise.
Then:
[]
[]
xAx zy z
yAy zx z
z
A
nn
nn
n
n
i
n
=+
=+
=
=−≈
−
∏
0000
0000
2
0
12 080
cosh sinh
cosh sinh
.
In vectoring mode (d
i
= +1 if y
i
< 0, -1 otherwise) the
rotation produces:
xAxy
y
zz
y
x
A
nn
n
n
n
i
n
=−
=
=+
=−
−
−
∏
0
2
0
2
0
1
0
0
2
0
12
tanh
The elemental rotations in the hyperbolic coordinate system
do not converge. However, it can be shown[1] that
convergence is achieved if certain iterations (I=4, 13, 40,...,
k, 3k+1,...) are repeated.
The hyperbolic equivalents of all the functions discussed
for the circular coordinate system can be computed in a
similar fashion. Additionally, as Walther[1] points out, the
following functions can be derived from the CORDIC
functions:
tan
α
= sin
α/
cos
α
tanh
α
= sinh
α/
cosh
α
exp
α
= sinh
α
+ cosh
α
ln
α
= 2tanh
-1
[y/x] where x=
α
+1 and y=
α
-1
(
α
)
1/2
= (x
2
-y
2
)
1/2
where x=
α
+1/4 and y=
α
-1/4
It is worth noting the similarities between the CORDIC
equations for circular, linear, and hyperbolic systems. The
selection of coordinate system can be made by introducing
a mode variable that takes on values 1,0, or -1 for circular,
linear and hyperbolic systems respectively. The unified [1]
CORDIC iteration equations are then:
xxmyd
yyxd
zzde
ii ii
i
iiii
i
iiii
+
−
+
−
+
=−⋅⋅⋅
=+⋅⋅
=−⋅
1
1
1
2
2
where e
i
is the elementary angle of rotation for iteration i in
the selected coordinate system. Specifically,
e
i
= tan
-1
(2
-i
)
for m=1, e
i
= 2
-i
for m=0, and e
i
= tanh
-1
(2
-i
) for m=-1.
This unification, due to Walther, permits the design of a
general purpose CORDIC processor.
