CORDIC

The CORDIC algorithm ( English co ordinate R otation Di gital C omputer ) is an efficient iterative algorithm , many with the aid of functions can be implemented such. B. trigonometric functions , exponential functions and logarithms as well as simple multiplication or division .

motivation

In computer technology, especially in digital signal processing , fast methods are required for the calculation of, for example, trigonometric functions . Conventional series developments such as B. the Taylor series often show only mediocre (i.e. slow, or even argument-dependent) convergence and poor numerical stability . A series expansion also consists mainly of a sum of products that are difficult to calculate.

historical development

The CORDIC algorithm was presented in 1959 by Jack E. Volder. So that in the original version it was possible trigonometric functions such as sine, cosine and tangent and the multiplication and division of numbers (Engl. Only by the easy to implement in digital circuits additions and shift operations shift-and-add operations ) to form. Shift operations to number base 2 can be implemented very easily in digital circuits through appropriate interconnection.

Volder's motivation was to replace the usual and error-prone analog navigation computers in Convair B-58 bombers with digital computers for precise position determination. The requirement was to calculate the position of the bombers flying at supersonic speed in real time over a surface of the earth assumed to be spherical.

In the mid-1960s, the CORDIC algorithm was also used in civil applications. Forerunner of today's calculators like the desktop computer 9100 by Hewlett-Packard from 1968 it put an to calculate the trigonometric functions.

It was not until 1971, however, that J. S. Walther's CORDIC algorithm was expanded to the form that is common today, making it possible to efficiently calculate logarithms , the exponential function and the square root in digital circuits.

Application examples

CORDIC digital circuit

CORDIC algorithms are used to calculate the most important elementary functions in microcontroller arithmetic units such as pocket calculators. The CORDIC algorithm for calculating mathematical operations can also be found in Intel x87 arithmetic coprocessors . Further application examples are in message transmission. This can be used to efficiently determine the amount and phase of a complex signal, for example.

Since multipliers are extensive and therefore expensive to implement, especially in digital circuits, CORDIC is often used precisely where multipliers are not efficiently available. This mainly includes the area of ​​digital circuit technologies such as FPGAs or ASICs .

While CORDIC is not the fastest algorithm, it is widely used because of its simplicity and versatility.

functionality

Illustration by CORDIC

CORDIC can be viewed in , but also only in the two-dimensional plane. In the following, the description covers the simpler, two-dimensional case. ${\displaystyle \mathbb {R} ^{3}}$

If you rotate a coordinate system by the angle , the vector appears rotated by the angle ; its end point in the new system is at and . ${\displaystyle \Theta }$ ${\displaystyle (1,\;0)^{T}}$ ${\displaystyle \Theta }$ ${\displaystyle x=\cos \Theta }$ ${\displaystyle y=\sin \Theta }$

The rotation around the angle corresponds to the matrix-vector product : ${\displaystyle \Theta }$

${\displaystyle {\begin{pmatrix}x_{n}\\y_{n}\end{pmatrix}}={\begin{pmatrix}\cos \Theta &-\sin \Theta \\\sin \Theta &\cos \Theta \end{pmatrix}}\cdot {\begin{pmatrix}1\\0\end{pmatrix}}}$

In other words, in order to get to the actual function value, the unit vector must be reversed. This can be done more easily if, within the transformation matrix, there is only a dependency on an angle function, e.g. B. consists of: ${\displaystyle (1,\;0)^{T}}$ ${\displaystyle \Theta }$ ${\displaystyle \tan }$

${\displaystyle {\begin{pmatrix}x_{n}\\y_{n}\end{pmatrix}}=\cos \Theta \cdot {\begin{pmatrix}1&-\tan \Theta \\\tan \Theta &1\end{pmatrix}}\cdot {\begin{pmatrix}1\\0\end{pmatrix}}}$

The rotation by is implemented in a tricky way as a linear combination of partial rotations by cleverly selected partial angles . ${\displaystyle \Theta }$ ${\displaystyle \alpha _{i}}$

${\displaystyle \Theta =\sum _{i}\sigma _{i}\cdot \alpha _{i}\qquad \sigma _{i}\in \{-1,1\}}$

Too much rotation in the step is compensated for by changing the sign . The method shown converges and is numerically stable for all that can result from the above sum. An auxiliary variable is now introduced that is responsible for the sense of rotation: ${\displaystyle i}$ ${\displaystyle \sigma _{i+1}:=-\sigma _{i}}$ ${\displaystyle \Theta }$ ${\displaystyle z}$

