CORDIC

Der CORDIC-Algorithmus (COordinate Rotation DIgital Computer) ist ein effizienter iterativer Algorithmus, mit dessen Hilfe sich viele Funktionen implementieren lassen, wie z.B.: trigonometrische, exponential und logarithmische sowie auch die einfache Multiplikation oder Division.

Motivation

In der Rechentechnik vornehmlich in der digitalen Signalverarbeitung benötigt man schnelle Verfahren für die Berechnung von bspw. trigonometrischen Funktionen. Herkömmliche Reihenentwicklungen wie z.B. die Taylorreihe zeigen oft nur mittelmäßige (d.h. langsame, oder gar von den Argumenten abhängige) Konvergenz und schlechte numerische Stablilität. Eine Reihenentwicklung besteht außerdem hauptsächlich aus einer Summe von Produkten. Da Muliplizierwerke relativ groß und damit teuer sind, wird CORDIC oft genau da eingesetzt wo Multiplizierer fehlen.

Anwendung

Cordic-Algorithmen werden zur Berechnung der wichtigsten Elementarfunktionen in Mikrocontroller-Rechenwerken eingesetzt und in der Nachrichtenübertragung. Mit ihnen lassen sich dort beispielsweise Betrag und Phase eines komplexen Signals bestimmen.

Funktionsweise (Zweidimensional)

CORDIC kann man im $\mathbb {~} ^{R^{3}}$ , aber auch nur in der zweidimensionalem Ebene betrachten.

Dreht man ein Koordinatensystem um den Winkel $-\Theta$ , erscheint der Vektor $(1,\;0)^{T}$ um den Winkel $\Theta$ gedreht; sein Endpunkt liegt im neuen System bei $y=\sin \,\Theta$ und $x=\cos \,\Theta$ .

Die Rotation um den Winkel $\Theta$ entspricht dem Matrix·Vektor Produkt:

${\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}}$

D.h. um auf den eigentlichen Funktionswert zu kommen, muss der Einheitsvektor $(1,\;0)^{T}$ um $\Theta$ gedreht werden. Dies läßt sich leichter bewerkstelligen, wenn innerhalb der Transformatrix nur noch eine Abhängigkeit von einer Winkelfunktion z.B. $\tan$ besteht:

${\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}}$

Die Drehung um $\Theta$ wird trickreich realisiert als Linearkombination von Teildrehungen um geschickt gewählte Teilwinkel $\alpha _{i}$ .

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

Eine zu weite Drehung im Schritt $i$ wird kompensiert durch einen Vorzeichenwechsel $\sigma _{i+1}:=-\sigma _{i}$ . Das gezeigte Verfahren konvergiert und ist numerisch stabil für alle $\Theta$ , die sich aus obiger Summe ergeben können. Man führt nun noch eine Hilfsvariable $z$ ein, die für den Drehsinn Verantwortung trägt:

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

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

Wenn nur einfachste Bauteile verwendet werden sollen und daher keine Multiplizierer vorhanden sind, muss man alles über Schiebe- und Addieroperationen bewerkstelligen. Dieses wird erreicht durch den Ansatz

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

Man erhält damit den folgenden Algorithmus:

${\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_{i}\\y_{i}\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_{i}\\y_{i}\end{pmatrix}}$

mit dem Skalierungsfaktor $K=\prod _{i=0}^{n-1}\cos \alpha _{i}$ , der während der Initialisierungsphase implizit berechnet wird.

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

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

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

Initialisierung

Vorweg wird eine Tabelle $T$ fester Länge $L$ angelegt mit $T_{1}=\pi /4=\arctan \,\delta _{1}$ wobei $\delta _{1}=1$ ist. Die folgenden Werte sind: $T_{j}=\arctan \,\delta _{j}$ mit $\delta _{j}=\delta _{j-1}/2$ . (Die Werte des Arcustangens lassen sich mit der hier gut konvergierenden Potenzreihenentwicklung bestimmen.)