${\displaystyle z_{0}=\Theta \qquad \quad z_{i}=z_{i-1}-\sigma _{i-1}\cdot \alpha _{i-1}}$

${\displaystyle \sigma _{i}={\begin{cases}-1&{\mbox{falls}}\quad z_{i}\leq 0\\1&{\mbox{sonst}}\end{cases}}}$

If only the simplest components are to be used and therefore no multipliers are available, everything has to be done with shift and add operations. This is achieved through the approach

${\displaystyle \tan \alpha _{i}=2^{-i}}$.

This gives the following algorithm:

${\displaystyle {\begin{pmatrix}x_{n}\\y_{n}\end{pmatrix}}=\prod _{i=0}^{n-1}\cos \alpha _{i}{\begin{pmatrix}1&-\sigma _{i}2^{-i}\\\sigma _{i}2^{-i}&1\end{pmatrix}}\cdot {\begin{pmatrix}x_{0}\\y_{0}\end{pmatrix}}=K\cdot \prod _{i=0}^{n-1}{\begin{pmatrix}1&-\sigma _{i}2^{-i}\\\sigma _{i}2^{-i}&1\end{pmatrix}}\cdot {\begin{pmatrix}x_{0}\\y_{0}\end{pmatrix}}}$

with the scaling factor 0.60725 ... for which is calculated implicitly during the initialization phase . ${\displaystyle K=\prod _{i=0}^{n-1}\cos \alpha _{i}\approx }$ ${\displaystyle n\to \infty }$

${\displaystyle x_{i+1}\,:=\,x_{i}-\sigma _{i}2^{-i}y_{i}}$

${\displaystyle y_{i+1}\,:=\,y_{i}+\sigma _{i}2^{-i}x_{i}}$

${\displaystyle z_{i+1}\,:=\,z_{i}-\sigma _{i}\arctan 2^{-i}}$

initialization

A table of fixed length is created beforehand with , where is. The following values ​​are: with . (The values ​​of the arctangent can be determined with the power series expansion, which converges well here.) ${\displaystyle T}$ ${\displaystyle L}$ ${\displaystyle T_{1}=\pi /4=\arctan \delta _{1}}$ ${\displaystyle \delta _{1}=1}$ ${\displaystyle T_{j}=\arctan \delta _{j}}$ ${\displaystyle \delta _{j}=\delta _{j-1}/2}$

The length of the table determines the accuracy that can be achieved. If you execute all rotations of a unit vector with the values ​​calculated in this way one after the other in the same direction of rotation, you achieve a total rotation of slightly more than . The scale factor is calculated with a call in vector mode (see below) by calculating the extension of the unit vector without scaling. ${\displaystyle L}$ ${\displaystyle (0,\;1)^{T}}$ ${\displaystyle \pm \pi /2}$ ${\displaystyle K}$ ${\displaystyle (0,1)^{T}}$

Rotation mode

The output vector is rotated in each of the steps so that the angle approaches zero. All partial rotations are always carried out, with a possibly changing sign . Since the cosine is an even function, the sign does not matter when scaling. After rescaling, the components of the final vector obtained are and . The convergence area results in , i.e. if it is large enough, approximately to , i.e. that is, it spans more than the fourth and first quadrants. ${\displaystyle (1,0)^{T}}$ ${\displaystyle z=\Theta }$ ${\displaystyle L}$ ${\displaystyle \cos \Theta }$ ${\displaystyle \sin \Theta }$ ${\displaystyle -\sum _{i=0}^{n-1}\arctan 2^{-i}\leq \Theta \leq +\sum _{i=0}^{n-1}\arctan 2^{-i}}$ ${\displaystyle n}$ ${\displaystyle -1{,}74\leq \Theta \leq +1{,}74}$

Vector mode

The specified vector, whose polar coordinates are searched for, is always rotated in such a way that the amount of its component is reduced. The angles of rotation are added with the correct sign. The component of the final vector is the output vector after rescaling of the amount. This mode is also used to calculate the arctangent from two arguments, start with . The convergence area is the same as above. The functions can be easily derived from and with the aid of . ${\displaystyle y}$ ${\displaystyle \Theta }$ ${\displaystyle x}$ ${\displaystyle c(\cos \Theta ,\sin \Theta )^{T}}$ ${\displaystyle \operatorname {arctan2} (y_{0},x_{0})}$ ${\displaystyle \arcsin y_{0}}$ ${\displaystyle \arccos x_{0}}$ ${\displaystyle x_{0}^{2}+y_{0}^{2}=1}$

Range outside of ± π / 2

The start vector or corresponds to a forward rotation of or (for the rotation mode ). In the case of a start vector with a negative component in vector mode, corresponding rotations are effected by exchanging the components and changing the signs. ${\displaystyle (0,1)^{T}}$ ${\displaystyle (0,-1)^{T}}$ ${\displaystyle \pi /2}$ ${\displaystyle -\pi /2}$ ${\displaystyle x}$

generalization

The iteration formulas used above

${\displaystyle {\begin{pmatrix}x_{i}\\y_{i}\end{pmatrix}}=\cos \Theta \cdot {\begin{pmatrix}1&-\tan \Theta \\\tan \Theta &1\end{pmatrix}}\cdot {\begin{pmatrix}x_{i-1}\\y_{i-1}\end{pmatrix}}}$

are a special case of the more general rule

${\displaystyle {\begin{pmatrix}x_{i}\\y_{i}\\z_{i}\end{pmatrix}}=k_{i}\cdot {\begin{pmatrix}1&-m\sigma _{i}\delta _{i}&0&0\\\sigma _{i}\delta _{i}&1&0&0\\0&0&1&-\sigma _{i}t_{i}\end{pmatrix}}\cdot {\begin{pmatrix}x_{i-1}\\y_{i-1}\\z_{i-1}\\1\end{pmatrix}}}$

with and as well . ${\displaystyle m=1}$ ${\displaystyle \delta _{i}=2^{-i}}$ ${\displaystyle t_{i}=\arctan \delta _{i}}$

Linear modes

For , and you get ${\displaystyle m=0}$ ${\displaystyle \delta _{i}=t_{i}=[1,1/2,1/4,1/8,\dotsc ]}$ ${\displaystyle k=1}$

${\displaystyle {\begin{pmatrix}x_{i}\\y_{i}\\z_{i}\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\\sigma _{i}\delta _{i}&1&0&0\\0&0&1&-\sigma _{i}\delta _{i}\end{pmatrix}}\cdot {\begin{pmatrix}x_{i-1}\\y_{i-1}\\z_{i-1}\\1\end{pmatrix}}}$,

with which multiplication and division can be carried out. A table is not necessary here. ${\displaystyle T}$

multiplication

${\displaystyle x_{0}=a}$, results in rotation mode ( towards 0) for all . ${\displaystyle z_{0}=b}$ ${\displaystyle z}$ ${\displaystyle y_{n}=y_{0}+a\cdot b}$ ${\displaystyle -2<b<2}$

division

${\displaystyle x_{0}=b}$, results in vector mode ( towards 0) for all . ${\displaystyle y_{0}=a}$ ${\displaystyle y}$ ${\displaystyle z_{n}=z_{0}+a/b}$ ${\displaystyle -2<a/b<2}$

Hyperbolic modes

With the hyperbolic functions, their inversions (area functions ), exponential function and logarithm as well as the square root can be calculated. The unit circle or hyperbola are described by with or . The angle or area argument belonging to a vector is given by, that is ${\displaystyle m=-1}$ ${\displaystyle {\sqrt {x^{2}+my^{2}}}=1}$ ${\displaystyle m=+1}$ ${\displaystyle m=-1}$ ${\displaystyle (x,y)^{T}}$ ${\displaystyle \Theta =\arctan({\sqrt {m}}\cdot y/x)/{\sqrt {m}}}$

${\displaystyle m=1}$, Trigonometric functions (see above): and ${\displaystyle \Theta =\arctan(y/x)}$

${\displaystyle m=-1}$, hyperbolic fct .:

${\displaystyle \Theta =\arctan(iy/x)/i}$

, here ; ; and because of too . ${\displaystyle i^{2}=-1}$ ${\displaystyle i\Theta =\arctan(iy/x)}$ ${\displaystyle \operatorname {tanh} \alpha =-i\tan(i\alpha )}$ ${\displaystyle \Theta =\operatorname {Atanh} (y/x)}$