Die Länge $L$ der Tabelle bestimmt die erreichbare Genauigkeit. Führt man alle Drehungen eines Einheitsvektors $(0,\;1)^{T}$ mit den so berechneten Werten hintereinander in gleichem Drehsinn aus, erzielt man eine Gesamtdrehung von etwas mehr als $\pm$ $\pi /2$ . Der Skalenfaktor $K$ wird mit einem Aufruf im Vektormodus (s.u.) berechnet, indem man die Verlängerung des Einheitsvektors $(0,\;1)^{T}$ ohne Skalierung berechnet.

Rotationsmodus

Der Ausgangsvektor $(1,\;0)^{T}$ wird in jedem der Schritte so gedreht, dass der Winkel $z=\Theta$ gegen Null geht. Es werden stets alle $L$ Teildrehungen ausgeführt, mit ggfls. wechselndem Vorzeichen. Da der Kosinus eine gerade Funktion ist, spielt das Vorzeichen bei der Skalierung keine Rolle. Nach Reskalierung sind die Komponenten des erhaltenen Endvektors $\cos \,\Theta$ und $\sin \,\Theta$ . Der Konvergenzbereich ergibt sich zu $-\sum _{i=0}^{n-1}\;\arctan \,2^{-i}\;\leq \Theta \leq \;+\sum _{i=0}^{n-1}\;\arctan \,2^{-i}$ , also bei genügend großem $n$ etwa zu $-1,74\;\leq \theta \leq \;+1,74$ , d.h. er erstreckt sich über mehr als den vierten und ersten Quadranten.

Vektormodus

Der vorgegebene Vektor, dessen Polarkoordinaten gesucht werden, wird immer so gedreht, dass sich der Betrag seiner $y$ -Komponente verringert. Der Drehwinkel $\Theta$ wird dabei vorzeichenrichtig aufaddiert. Die $x$ -Komponente des Endvektors ist nach Reskalierung der Betrag des Ausgangsvektors. Dieser Modus wird auch benutzt zu Berechnung des Arcustangens aus zwei Argumenten, Start mit $c\,(\cos \,\Theta ,\quad \sin \,\Theta )^{T}$ . Der Konvergenzbereich ist derselbe wie oben. Aus $\operatorname {arctan2} \,(y_{0},\,x_{0})$ lassen sich die Funktionen $\arcsin \,y_{0}$ und $\arccos \,x_{0}$ unter Zuhilfenahme von $x_{0}^{2}+y_{0}^{2}=1$ leicht ableiten.

Bereich außerhalb von $\mp$ $\pi /2$

Der Startvektor $(0,\;1)^{T}$ bzw. $(0,\;-1)^{T}$ entspricht einer Vorwegdrehung von $\pi /2$ bzw. $-\pi /2$ (für den Rotationsmodus). Bei einem Startvektor mit negativer $x$ -Komponente im Vektormodus bewirkt man entsprechende Drehungen durch Vertauschen der Komponenten und Änderungen der Vorzeichen.

Verallgemeinerung

Die oben benutzten Iterationsformeln ${\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}}$ sind ein Sonderfall der allgemeineren Vorschrift

${\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}}$ mit $m=1$ und $\delta _{i}=2^{-i}$ sowie $t_{i}=\arctan \,\delta _{i}$ .

Lineare Modi

Für $m=0$ , $\delta _{i}=t_{i}=[1,\;1/2,\;1/4,\;1/8\;...]$ und $k=1$ erhält man

${\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}}$ , womit sich Multiplikation und Division durchführen lassen. Eine Tabelle $T$ erübrigt sich hier.

Multiplikation: $x_{0}=a$ , $z_{0}=b$ , ergibt im Rotationsmodus ( $z$ gegen 0) $y_{n}=y_{0}+a\cdot b$ für alle $-2<b<2$ .

Division: $x_{0}=b$ , $y_{0}=a$ , ergibt im Vektormodus ( $z$ gegen 0) $z_{n}=z_{0}+a/b$ für alle $-2<b<2$ .

Hyperbolische Modi

Mit $m=-1$ werden die Hyberbelfunktionen, ihre Umkehrungen (Area-Funktionen), Exponentialfunktion und Logarithmus sowie die Quadratwurzel berechenbar. Einheitskreis bzw. -Hyberbel werden durch ${\sqrt {x^{2}+my^{2}}}=1$ mit $m=+1$ bzw. $m=-1$ beschrieben. Das zu einem Vektor $(x,\;y)^{T}$ gehörende Winkel- bzw. Areaargument ist durch $\Theta =\arctan({\sqrt {m}}\cdot y/x)/{\sqrt {m}}$ gegeben, also

$m=1$ , Winkelfunktionen (s.o): $\Theta =\arctan \,(y/x)$ und

$m=-1$ , hyperbolische Fkt.: $\Theta =\arctan \,(iy/x)/i$ , hier $i^{2}=-1$ ; $i\Theta =\arctan \,(iy/x)$ ; und wegen $\operatorname {tanh} \,\alpha =-i\tan \,(i\alpha )$ auch $\Theta =-\operatorname {Atanh} \,(y/x)$ .

Das Verfahren ist ganz analog zu dem eingangs gezeigten für die Winkelfunktionen. Erforderlich sind nur eine weitere Tabelle mit $t_{i}=\operatorname {Atanh} \,\delta _{i}$ , $\delta _{i}=[1/2,\;1/4,\;1/8\;...]$ und die einmalige Berechnung des Skalenfaktors $K$ .

${\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}}$ Die Iterationen $i=$ 4, 13, 40, $k$ , $3k+1$ müssen immer wiederholt werden, da sonst $\sum \sigma _{i}t_{i}\leq 0$ nicht erreicht werden kann.

Rotation mit $[x_{0},\;y_{0}]=[1,\;0]$ liefert: $x_{n}=\operatorname {cosh} \,z_{0},\quad y_{n}=\operatorname {sinh} \,z_{0}$ ,

davon abgeleitet: $\operatorname {tanh} \,z=\operatorname {sinh} \,z/\operatorname {cosh} \,z$ und $e^{z}=\operatorname {cosh} \,z+\operatorname {sinh} \,z$

Vektormodus mit $z_{0}=0$ berechnet: $z_{n}=\operatorname {Atanh} \,(y_{0}/x_{0})$ und den hyperbolischen Betrag $x_{n}={\sqrt {x_{0}^{2}-y_{0}^{2}}}$

davon abgeleitet: $\operatorname {ln} \,w=2\,\operatorname {Atanh} \,{\frac {w-1}{w+1}}$ sowie ${\sqrt {w}}$ aus dem Betrag des Startvektors $[w+1/4,\;w-1/4]$

Der Konvergenzbereich ist in beiden Modi beschränkt durch die maximal mögliche Änderung von $z$ . Alle mathematisch erlaubten Argumente können jedoch durch einfache Umstellungen und Shift-Operationen auf ihn abgebildet werden.

Alternativen

Sind hauptsächlich schnelle Tablelookup-Verfahren wie sie in DSPs stattfinden.

References

Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, September 1959
Daggett, D. H., Decimal-Binary conversions in CORDIC, IRE Transactions on Electronic Computers, Vol. EC-8 #5, pp335-339, IRE, September 1959.
J. E. Meggitt, Pseudo Division and Pseudo Multiplication Processes, 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
Schmid, Hermann, Decimal computation. New York, Wiley, 1974
V.D.Baykov,V.B.Smolov, Hardware implementation of elementary functions in computers, Leningrad State University, 1975, 96p.
Senzig, Don, Calculator Algorithms, IEEE Compcon Reader Digest, IEEE Catalog No. 75 CH 0920-9C, pp139-141, IEEE, 1975.
V.D.Baykov,S.A.Seljutin, Elementary functions evaluation in microcalculators, Moscow, Radio & svjaz,1982,64p.
Vladimir D.Baykov, Vladimir B.Smolov, Special-purpose processors: iterative algorithms and structures, Moscow, Radio & svjaz, 1985, 288 pages
M. E. Frerking, Digital Signal Processing in Communication Systems, 1994
Vitit Kantabutra, On hardware for computing exponential and trigonometric functions, IEEE Trans. Computers 45 (3), 328-339 (1996).
Andraka, Ray, A survey of CORDIC algorithms for FPGA based computers

Weblinks

[1]