The method is analogous to that shown at the beginning for the angle functions. All that is required is another table with , and the one-time calculation of the scale factor . ${\displaystyle t_{i}=\operatorname {Atanh} \delta _{i}}$ ${\displaystyle \delta _{i}=[1/2,1/4,1/8,\dotsc ]}$ ${\displaystyle K}$

${\displaystyle {\begin{pmatrix}x_{i}\\y_{i}\\z_{i}\end{pmatrix}}=k_{i}\cdot {\begin{pmatrix}1&\sigma _{i}\delta _{i}&0&0\\\sigma _{i}\delta _{i}&1&0&0\\0&0&1&-\sigma _{i}t_{i}\end{pmatrix}}\cdot {\begin{pmatrix}x_{i-1}\\y_{i-1}\\z_{i-1}\\1\end{pmatrix}}}$

The iterations must always be repeated, since the hyperbolic areatangens does not meet the condition that would not converge for the series . ${\displaystyle i=4,13,40,k,3k+1}$ ${\displaystyle t_{i+1}\leq {\tfrac {1}{2}}t_{i}}$ ${\displaystyle \sum \sigma _{i}t_{i}}$

Rotation with returns: _ , ${\displaystyle [x_{0},y_{0}]=[1,0]}$ ${\displaystyle x_{n}=\operatorname {cosh} z_{0},\quad y_{n}=\operatorname {sinh} z_{0}}$

derived from:

${\displaystyle \operatorname {tanh} z=\operatorname {sinh} z/\operatorname {cosh} z}$ and ${\displaystyle e^{z}=\operatorname {cosh} z+\operatorname {sinh} z}$

Vector mode with calculated: ${\displaystyle z_{0}=0}$

${\displaystyle z_{n}=\operatorname {Atanh} (y_{0}/x_{0})}$

and the hyperbolic amount

${\displaystyle x_{n}={\sqrt {x_{0}^{2}-y_{0}^{2}}}}$

derived from:

${\displaystyle \ln w=2\operatorname {Atanh} {\frac {w-1}{w+1}}}$

such as

${\displaystyle {\sqrt {w}}}$

from the amount of the start vector ${\displaystyle [w+1/4,w-1/4]}$

The convergence range is limited in both modes by the maximum possible change of . However, all mathematically permitted arguments can be mapped to it by simple adjustments and shift operations. ${\displaystyle z}$

Alternatives

The main alternatives are fast table look-up methods, such as B. in DSPs , and bit algorithms that use a similar approach to CORDIC to perform the calculation.

literature

• Jack E. Volder: The CORDIC Trigonometric Computing Technique. In: IRE Transactions on Electronic Computers . September 1959.
• DH Daggett: Decimal-Binary conversions in CORDIC. In: IRE Transactions on Electronic Computers. Vol. EC-8 # 5, pp 335-339, IRE, September 1959.
• JE Meggitt, Pseudo Division and Pseudo Multiplication Processes . In: IBM Journal. April 1962.
• Vladimir Baykov: Problems of Elementary Functions Evaluation Based on Digit by Digit (CORDIC) Technique . PhD thesis, Leningrad State Univ. of Electrical Eng., 1972.
• Hermann Schmid: Decimal computation. Wiley, New York 1974.
• VD Baykov, VB Smolov: Hardware implementation of elementary functions in computers . Leningrad State University, 1975, 96p.
• Don Senzig, Don: Calculator Algorithms . In: IEEE Compcon Reader Digest. IEEE Catalog No. 75 CH 0920-9C, pp 139-141, IEEE, 1975.
• VD Baykov, SA Seljutin: Elementary functions evaluation in microcalculators . Radio & svjaz, Moscow 1982, 64p.
• Vladimir D. Baykov, Vladimir B. Smolov: Special-purpose processors: iterative algorithms and structures. Radio & svjaz, Moscow 1985, 288 pages.
• ME Frerking: Digital Signal Processing in Communication Systems. 199
• Vitit Kantabutra: On hardware for computing exponential and trigonometric functions . In: IEEE Trans. Computers. 45 (3): 328-339 (1996).
• Ray Andraka: A survey of CORDIC algorithms for FPGA based computers.
• Jean-Michel Muller: Elementary Functions . Verlag Birkhäuser, 2006, ISBN 0-8176-4372-9 .

Web links

• CORDIC Bibliography Site
• http://www.andraka.com/cordic.php
• BASIC Stamp II math note #C
• VHDL implementation of the CORDIC algorithm on Opencores.org