CORDIC-Algorithmus zur Auswertung elementarer Funktionen in Hardware

• Verfahren • Implementierung in VHDL

FH-Report

Marcel Beuler

II

Dipl.-Ing.(FH) Marcel Beuler Fachhochschule Gießen-Friedberg FB Elektro- und Informationstechnik Wiesenstraße 14 35390 Gießen

CORDIC-Algorithmus zur Auswertung elementarer Funktionen in Hardware

Version 1.2

Juni 2008

III

Inhaltsverzeichnis 1 Einführung......................................................1

2 Der CORDIC-Algorithmus..........................................2

2.1 Herleitung.................................................2 2.2 Betriebsarten des CORDIC-Algorithmus.......................5 2.3 Verallgemeinerung des CORDIC-Algorithmus...................7 2.3.1 Hyperbolischer Modus (m=-1).........................10 2.3.2 Linearer Modus (m=0)................................13 2.3.3 Berechnung der Exponentialfunktion..................14 2.3.4 Erweiterung des zulässigen Argument-Bereiches für artanh(), ln() und sqrt().......................15 2.3.5 Zusammenfassung CORDIC-Algorithmus..................18

3 Prinzipielle Möglichkeiten der Implementierung.................19

3.1 Bitparallele iterative CORDIC-Struktur....................19 3.2 Bitserielle iterative CORDIC-Struktur.....................19 3.3 Bitparallele CORDIC-Struktur, rein kombinatorisch.........20 3.4 Bitparallele CORDIC-Struktur mit Pipeline-Registern.......21

4 Implementierung in VHDL........................................22

4.1 Funktionaler Überblick....................................22 4.2 Vorstufe..................................................23 4.3 Erweiterung des Argument-Bereiches der exp-Funktion.......25 4.4 CORDIC-Iterationen........................................27 4.5 Skalierung – Konstantenmultiplizierer.....................30 4.6 Zusätzliche Operationen für ln-Berechnung.................34 4.7 Zusätzliche Operationen für sqrt-Berechnung...............35 4.8 Endstufe..................................................35 4.9 Simulation................................................37

A Anhang.........................................................40

A1 C-Quellcode zur Berechnung von sin() und cos()............40 A2 C-Quellcode zur Berechnung von arctan()...................40 A3 C-Quellcode zur Berechnung von arcsin()...................41 A4 C-Quellcode zur Berechnung von sinh() und cosh()..........41 A5 C-Quellcode zur Berechnung von sinh() und cosh()..........42 A6 C-Quellcode zur Berechnung von artanh()...................42 A7 C-Quellcode zur Multiplikation im linearen Modus..........43 A8 C-Quellcode zur Division im linearen Modus................43 A9 C-Quellcode zur Berechnung von exp()......................44 A10 VHDL-Quellcode für CORDIC-Algorithmus.....................45

Literaturverzeichnis..............................................66

1

1 Einführung Taylor-Polynome werden traditionell als die erste und häufig auch alleinige Möglichkeit dargestellt, elementare Funktionen auszuwerten [Riss04]. Dabei zeigen schon einfachste, ältere Taschenrechner, dass eine Potenzreihen-Entwicklung eher nicht die Methode dieser Berechnung sein kann. Eigenschaften von Intel x86 (Co-) Prozessoren deuten in dieselbe Richtung. Zunächst sei an die Taylor-Reihe einer Exponentialfunktion erinnert. Nach [Papu01] gilt hierfür:

1 2 3 4

0

11! 2! 3! 4! !

nx

n

x x x x xe

n

∞

== + + + + + =∑… (1.1)

Um zu zeigen, dass die Konvergenz-Geschwindigkeit dieser Reihe vom Argument abhängig ist, betrachten wir folgendes Zahlenbeispiel (die Taylor-Reihe wird bei n=8 abgebrochen):

0.5

5.7

1.6487212707 (Taschenrechner)

e 298.8674009671 (Taschenrechner)

e =

=

8

0.5

0

0.51.6487212650

!

n

n

en=

≈ =∑

8

5.7

0

5.7261.9926451175

!

n

n

en=

≈ =∑

Bei einer Hardware-Realisierung könnten die Nenner-Fakultäten als Konstanten in einem Speicher (Look-Up Table) abgelegt werden. Dann hat man immer noch das Problem, dass mit größer werdendem x immer mehr Glieder zur Erzielung einer vorgegebenen Genauigkeit erforderlich sind. Um die Berechnungszeit in Grenzen zu halten, müssten die Potenzen durch den Einsatz vieler Multiplizierer parallel ermittelt werden, und deren Umsetzung ist in digitalen Schaltungen bekanntlich aufwändig. Es wird deutlich, dass die Berechnung elementarer Funktionen mittels Taylor-Reihe für eine Hardware-Realisierung aus Aufwandsgründen und schlechter Abschätzung des Konvergenzverhaltens (Berechnungszeit der Exponentialfunktion würde mit dem Argument schwanken) nicht in Frage kommt. Abhilfe schaffen iterative Verfahren wie CORDIC- und Bitalgorithmus, die auf einem Tripel (x,y,z) bzw. einem Tupel (x,y) operieren. Diese Algorithmen werden üblicherweise in der ingenieurwissenschaftlichen Ausbildung nicht behandelt, obwohl sie in vielen Anwendungen wie beispielsweise Signal- oder Bildverarbeitung eine wichtige Rolle spielen. Grund genug, sich mit diesen Verfahren auseinanderzusetzen! Eine ganz besondere Bedeutung hat der hier behandelte CORDIC-Algorithmus erlangt. Er ist zwar nach [WEB01] nicht der schnellste Algorithmus, wird aber aufgrund des Bestrebens nach einer regulär und einfach strukturierten Hardware oft eingesetzt. Durch Änderung eines Parameters können hier trigonometrische und hyperbolische Funktionen sowie die Exponential-, Logarithmus- und die Wurzelfunktion berechnet werden.

2

β

( , )n nx y

0 0( , )x yϕ

y

x

r

r

β

( , )n nx y

0 0( , )x yϕ

y

x

r

r

2 Der CORDIC-Algorithmus Der CORDIC-Algorithmus wurde 1959 von J. E. Volder entwickelt und 1971 von J. S. Walther auf die heute übliche Form erweitert. CORDIC ist die Abkürzung für Coordinate Rotation Digital Computer [Pirs96]. Dieser Algorithmus ermöglicht es, elementare trigonometrische und hyperbolische Funktionen unter fast ausschließlicher Verwendung von schnellen Operationen effizient zu berechnen. Mit schnellen Operationen sind hier Additionen (adds) und Multiplikationen mit Zweier-Potenzen (shifts) gemeint. 2.1 Herleitung Ausgangspunkt für die Herleitung des CORDIC-Algorithmus ist die Drehung eines Vektors [x0,y0]

T um den Winkel φ in einem kartesischen Koordinatensystem, siehe Bild 2.1. Durch die Drehung wird der Betrag des Vektors nicht verändert.

Bild 2.1: Vektordrehung im kartesischen Koordinatensystem

00

00

cos( ) cos( )

sin( ) sin( )

xx r

ry

y rr

β β

β β

= ⇒ = ⋅

= ⇒ = ⋅

Unter Anwendung der Additionstheoreme nach [Papu01] ergibt sich:

0 0

cos( ) cos( ) cos( ) sin( ) sin( )nx y

x r r rβ ϕ β ϕ β ϕ= ⋅ + = ⋅ ⋅ − ⋅ ⋅����� �����

0 0cos( ) sin( )nx x yϕ ϕ= ⋅ − ⋅ (2.1a)

0 0

sin( ) sin( ) cos( ) cos( ) sin( )ny x

y r r rβ ϕ β ϕ β ϕ= ⋅ + = ⋅ ⋅ + ⋅ ⋅����� �����

0 0sin( ) cos( )ny x yϕ ϕ= ⋅ + ⋅ (2.1b) Die Gln. (2.1a) und (2.1b) lassen sich zur sog. Drehmatrix zusammenfassen. Der Ergebnisvektor [xn,yn]

T berechnet sich demnach wie folgt:

0

0

cos( ) sin( )

sin( ) cos( )n

n

x x

y y

ϕ ϕϕ ϕ

− = ⋅

(2.2)

Nach [Papu01] gilt die Identität:

3

2

1cos( )

1 tan ( )ϕ

ϕ=

+ (2.3)

Durch Ausklammern von cos(φ) lässt sich die Beziehung nach Gl. (2.2) derart modifizieren, dass innerhalb der Drehmatrix nur noch eine Abhängigkeit von einer Winkelfunktion besteht:

0 020 0

1 tan( ) 1 tan( )1cos( )

tan( ) 1 tan( ) 11 tan ( )

n

n

x x x

y y y

ϕ ϕϕ

ϕ ϕϕ

− − = ⋅ ⋅ = ⋅ ⋅

+ (2.4)

Bei dem CORDIC-Algorithmus wird die Drehung um den Winkel φ durch mehrere Rotationen mit bekanntem Teilwinkel αi realisiert [Pirs96]:

1

0

{ 1,1}n

i i ii

ϕ σ α σ−

== ⋅ ∈ −∑ (2.5)

Durch festzulegende Vorzeichen σi wird die Summe der n Teilwinkel αi den gegebenen Winkel φ annähern. Um n für eine vorgegebene Genauigkeit möglichst gering zu halten, wird der Betrag der Teilwinkel αi mit größerem Index abnehmen. Sollen nur einfachste Bauelemente zum Einsatz kommen, dann kann dies durch eine Folge von Zweierpotenzen erreicht werden:

tan( ) 2 iiα−= (2.6)

Weiterhin wird nach [Pirs96] eine Hilfsvariable zi eingeführt, die die aufakkumulierten Teilwinkel enthält und sich zur Steuerung der Vorzeichen der Teilwinkel verwenden lässt (Rotating-Mode):

1 arctan(2 ) sgn( ) arctan(2 )i i

i i i i iz z z zσ− −

+ = − ⋅ = − ⋅ (2.7)

mit: 1 0

sgn( )1 0

ii i

i

zz

zσ

+ ≥= = −

4

( )

1

1

1 0

1 2 2

2 1 2

0 sgn( ) arctan 2

i ii i i i i i

i ii ii i i i

ii i i

x x x y

y y y x

z z z z

σ σσ σ

ϕ

− −+

− −+

−+

− ⋅ − ⋅ ⋅ = ⋅ =

⋅ + ⋅ ⋅

← = − ⋅ =

(2.12)

Die in [Riss04] angegebenen SHIFT- und ADD-Operationen berücksichtigen jedoch noch nicht den zu bestimmenden Arkustangens der Hilfsvariablen zi. Hierzu wird nach [WEB01] vor der eigentlichen Berechnung eine Tabelle T fester Länge L angelegt, wobei gilt:

1 1 1arctan( ) mit: 14T

πδ δ= = = (2.13a)

1arctan( ) mit: 2 ; 2j j j jT j Lδ δ δ −= = ≤ ≤ (2.13b) Die Länge L der Tabelle bestimmt die erreichbare Genauigkeit. Führt man alle CORDIC-Rotationen mit den so berechneten Tabellenwerten hintereinander im gleichen Drehsinn aus, so erzielt man nach Gl. (2.14) eine Gesamtdrehung von etwas mehr als ±π/2. Der Konvergenzbereich des Verfahrens erstreckt sich somit über mehr als den vierten und ersten Quadranten.

( )1

0

arctan 2 1.74n

i

i

ϕ−

−

== ± ≈ ±∑ (2.14)

Beispiel 2.1: Der Vektor [x0,y0]

T=[4,1]T soll um φ=52°=0,907571 gedreht werden (vgl. mit Bild 2.1). Mittels CORDIC-Algorithmus ist der skalierte Ergebnisvektor [x8,y8]

T mit K nach Gl. (2.11) zu bestimmen und mit der exakten Lösung zu vergleichen. Da φ>0 ist, ist das Start-Vorzeichen σ0 positiv. Mit Gl. (2.12) ergeben sich folgende Iterationsschritte:

( )

01 0 0 0

01 0 0 0

01 0 0 1

2 3

2 5

arctan 2 0.122173 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

( )

12 1 1 1

12 1 1 1

12 1 1 2

2 0.5

2 6.5

arctan 2 0.341475 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

23 2 2 2

23 2 2 2

23 2 2 3

2 2.125

2 6.375

arctan 2 0.096496 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

34 3 3 3

34 3 3 3

34 3 3 4

2 2.92188

2 6.10938

arctan 2 0.027859 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

( )

45 4 4 4

45 4 4 4

45 4 4 5

2 2.54004

2 6.292

arctan 2 0.03456 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

56 5 5 5

56 5 5 5

56 5 5 6

2 2.73667

2 6.21262

arctan 2 0.00332 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

67 6 6 6

67 6 6 6

67 6 6 7

2 2.83374

2 6.16986

arctan 2 0.012304 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

( )

78 7 7 7

78 7 7 7

78 7 7 8

2 2.78554

2 6.192

arctan 2 0.004492 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

Skalierungsfaktor nach Gl. (2.11): 8

0.607259n

K = = ⇒ 8 8

8 8

K 1.69154

3.76015

x x

K y y

⋅ → =⋅ → =

Für den exakten Ergebnisvektor [ , ]Tn nx y gilt nach Gl. (2.2): 1.67464

3.7677n

n

x

y

=

5

Beispiel 2.2: Mit dem CORDIC-Algorithmus sollen die trigonometrischen Funktionen sin(φ) und cos(φ) für φ im Bogenmaß mit |φ|

6

xK e⋅1xe =nx

1r = argny =

xK e⋅1xe =nx

1r = argny =

( )

1

1

1

2

2

arctan 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= − ⋅ ⋅

= + ⋅ ⋅

= − ⋅

mit: 01 0

0 ; 1 0

ii

i

yz

yσ

+

(2.15)

Daraus ergibt sich schließlich:

�

1 2 20 0

0 00

0 00

0

arctan arctan

n

n

n

x K x y

y

y yz z

x x

−= ⋅ +=

= + =

(2.16)

Die Drehrichtung des Vektors hängt vom Quadranten ab. Mit der in Gl. (2.15) angegebenen Vorzeichensteuerung lässt sich arctan(y0/x0) für den 1. und 4. Quadranten bestimmen. Der C-Quellcode im Anhang A2 berechnet wieder n=15 Iterationen mit den Initialwerten x=1.8 und y=0.7. Gemäß Taschenrechner gilt hierfür arctan(0.7/1.8)=0.370891. Die mittels Programm berechneten Ergebnisse lauten:

x = 1.93132079 // Betrag des Vektors z = 0.37086993 // arctan(0.7/1.8)

Arkussinusfunktion: Der Vectoring-Mode eignet sich nach [WEB02] und [Andr98] auch zur Berechnung der inversen trigonometrischen Funktion φ=arcsin(arg). Hierzu muss der skalierte Einheitsvektor K·ex einer Folge von Rotationen unterzogen werden, bis die y-Komponente des gedrehten Vektors mit arg übereinstimmt. Durch die anfängliche Skalierung wird erreicht, dass nach der Rotation der Betrag des Vektors gleich Eins ist, siehe Bild 2.3:

2 2 2 2 2n1 arg x 1 argn n nr x y x= = + = + ⇒ = − (2.17)

Bild 2.3: Berechnung von arcsin(arg) Mit der Hilfsvariablen zi werden nach Bild 2.3 alle Teilwinkel zum Winkel arctan(yn/xn) aufakkumuliert. Nach Gl. (2.17) kann man hierfür auch schreiben:

2

argarctan arctan

1 arg

n

n

y

x = −

Nach [Papu01] gilt zwischen dem Arkussinus und dem Arkustangens folgender Zusammenhang:

2 2

argarcsin( ) arctan arcsin(arg) arctan

1 1 arg

xx

x

= ⇒ = − − (2.18)

Somit entsprechen die aufakkumulierten Teilwinkel der inversen trigonometrischen Funktion φ=arcsin(arg). Es ergeben sich folgende Iterationsgleichungen:

7

( )

1

1

1

2

2

arctan 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= − ⋅ ⋅

= + ⋅ ⋅

= + ⋅

mit: 0 0 01 (arg ) 0

; 0 ; 1 (arg ) 0i

yx K y z

yσ

+ − >= = = = − − 0 ein Anwachsen von yi+1, d.h. es gilt hier yi+1>yi. Der im Anhang A3 dargestellte C-Quellcode berechnet wieder n=15 Iterationen mit dem Initialwert arg=0.7. Hierfür gilt gemäß Taschenrechner arcsin(0.7)=0.775397. Das Programm liefert folgende Ergebnisse: x = 0.71415968 y = 0.69998283 // arg = 0.7 z = 0.77537345 // arcsin(0.7) Das Dokument [WEB02] ermöglicht auch hier wieder eine interaktive Untersuchung. Die Ergebnisse für arg=0.7 nach n=15 Iterationen zeigt Bild 2.4, sie stimmen erwartungsgemäß mit denen des C-Quellcodes (gerundet) überein. [WEB02] berechnet zusätzlich noch die arccos-Funktion, für die folgender Zusammenhang gilt:

arccos(arg) arcsin(arg)2

π= − (2.20)

Bild 2.4: Berechnung von arcsin(0.7) und arccos(0.7) mittels interaktiven Dokuments [WEB02]

2.3 Verallgemeinerung des CORDIC-Algorithmus Der CORDIC-Algorithmus in seiner bisher definierten Form erlaubt mit den beiden Betriebsarten Rotating und Vectoring die Bestimmung der trigonometrischen Funktionen sin und cos (damit auch tan) sowie der Arkusfunktionen arctan und arcsin (damit auch arccos). Durch die von J. S. Walther formulierte Verallgemeinerung des Verfahrens können zusätzliche Elementarfunktionen wie beispielsweise die hyperbolischen Funktionen und die Exponentialfunktion ermittelt werden [Walt71]. Um die Verallgemeinerung des CORDIC-Algorithmus zu verstehen, betrachten wir die hyperbolischen Additionstheoreme:

0 0

0 0

cosh( ) cosh( ) cosh( ) sinh( ) sinh( )

sinh( ) sinh( ) cosh( ) cosh( ) sinh( )

x y

y x

β ϕ β ϕ β ϕ

β ϕ β ϕ β ϕ

± = ⋅ ± ⋅

± = ⋅ ± ⋅

����� �����

����� �����

(2.21)

Diese Additionstheoreme lassen sich wieder zu einer Drehmatrix zusammenfassen:

8

0 0

0 0

cosh( ) cosh( ) sinh( ) 1 tanh( )cosh( )

sinh( ) sinh( ) cosh( ) tanh( ) 1

x x

y y

β ϕ ϕ ϕ ϕϕ

β ϕ ϕ ϕ ϕ± ± ±

= ⋅ = ⋅ ⋅ ± ± ± (2.22)

Diese Drehmatrix nach [WEB02] legt es nahe, die Funktionen cosh(φ) und sinh(φ) für ein gegebenes Argument φ durch eine Folge von hyperbolischen Rotationen zeitgleich zu berechnen. Hierzu muss φ nach [WEB02] wieder als Summe geeigneter Dreh- bzw. Teilwinkel αi und Drehrichtungen σi dargestellt werden:

( )tanh( ) 2 bzw. artanh 2i ii iα α− −= = (2.23)

( )1 1

artanh 2 { 1,1}n n

ii i i i

i i

ϕ σ α σ σ−= =

= ⋅ = ⋅ ∈ −∑ ∑ (2.24)

Die Summation beginnt mit i=1, da artanh(0) nicht definiert ist. Somit lassen sich auch hier die einzelnen Rotationen – bis auf die Skalierung am Ende – ausschließlich mit SHIFT- und ADD-Operationen ausführen. Da der Kosinus hyperbolicus eine gerade Funktion ist, spielt das Vorzeichen σi bei der Bestimmung des Skalierungsfaktors keine Rolle:

( )1 1

cosh( ) cosh artanh 2n n

ii

i i

K α −= =

= = ∏ ∏ (2.25) Über den hyperbolischen Pythagoras kann nach [WEB02] folgende allgemeine Beziehung angegeben werden:

2 2

22

1 cosh ( ) sinh ( )1 tanh ( )

cosh( ) cosh ( )

x xx

x x

−= = −

Für den Skalierungsfaktor ergibt sich jetzt:

( ) ( )artanh 2 tanh( ) tanh artanh 2 2i i ix x− − − = ⇒ = =

2

1

1

1 2

n

ii

K− ⋅=

⇒ =−

∏ (2.26) Für n → ∞ liefert K folgenden Wert:

� 21

1lim 1.20514

1 2

n

in i

K− ⋅→∞ =

= ≈−

∏ (2.27) Mit β=0 erhält man: 0 00 0cosh( ) 1 ; sinh( ) 0x yβ ββ β= == = = = Die Hilfsvariable zi enthält wieder die aufakkumulierten Teilwinkel und steuert (im Rotating-Mode) das Vorzeichen σi=sgn(zi). Damit gilt für die Folge von hyperbolischen Rotationen zur Berechnung von cosh(φ) und sinh(φ):

( )

21 1

1 0

1 2 1 2cosh( ) 1 11

sinh( ) 0 02 1 2 11 2

0 sgn( ) artanh 2

i in ni i

i iii ii i

ii i i

K

z z z z

σ σϕϕ σ σ

− −

− −− ⋅= =

−+

⋅ ⋅ = ⋅ ⋅ = ⋅ ⋅

⋅ ⋅ −

← = − ⋅

∏ ∏

ϕ=

(2.28)

9

Gl. (2.28) zeigt, dass sich auf Grundlage der hyperbolischen Additionstheoreme die hyperbolischen Funktionen cosh(φ) und sinh(φ) mittels CORDIC-Verfahren ermitteln lassen. Nun sollen die Berechnungsvorschriften auf Basis der trigonometrischen und hyperbolischen Additionstheoreme gemeinsam durch eine allgemeine Iterationsgleichung beschrieben werden. Betrachten wird zunächst die Iterationsgleichungen für beide Fälle getrennt.

trigonometrisch (Gl. (2.12)):

( )

1

1

1

2

2

arctan 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= − ⋅ ⋅

= + ⋅ ⋅

= − ⋅

hyperbolisch (Gl. (2.28)):

( )

1

1

1

2

2

artanh 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= + ⋅ ⋅

= + ⋅ ⋅

= − ⋅

Wir führen nun nach [Pirs96] folgende Größen ein: {1,0, 1}m∈ − Parameter für den Modus zirkular (m=1), linear (m=0) und

hyperbolisch (m=-1)

, 2i

m iδ−= Schrittweite des i-ten Iterationsschrittes (0

10

2.3.1 Hyperbolischer Modus (m = -1) Rotating: Im letzten Abschnitt wurde durch Gl. (2.28) bereits gezeigt, dass mit dem Ausgangsvektor [x0,y0]

T=[1,0]T die hyperbolischen Funktionen cosh(φ) und sinh(φ) approximiert werden können. Der im Anhang A4 angegebene C-Quellcode führt dies mit dem Argument φ=0.8 für n=15 Iterationen aus. Gemäß Taschenrechner erhält man folgende Werte: cosh(0.8) 1.3374349463 ; sinh(0.8) 0.8881059822= = Die berechneten Ergebnisse nach Anhang A4 lauten: x = 1.33741188 // cosh(0.8) y = 0.88807125 // sinh(0.8) Gemäß [Papu01] weisen die hyperbolischen Funktionen cosh und sinh einen Definitionsbereich von −∞ bis +∞ auf. Werden alle Mikrorotationen artanh(δ-1,i) nach Tabelle 2.1 mit dem gleichen Drehsinn aufaddiert, dann erhält man allerdings folgenden Konvergenzbereich des C-Quellcode nach A4 für φ:

( )1

artanh 2 1.05547n

i

i

ϕ −=

= ± ≈ ±∑ (2.30)

Das bedeutet, dass der in Gl. (2.30) angegebene Wert das maximal zulässige Argument für den Quellcode nach A4 darstellt. Im Zusammenhang mit der Exponentialfunktion wird in Abschnitt 2.3.3 gezeigt, wie mit unzulässigen Argumenten zu verfahren ist. Konvergenzproblem: Zur einfachen Realisierung des CORDIC-Verfahrens wurden bisher für die Iterationsschrittweiten δm,i Zweierpotenzen eingesetzt, wodurch sich die Matrixmultiplikation durch eine Shiftoperation implementieren lässt (siehe Tabelle 2.1):

( , ) , 2 ( , ) {0,1, , 1} mit: {1,0}S m i

m i S m i n mδ−= = − ∈… (2.31a)

( , ), 2 ( , ) {1,2, , } mit: { 1}S m i

m i S m i n mδ−= = ∈ −… (2.31b)

Die ganzzahlige Folge S(m,i) wird nach [Pirs96] als Shiftfolge des CORDIC-Algorithmus bezeichnet. Sie beeinflusst die Konvergenz, den Konvergenzbereich und den Skalierungsfaktor. Die aus ihr resultierenden Teilwinkel αm,i sollen zunächst einmal eine monoton fallende Folge bilden. Darüber hinaus muss jeder Teilwinkel einer Iteration durch eine nachfolgende Richtungsumkehr durch alle folgenden Teilwinkel bis auf einen Restfehler kompensiert werden können. Dieser ist durch den Teilwinkel der letzten Iteration definiert. Somit muss für die Winkel gelten:

1

, , , 11

{0,1, , 2} mit: {1,0}n

m i m j m nj i

i n mα α α−

−= +

− < ∈ − ∈∑ … (2.32a)

, , ,1

{1,2, , 1} mit: { 1}n

m i m j m nj i

i n mα α α= +

− < ∈ − ∈ −∑ … (2.32b)

Für eine Folge, in der der Nachfolgewinkel mindestens die Hälfte des vorhergehenden beträgt, ist Gl. (2.32) erfüllt:

11 2i iα α+ ≥ ⋅ (2.33) Bei der zirkularen Drehung ist diese Bedingung erfüllt, denn für i=3 gilt beispielsweise:

11

( ) ( )4 312arctan 2 0.062419 arctan 2 0.062177− −= > ⋅ =

Die Verhältnisse im linearen Fall werden in Abschnitt 2.3.2 behandelt. Im hyperbolischen Modus wird Gl. (2.33) für eine ganzzahlige Shiftfolge nach Gl. (2.31b) nicht erfüllt. Mit dem C-Quellcode nach A4 soll nun gezeigt werden, wie sich diese Verletzung der Konvergenzbedingung auf die Bestimmung der Funktionen cosh(0.01) und sinh(0.01) auswirkt. Tabelle 2.2 zeigt die exakten Werte sowie die Ergebnisse für eine steigende Anzahl von Iterationsschritten (10 Stellen nach dem Komma). Man sieht, dass das Verfahren mit der gewählten Shiftfolge bei einem kleinen Argument völlig versagt.

n x (cosh(0.01)) y (sinh(0.01)) 8 1.0011095622 0.0470666476 16 1.0009314586 0.0431715729 24 1.0009308025 0.0431563596 32 1.0009307999 0.0431563001 40 1.0009307999 0.0431562999

exakt: 1.0000500004 0.0100001667 Tabelle 2.2: Ergebnisse für 0.01ϕ = bei Nichteinhaltung der

Konvergenzbedingung

Um die Konvergenzbedingung im hyperbolischen Fall zu erfüllen, muss eine Winkeldrehung an bestimmten Stellen wieder vollständig rückgängig gemacht werden können. Dies wird durch die Wiederholung der Iterationen i=4,13,40,k,3k+1 erreicht, siehe auch [Andr98, Pirs96, WEB01]. Man erhält somit die komplexere Shiftfolge nach Gl. (2.34). ( 1, ) {1,2,3,4,4,5,6, ,13,13,14,15, ,40,40,41, , }S i n− = … … … (2.34) Da die Iterationen i=4,13,40,k,3k+1 doppelt ausgeführt werden, ändert sich auch der Skalierungsfaktor. Für ihn gilt jetzt:

4 13 40 121

2 2 2 21 4 13 40

1 1 1 11.207497

1 2 1 2 1 2 1 2i i i ii i i iK

− ⋅ − ⋅ − ⋅ − ⋅= = = =

= ⋅ ⋅ ⋅ ⋅ ≈ − − − − ∏ ∏ ∏ ∏ … (2.35)

Der C-Quellcode im Anhang A5 berücksichtigt die Shiftfolge nach Gl. (2.34) und den neuen Skalierungsfaktor g=K-1 nach Gl. (2.35). Tabelle 2.3 zeigt die berechneten Ergebnisse für das Argument φ=0.01. Das CORDIC-Verfahren konvergiert jetzt gegen die exakten Werte.

n x (cosh(0.01)) y (sinh(0.01)) 8 1.0000337989 0.0079055564 16 1.0000499596 0.0099960772 24 1.0000500007 0.0100001902 32 1.0000500004 0.0100001667

exakt: 1.0000500004 0.0100001667 Tabelle 2.3: Ergebnisse für 0.01ϕ = bei Einhaltung der Konvergenzbedingung

Der Konvergenzbereich für φ erhöht sich nun gegenüber Gl. (2.30) etwas:

( ) ( ) ( ) ( )4 13 401

artanh 2 artanh 2 artanh 2 artanh 2 1.118173n

i

i

ϕ − − − −=

= ± + + + + ≈ ±

∑ … (2.36)

Ist φ vom Betrag her größer als der in Gl. (2.36) angegebene Wert, so liefert der Quellcode nach A5 ebenfalls falsche Ergebnisse. Wie bereits erwähnt, wird in Abschnitt 2.3.3 gezeigt, wie mit solchen unzulässigen Argumenten im Zusammenhang mit der Exponentialfunktion zu verfahren ist.

12

Vectoring: Nach Abschluss der Iteration liegt der Ausgangsvektor wieder auf der positiven x-Achse (y→0). Die aufakkumulierten Teilwinkel repräsentieren jetzt den Winkel artanh(y0/x0) und die x-Koordinate den noch zu skalierenden hyperbolischen Betrag des Vektors. Die einzelnen Teilschritte entsprechen bis auf die Initialisierung von z0 und der Vorzeichenbestimmung σi gemäß [Andr98] dem Rotating-Mode nach Gl. (2.28):

( )

1

1

1

2

2

artanh 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= + ⋅ ⋅

= + ⋅ ⋅

= − ⋅

mit: 01 0

0 ; 1 0

ii

i

yz

yσ

+

(2.37)

Daraus ergibt sich:

�

1 2 2 1 2 20 0 0 0

0 00

0 00

0

artanh artanh

n

n

n

x K x m y K x y

y

y yz z

x x

− −= ⋅ + ⋅ = ⋅ −=

= + =

(2.38)

Die Vorzeichensteuerung nach Gl. (2.37) ermöglicht es, artanh(y0/x0) für den 1. und 4. Quadranten zu bestimmen. Gemäß [Papu01] hat die Funktion artanh einen Definitionsbereich von 1− bis 1+ und einen Wertebereich von −∞ bis +∞ . Beim CORDIC-Algorithmus ist dieser Wertebereich für φ auf den Konvergenzbereich nach Gl. (2.36) beschränkt. Der im Anhang A6 angegebene C-Quellcode berechnet n=15 Iterationen für die Werte x=3.9 und y=2.5. Der Taschenrechner liefert hierfür artanh(2.5/3.9)=0.759913. Die Ergebnisse nach A6 lauten: x = 2.99332591 // hyperbolischer Betrag z = 0.75993257 // artanh(2.5/3.9) Der maximal mögliche Quotient der Werte x0 und y0 für den Quellcode A6 beträgt:

0 0

0 0

artanh 1.118173 tanh( 1.118173) 0.806932y y

x xϕ = = ± ⇒ = ± = ±

(2.39)

Logarithmus naturalis: Über den Areatangens hyperbolicus lässt sich auch der natürliche Logarithmus mittels CORDIC-Verfahren berechnen:

( ) 12ln ln

12ln ln

1 1tanh ln ; ln ln ln

11

x x

x x

e e x x xx x x x

xx xe e

−

−

− − −= = = = = ⋅+++

1

ln 2 artanh1

xx

x

− ⇒ = ⋅ +

(2.40)

Quadratwurzel: Für r kann man auch schreiben:

( ) ( )2 2 2 21 1 1 1 1 1 1 14 4 2 16 2 16 2 2r r r r r r r r r r r+ − − = + + − + − = + = (2.41)

mit: 1 10 04 4 ; x r y r= + = −

13

2.3.2 Linearer Modus (m = 0) Der CORDIC-Algorithmus im linearen Modus multipliziert (Rotating) oder dividiert (Vectoring). Auch wenn diese Implementierung den “herkömmlichen“ Implementierungen für Multiplizierer bzw. Dividierer nicht ebenbürtig ist, so bietet sich hier doch eine attraktive Variante an, beispielsweise wenn es um die Einsparung zusätzlicher funktionaler Baugruppen geht. Im linearen Fall entspricht die Schrittweite δ0,i dem Teilwinkel α0,i im i-ten Iterationsschritt:

0, 0, 2 {0,1, , 1}i

i i i nα δ−= = ∈ −… (2.42)

Der Skalierungsfaktor K0 ist 1. Aus Gl. (2.42) ist ersichtlich, dass sich im linearen Fall das System am Rande des Konvergenzbereiches befindet, da die Gleichheitsbeziehung in Gl. (2.33) gilt. Für Gl. (2.29) erhält man jetzt:

1

1 0,

1 0,

1 0 0 0

1 0 0

0 0 11

ii

ii i i

ii i i

xx

yy

zz

σ δσ δ

+

+

+

= ⋅ −

(2.43)

Rotating: In dieser Betriebsart gilt für die einzelnen Iterationsschritte:

1

1

1

2

0 2

i i

ii i i i

ii i i

x x

y y x

z z

σ

σ

+−

+−

+

=

= + ⋅ ⋅

← = − ⋅

mit: 01 z 0

0 ; 1 z 0

ii

i

y σ+ >

= = − ∑ (2.45)

Der Betrag von z0 muss somit

14

Daraus ergibt sich:

�

0

0 00

0 00

0n

n

n

x x

y

y yz z

x x

==

= + =

für: 1

0

00

lim 2 2n

i

n i

y

x

−−

→∞ == >∑ (2.47)

Die Vorzeichensteuerung in Gl. (2.46) ermöglicht es, den Quotienten y0/x0 für den 1. und 4. Quadranten zu bestimmen. Der C-Quellcode im Anhang A8 führt die Division für die Initialwerte x=4.1 und y=-5.8 aus. Das Ergebnis nach n=24 Iterationen lautet: z = -1.41463411 // y/x 2.3.3 Berechnung der Exponentialfunktion Nach [Papu01] gilt für cosh(x) und sinh(x):

cosh( ) ; sinh( )2 2

x x x xe e e ex x

− −+ −= =

Durch Addition dieser beiden hyperbolischen Funktionen erhält man die Exponentialfunktion:

cosh( ) sinh( )2 2

x x x xx xe e e e e e x x

− −+ −+ = ⇒ = + (2.48)

Die Berechnung dieser hyperbolischen Funktionen mittels CORDIC-Algorithmus wurde in Abschnitt 2.3.1 behandelt. Als maximal zulässiges Argument für φ gilt nach Gl. (2.36):

1.118173ϕ = ± Hierfür wollen wir nachfolgend ein beliebiges Argument x ∈ℝ einsetzen, das sich als Summe eines zulässigen Argumentes z mit |z|

15

* **

sinh( ) sinh( ln 2) sinh( ) cosh( ln 2) cosh( ) sinh( ln 2)x z m z m z m= + ⋅ = ⋅ ⋅ + ⋅ ⋅��������� ���������

(2.51a)

( )( )

1 1 1

1 1 1

sinh( ) 2 2 2 2 sinh( ) 2 sinh( )

cosh( ) 2 2 2 2 cosh( ) 2

*

cosh** ( )

m m m m

m m m m

z z z

z z z

− − − − −

− − − − −

⋅ ⋅ + = ⋅ + ⋅

⋅ ⋅ − = ⋅ − ⋅

( )1 2 sinh( ln 2) 2 cosh( ) sinh( ) 2 cosh( ) sinh( )m mz m z z z z− − ⇒ + ⋅ = ⋅ + − ⋅ − (2.51b) Die Gln. (2.50b) und (2.51b) ermöglichen es, die hyperbolischen Funktionen cosh(x) und sinh(x) für ein beliebiges Argument x ∈ℝ zu berechnen. Nach Gl. (2.48) erhält man daher für exp(x):

[ ] [ ]1exp( ) exp( ln 2) 2 2 cosh( ) sinh( ) 2 cosh( ) sinh( )m mx z m z z z z−= + ⋅ = ⋅ ⋅ + = ⋅ + (2.52) Der im Anhang A9 angegebene C-Quellcode berechnet die Exponentialfunktion exp(arg) mit arg ∈ℝ nach Gl. (2.52). Zunächst wird die Variable arg mittels Cast-Operator (int) durch ln_2 dividiert, um den ganzzahligen Faktor m für die spätere Stellenverschiebung zu erhalten. Das Argument z, für welches die hyperbolischen Funktionen direkt bestimmt werden, geht aus Gl. (2.49) hervor: m = (int)(arg / ln_2); z = arg – m * ln_2; Bei einer Hardware-Realisierung ist z der echt-gebrochenrationale Divisionsrest von arg/ln_2. Für die Initialwerte arg=0.5, arg=5.7 und arg=-5.7 berechnet A9 nach jeweils n=32 Iterationen folgende Funktionswerte: exp_fkt = 1.6487212710 // für arg = 0.5 exp_fkt = 298.8674009805 // für arg = 5.7 exp_fkt = 0.0033459655 // für arg = -5.7 Vergleicht man die Ergebnisse für arg=0.5 und arg=5.7 mit denen des Zahlenbeispiels aus Abschnitt 1, dann erkennt man, dass das Konvergenzverhalten der Exponentialfunktion hier nicht vom Argument abhängig ist! 2.3.4 Erweiterung des zulässigen Argument-Bereiches für artanh(), ln() und sqrt() artanh(): Das betragsmäßig größte Argument für den Areatangens hyperbolicus innerhalb des Konvergenz-Bereiches beträgt nach Gl. (2.39) ≈0.81, die Funktion weist aber einen Definitionsbereich von -1 bis +1 auf. Um auch den unzulässigen Bereich abzufangen, führen wir eine Beschreibung mittels Mantisse M und Exponent E ein, wobei gelten soll:

arg 1 2 bzw. 2 1 argE EM M− −= − ⋅ ⋅ = − (2.53)

mit: 0.5 1.0 ; 0M E≤ < ≥ Die Normalisierung der Mantisse entspricht nicht der Norm ANSI/IEEE 754, sondern ist der Original-Veröffentlichung von J. S. Walther entnommen [Walt71]. Da der ganzzahlige Exponent E nicht

16

Gesucht ist somit folgende Lösung:

( )artanh 1 2 ?EM −− ⋅ = Diese Funktion wird nun so umgeformt, dass ein Term artanh() entsteht, dessen Argument sich im Konvergenzbereich befindet. Nach [Papu01] gilt für artanh(x):

121

artanh( ) ln1

xx

x

+ = ⋅ −

( ) 1 12 21 1 2 2 2 2 2

artanh 1 2 ln ln 1 1 2 2 2 2

E E EE

E E E

M MM

M M

− −−

− −

+ − ⋅ − ⋅ ⋅⇒ − ⋅ = ⋅ = ⋅

− + ⋅ ⋅ ⋅

( )( )

1 12 2

4 2 2 4 2 2 2 artanh 1 2 ln ln

2 2

E E EE M MM

M M

−− ⋅ − − ⋅ ⋅⇒ − ⋅ = ⋅ = ⋅

( ) ( )

( )( )

12

1 1 2 2

2 2 2 2 2artanh 1 2 ln

2 2 2 2 ln l

2 2 2 2

n 22 2 2 2

E E EE

E EE

E E

E E

M MM

M M

M M M M

M M

M

M

M

M M

M

− −−

− −

− −

− −

+ − ⋅ − ⋅ ⋅− ⋅ = ⋅ +

+ − ⋅ + − − ⋅= ⋅ + ⋅ + − ⋅ − −

+ −+ − + ⋅

⋅

−

− ⋅

( ) 12 2

2 21

2 2artanh 1 2 ln ln 22 2

12 2

E

EE EE

E

M M

M MMM M

M M

−

−−−

−

− − ⋅+ + − ⋅ − ⋅ = ⋅ + ⋅

− − ⋅− + − ⋅

( ) 2 artanh 1 2 artanh( ) ln 2E EM T−⇒ − ⋅ = + ⋅ (2.55)

mit: 2 2

; 0.5 1.0 ; 02 2

E

E

M MT M E

M M

−

−− − ⋅= ≤ < ≥+ − ⋅

Um den Wertebereich von T zu ermitteln, betrachten wir folgende Zahlenkombinationen von M und E: M → 1; E = 0 => T = 0 M = 0.5; E = 0 => T = 0.5 M → 1; E = 1 => T = 0.2 M = 0.5; E = 1 => T = 0.555556 M → 1; E = 100 => T = 0.3 M = 0.5; E = 100 => T = 0.6 0 0.6T⇒ < ≤ (2.56) Der Wertebereich von T liegt somit im Konvergenzbereich der artanh-Funktion nach Gl. (2.39).

Beispiel 2.3: Es ist der Areatangens hyperbolicus für die folgenden Argumente mittels Gl. (2.55) zu bestimmen:

a) artanh(0.99609375) = 3.11818480 (Taschenrechner) b) artanh(0.00390625) = 0.00390627 “ c) artanh(0.50390625) = 0.55452812 “

a)

7

72

arg 1 2 0.99609375 0.11111111

2 0.00390625 0.00000001 0.1 2

0.5 ; 7 0.59937402

artanh(0.59937402) ln 2 3.11818479

E

E

M b

M b b

M E T

−

− −

= − ⋅ = =

⋅ = = = ⋅= = ⇒ =

⇒ + ⋅ =

17

b)

02

arg 1 2 0.00390625 0.00000001

2 0.99609375 0.11111111

0.99609375 ; 0 0.00390625

artanh(0.00390625) ln 2 0.00390627

E

E

M b

M b

M E T

−

−

= − ⋅ = =

⋅ = == = ⇒ =

⇒ + ⋅ =

c)

1

12

arg 1 2 0.50390625 0.10000001

2 0.49609375 0.01111111 0.11111110 2

0.9921875 ; 1 0.20500782

artanh(0.20500782) ln 2 0.55452812

E

E

M b

M b b

M E T

−

− −

= − ⋅ = =

⋅ = = = ⋅= = ⇒ =

⇒ + ⋅ =

N

ln(): Die Logarithmus-Funktion lässt sich nach Gl. (2.40) berechnen:

1

ln 2 artanh1

xx

x

− = ⋅ +

Für das Argument der artanh-Funktion gilt wieder der Konvergenzbereich nach Gl. (2.39), d.h. für den Konvergenzbereich des ln ergibt sich zunächst:

maxmax max maxmax

1 arg10.806932 arg ; 1 arg arg

1 1 arg

xx x x

x

±− = ± = ± − = ± ⋅ ± ⇒ =+ ∓

max maxmin maxmax max

1 arg 1 arg0.11 ; 9.36 0.11 9.36

1 arg 1 argx x x

− += ≈ = ≈ ⇒ ≤ ≤

+ − (2.57)

Der Numerus lässt sich auch hier durch Mantisse und Exponent ausdrücken, wobei im Hinblick auf eine spätere Umsetzung in VHDL eine von [Walt71] abweichende, jedoch ANSI/IEEE-konforme Darstellung gewählt wird. Die Mantisse liegt im Konvergenzbereich von Gl. (2.57):

( )ln 2 ln( ) ln 2EM M E⋅ = + ⋅ (2.58)

mit: 1.0 2.0M≤ < Beispiel 2.4: Es ist der natürliche Logarithmus für folgende Werte mittels Gl. (2.58) zu bestimmen:

a) ln(0.00390625) = -5.54517744 (Taschenrechner) b) ln(11.8750000) = 2.47443535 “

a)

82 0.00390625 0.00000001 1.0 2

1.0 ; 8 ln(1.0) 8 ln 2 5.54517744

EM b b

M E

−⋅ = = = ⋅= = − ⇒ − ⋅ = −

b)

32 11.875 1011.111 1.011111 2

1.484375 ; 3 ln(1.484375) 3 ln 2 2.47443535

EM b b

M E

⋅ = = = ⋅= = ⇒ + ⋅ =

N

sqrt(): Zunächst bestimmen wir den Konvergenzbereich der Wurzelfunktion nach Gl. (2.41):

1

0 41 10 0 max4 4 1

0 4

; 0.806932 argry

x r y rx r

−= + = − ⇒ = = ± = ±

+

18

( ) ( )1 1 1max max max max4 4 4arg arg 1 arg 1 argr r r− = ± ⋅ ± ⋅ ⇒ ⋅ = ± ⋅∓

max14max

1 arg 0.03 2.34

1 argr r

±= ⋅ ⇒ < <

∓ (2.59)

Der Radikand wird auch hier ANSI/IEEE-konform durch Mantisse und Exponent beschrieben, wobei nach [Walt71] zwischen geraden und ungeraden Exponenten unterschieden werden muss. In beiden Fällen wird der Konvergenzbereich von Gl. (2.59) eingehalten (für E mod 2 = 1 wird die Erhöhung des Exponenten durch Division der Mantisse ausgeglichen):

( )1

122 22 2 2E E Er M M M= ⋅ = ⋅ = ⋅

2 2 für: 1.0 2.0

mod 2 0 (gerade)

Er M M

E

⇒ = ⋅ ≤ <=

(2.60a)

( 1) 2 2 2 für: 0.5 2 1.0

mod 2 1 (ungerade)

Er M M

E

+⇒ = ⋅ ≤ <=

(2.60b)

2.3.5 Zusammenfassung CORDIC-Algorithmus Nachfolgend ist der CORDIC-Algorithmus für die beiden Betriebsarten Rotating und Vectoring sowie für die drei Modi linear (m=0), zirkular (m=1) und hyperbolisch (m=-1) zusammengefasst.

Mode m Rotating

0 0 01 , 0 , 0 , sgn( )n i i

x y zz z

ϕσ

= = =→ =

Vectoring

0 0 , 0sgn( )

n

i i

z yyσ

= →= −

0

01

(0, ) {0,1,..., 1}

mKS i n

==

= −

1 0(0, )

1(0, )

1

mit: 0

2

0 2

i iS i

i i i iS i

i i i

x x y

y y x

z z

σσ

+−

+−

+

= == + ⋅ ⋅

← = − ⋅

⇒ 00 0

n

n

x xy x z

== ⋅

1(0, )

1(0, )

1

0 2

2

i iS i

i i i iS i

i i i

x x

y y x

z z

σσ

+−

+−

+

=← = + ⋅ ⋅

= − ⋅

⇒ 00 0

n

n

x xz y x

==

1

10.607252935

(1, ) {0,1,..., 1}

mKS i n

==

= −

( )

(1, )1

(1, )1

(1, )1

2

2

0 arctan 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

−+

−+

−+

= − ⋅ ⋅= + ⋅ ⋅

← = − ⋅

⇒1

11

1

cos( )

sin( )n

n

x K

y K

ϕϕ

−

−= ⋅= ⋅

( )

(1, )1

(1, )1

(1, )1

2

0 2

arctan 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

−+

−+

−+

= − ⋅ ⋅← = + ⋅ ⋅

= − ⋅

⇒1 2 2

1 0 0

0 0arctan( )n

n

x K x y

z y x

−= ⋅ +=

1

11.207497068

( 1, ) {1,2,3,4,4,5,..,13,13,14,..., }

mKS i n

−

= −=

− =

( )

( 1, )1

( 1, )1

( 1, )1

2

2

0 artanh 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

− −+

− −+

− −+

= + ⋅ ⋅= + ⋅ ⋅

← = − ⋅

⇒1111

cosh( )

sinh( )n

n

x K

y K

ϕϕ

−−−−

= ⋅= ⋅

( )

( 1, )1

( 1, )1

( 1, )1

2

0 2

artanh 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

− −+

− −+

− −+

= + ⋅ ⋅← = + ⋅ ⋅

= − ⋅

⇒1 2 21 0 0

0 0artanh( )n

n

x K x y

z y x

−−= ⋅ −

=

Tabelle 2.4: Zusammenfassung CORDIC-Algorithmus

19

3 Prinzipielle Möglichkeiten der Implementierung Die Implementierung des CORDIC-Algorithmus ist immer unter dem Gesichtspunkt bestimmter Entwurfsziele zu sehen. Diese sind beispielsweise:

• Geschwindigkeit • Hardwareaufwand (notwendige Chipfläche) • Genauigkeit • Verlustleistung

Von den zahlreichen Umsetzungsmöglichkeiten sollen 4 Verfahren erläutert werden. Aspekte wie die Umwandlung in ein anderes Datenformat (Gleit- in Festkommaformat bzw. umgekehrt) sowie die Skalierung bleiben in diesem Abschnitt unberücksichtigt. 3.1 Bitparallele iterative CORDIC-Struktur In der bitparallelen iterativen CORDIC-Struktur nach Bild 3.1 wird jede der drei Differenzengleichungen für x, y und z in Hardware realisiert. Das bisher σi genannte Vorzeichen wird dort di genannt und ist abhängig von der Betriebsart, dem Modus sowie dem Vorzeichen des y- oder z-Registers. Zu Beginn der Iteration werden die Initialwerte über Multiplexer in die x-, y- und z-Register geladen.

In jedem Iterationszyklus durchlaufen die xi und yi die Barrelshifter und Addierer/Subtrahierer und das Ergebnis wird im zugehörigen Register gespeichert. Die gewünschte Anzahl an Shiftoperationen wird, abhängig von dem Modus, in jedem Schritt angepasst. Um den durch z repräsentierten Gesamtwinkel bestimmen zu können, sind die erforderlichen Teilwinkel (arctan bzw. artanh) in jeweils einem ROM-Speicher abgelegt, dessen Adresse in jedem Iterationszyklus inkrementiert wird (in nebenstehendem Bild ist nur ein ROM-Speicher skizziert). In der letzten Iteration werden die Ergebnisse für xn, yn und zn direkt von den Addierern/Subtrahierern gelesen. Zur Ablaufsteuerung und ROM-Adressierung muss das in Bild 3.1 gezeigte Rechenwerk noch um einen Zustandsautomaten (Steuerwerk) erweitert werden. Die hier beschriebene CORDIC-Struktur führt die Shift- und arithmetischen Operationen mit Datenwortbreite innerhalb eines Taktes aus. Für die beiden Barrelshifter bedeutet diese variabel durchführbare parallele Schiebeoperation aller Bits, dass eine synthetisierte Hardware auf einem FPGA nach [Andr98] zu einer hohen Anzahl an hintereinander geschalteten Logikzellen und dementsprechend zu relativ hohen Signallaufzeiten führt.

Bild 3.1: Bitparallele iterative CORDIC- Struktur [Andr98] 3.2 Bitserielle iterative CORDIC-Struktur Der Hardwareaufwand für die Umsetzung des CORDIC-Verfahrens kann, ausgehend von der bitparallelen CORDIC-Struktur, durch Verwendung bitserieller Arithmetik weiter reduziert werden. Die einfachere Logik und Verdrahtung eines solchen bitseriellen Designs erlaubt eine wesentlich höhere Taktrate als ein vergleichbares bitparalleles Design, allerdings werden jetzt auch pro Iteration w Taktzyklen benötigt (w entspricht der Datenwortbreite). Die bitserielle iterative CORDIC-Struktur nach Bild 3.2 setzt sich zusammen aus drei bitseriellen Addierern/Subtrahierern, drei Shiftregistern und

20

einem (bzw. zwei) seriellen ROM-Speicher und benötigt somit ein Minimum an Hardware. Die Länge der Shiftregisters ist identisch mit der Wortbreite. Ihre Initialisierung kann entweder bitseriell über die linksseitig skizzierten Multiplexer in w Takten oder bitparallel in einem Takt erfolgen. Einmal initialisiert, werden die Daten nach rechts geshiftet, durchlaufen den seriellen Addierer/Subtrahierer und werden linksseitig des zugehörigen Shiftregisters wieder eingeschoben.

Bild 3.2: Bitserielle iterative CORDIC-Struktur [Andr98] Jede Iteration erfordert w Taktzyklen, um das Zwischenergebnis in das Register zurückzuschreiben. Die Funktionalität der Barrelshifter wird zum einen ersetzt durch den Rechtsshift der Register x und y in jedem Takt, zum anderen über jeweils einen Multiplexer, der in jeder Iteration entsprechend der Anzahl der Schiebeoperationen das höchstwertige Bit selektiert. Ist beispielsweise i=3, dann ist durch die beiden Multiplexer bei der Berechnung von x3 und y3 das 4. Bit von rechts das höchstwertige Bit der gekreuzten Zweige in Bild 3.2 und geht in die Addition/Subtraktion ein. Die niederwertigsten drei Bits bleiben in diesem Zweig unberücksichtigt, für die letzten drei Schiebeoperationen muss jeweils eine ’0’ verwendet werden. Diese Tatsache sowie die komplette Ablaufsteuerung werden wieder von einem (nicht skizzierten) Automaten übernommen. Während der n-ten Iteration lassen sich die Ergebnisse an den Ausgängen der seriellen Addierer/Subtrahierer lesen, während die nächsten Initialisierungsdaten in die Register geshiftet werden. 3.3 Bitparallele CORDIC-Struktur, rein kombinatorisch Die beiden bisherigen Hardware-Realisierungen des CORDIC-Algorithmus sind iterativ, d.h. für die Berechnung eines Ergebnisses werden n Iterationen benötigt. Es ist jedoch auch möglich, die für jede Iterationsstufe erforderliche Hardware zu implementieren. Dieser Ansatz führt auf eine parallele, rein kombinatorische CORDIC-Struktur, die nach [Andr98] als “unrolled CORDIC processor“ bezeichnet wird und in Bild 3.3 dargestellt ist. Zwei wesentliche Punkte sind hier zu nennen: Zum einen führt jeder Shifter eine konstante Shiftoperation durch, d.h. er ist hard-wired. Zum anderen sind die Look-Up-Werte der Teilwinkel als Konstanten für jeden Addierer/Subtrahierer im z-Zweig aufgeführt. Da die Logik rein kombinatorisch arbeitet, sind keine Register erforderlich. Die Verzögerungszeit bis zum Anliegen der gültigen Ergebnisse am Ausgang ist entsprechend hoch, jedoch fällt hier im Vergleich zu Bild 3.1 die Einhaltung der Setup- und Hold-Zeit der Register zur Zwischenspeicherung weg.

21

Bild 3.3: Bitparallele CORDIC-Struktur, rein kombinatorisch [Andr98] 3.4 Bitparallele CORDIC-Struktur mit Pipeline-Registern Beim Einsatz von FPGAs als Zielhardware macht es keinen Sinn, eine umfangreiche kombinatorische Logik nach Bild 3.3 zu synthetisieren. Besser ist es, zwischen den Addierern/Subtrahierern Pipeline-Register vorzusehen. Da in einem FPGA jede Logikzelle ohnehin ausgangsseitig über ein Flipflop verfügt, entsteht durch zusätzliche Pipeline-Register kein Mehraufwand an Hardware. Bei einem konstanten Strom an Eingangsdaten weist diese Architektur die höchste Durchsatzrate auf. Nach einer gewissen Latenz, die abhängig ist von der Stufenzahl der Pipeline und damit der Genauigkeit, liegt in jedem Takt ein Ergebnis am Ausgang vor! Ein Beispiel für eine solche Pipeline-Architektur ist der in [Timm93] vorgestellte und für Anwendungen in der digitalen Signalverarbeitung entwickelte Chip IMSCOR24 (die im nächsten Kapitel beschriebene VHDL-Implementierung orientiert sich teilweise hieran). Extern liegen die Daten im 32-Bit-Gleitkommaformat vor, intern arbeitet der Chip mit einem 32-Bit-Festkommaformat und einer 44-stufigen Pipeline (Iterationen einschließlich Skalierung und Datenformatumwandlung). Mit diesem Chip wurden erstmals alle möglichen CORDIC-Berechnungen durchgeführt.

22

XFP[31:0] YFP[31:0] ZFP[31:0]

VorstufeX_INT[31:0] Y_INT[31:0] Z_INT[31:0]

24 CORDIC-Iterationen

QX_INT[31:0] QY_INT[31:0] QZ_INT[31:0]

SkalierungQX_MULT[32:0] QY_MULT[32:0]

Endstufe

XFP_OUT[31:0] YZFP_OUT[31:0]

Operationen fürln-Berechnung

XFP[31:0] YFP[31:0] ZFP[31:0]

VorstufeX_INT[31:0] Y_INT[31:0] Z_INT[31:0]

24 CORDIC-Iterationen

QX_INT[31:0] QY_INT[31:0] QZ_INT[31:0]

SkalierungQX_MULT[32:0] QY_MULT[32:0]

Endstufe

XFP_OUT[31:0] YZFP_OUT[31:0]

Operationen fürln-Berechnung

4 Implementierung in VHDL Der in diesem Abschnitt beschriebene strukturelle VHDL-Entwurf berücksichtigt nicht den vollständigen CORDIC-Algorithmus nach Tabelle 2.4, sondern lediglich die Berechnung der trigonometrischen Funktionen cos und sin sowie Exponential-, Logarithmus- und Wurzelfunktion. Alle implementierten Funktionen sowie der zugehörige Wertebereich der Eingabedaten sind in Tabelle 4.1 zusammengefasst. OP_CODE[7:0]

Eingabe

XFP YFP ZFP Ergebnis

XFP_OUT YZFP_OUT Wertebereich der

Eingabedaten X“01“ 1 0 z cos(z) sin(z) 0 2z π≤ ≤ X“03“ 1 0 z exp(z) 0 16z < X“04“ r r 0 0 ln(r) 1.0842022E-19 1.8446743E19r≤ ≤ X“05“ r r 0 sqrt(r) 0 9.536743E-7 2.0971519E6r≤ ≤

Tabelle 4.1: Implementierte Funktionen des CORDIC-Verfahrens 4.1 Funktionaler Überblick Bild 4.1 zeigt das Blockschaltbild für die Hardware-Realisierung. Die Eingangsdaten liegen extern im 32-Bit-Gleitkommaformat vor (1 Vorzeichenbit, 8-Bit-Charakteristik und 23-Bit-Mantisse) und werden mit XFP, YFP und ZFP bezeichnet. In der Vorstufe werden diese Eingangsdaten in das interne 32-Bit-Festkommaformat umgewandelt, wobei zur möglichst einfachen Umsetzung der Addition/Subtraktion die Zweikomplement-Darstellung gewählt wird. Für XFP und YFP wird ein gemeinsamer 10-Bit-Referenzexponent QREF_EXP_XY gespeichert. Das in Bild 4.2 gezeigte interne Festkommaformat berücksichtigt neben der um das Hidden-Bit auf 24 Bit erweiterten Mantisse noch 2 Overflow- und 5 Addition-Bits.

Bild 4.1: Blockdiagramm der Hardware-Umsetzung Die eigentliche Berechnung des Cordic-Algorithmus erfolgt in 24 Iterationen, in denen entweder die einfache Shiftfolge nach Gl. (2.31a) für den zirkularen Modus oder die komplexere Shiftfolge nach Gl. (2.34) für den hyperbolischen Modus zugrunde gelegt wird. Die Funktionswerte für arctan(2-i) bzw. artanh(2-i) sind intern in einer Tabelle (ROM-Speicher) mit 28 Nachkommastellen abgelegt.

23

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

ManX[23:0]

24-Bit-Barrelshifter

ManXshft[23:0] OV_BITS_X[1:0]

ManY[23:0]

24-Bit-Barrelshifter

ManYshft[23:0] OV_BITS_Y[1:0]

CONV_FP_TC

XFP[31] YFP[31]

X_INT[31:0]

CONV_FP_TCY_INT[31:0]

QX_INT[31:0] QY_INT[31:0]

XFP[30:23] YFP[30:23]

Referenzexponent

QREF_EXP_XY[9:0]

CountY[7:0] CountX[7:0]

ManZ[23:0]

24-Bit-Barrelshifter

ManZshft[23:0] OV_BITS_Z[1:0]

CONV_FP_TCZ_INT[31:0]

QZ_INT[31:0]

ZFP[30:23]

Argument-Erweiterung(exp-Funktion)

QUOTIENT[4:0]REST[23:0]

Fallunterscheidung(cos/sin-Funktion):

z pi/2 and z pi and z 3pi/2 and z

24

Für ZFP sind die Verhältnisse etwas komplizierter, weil die Funktionswerte für arctan()bzw. artanh() intern in einer Tabelle im 32-Bit-Festkommaformat nach Bild 4.2 abgelegt sind, d.h. es handelt sich um denormalisierte Mantissen und der Exponent (ohne Bias) ist 0. Ist der Exponent von ZFP, d.h. ZFP(30 downto 23), größer als 127d, dann erfolgt ein Linksshift um die Differenzstellenzahl “ZFP(30 downto 23) – 127“, andernfalls ein Rechtsshift um die Differenzstellenzahl “127 – ZFP(30 downto 23)“. Um zwischen Links- und Rechtsshift unterscheiden zu können, besitzt der 24-Bit-Barrelshifter das Bit SHIFT_DIRECTION. Bei SHIFT_DIRECTION = ’0’ erfolgt ein Rechtsshift, bei ’1’ entsprechend ein Linksshift um die über COUNT angegebene Stellenzahl. Die geshiftete Mantisse wird mit ManZshft[23:0] bezeichnet. Ihr zugehöriger Referenzexponent ist immer 0, so dass dieser nicht explizit abgespeichert werden muss. Die 24-Bit-Barrelshifter liefern zusätzlich noch die notwendigen Overflow-Bits OV_BITS_X[1:0] OV_BITS_Y[1:0] und OV_BITS_Z[1:0]. Bei einem Rechtsshift sind diese OV-Bits “00“ (für XFP und YFP immer der Fall), bei einem Linksshift um eine bzw. zwei Stellen werden die geshifteten Bits in die OV-Bits übernommen. Argument-Erweiterung der exp-Funktion: Im hyperbolischen Modus erfolgt zur Erweiterung des Argument-Bereiches bei der Berechnung der Exponentialfunktion eine Division von ZFP durch ln(2), bei der ein ganzzahliger Quotient, QUOTIENT[4:0], und ein echt-gebrochenrationaler Rest, REST[23:0], entsteht. Eine ausführliche Erläuterung dieser Zusammenhänge erfolgt im nächsten Abschnitt. Die eigentliche Umwandlung der geshifteten Mantissen ManXshft[23:0], ManYshft[23:0] und ManZshft[23:0] in das interne 32-Bit-Festkommaformat in Zweikomplement-Darstellung erfolgt in den Blöcken CONV_FP_TC, siehe Anhang A10. Ist beispielsweise die Gleitkommazahl XFP positiv, d.h. XFP(31) = ’0’, dann werden XFP(31), OV_BITS_X[1:0] sowie ManXshft[23:0] direkt in das Festkommaformat X_INT bzw. QX_INT übernommen (der Index Q deutet an, dass es sich um ein Flipflop bzw. Register handelt). Die noch fehlenden Addition-Bits werden alle auf 0 gesetzt. Ist die Gleitkommazahl hingegen negativ, dann wird nur XFP(31) direkt übernommen. Alle Stellen der OV-Bits und der geshifteten Mantisse werden invertiert und die Addition-Bits alle auf 1 gesetzt. Um das Komplement zu erhalten, wird abschließend noch eine 1 addiert. Fallunterscheidung für cos und sin: Der Argument-Bereich für die trigonometrischen Funktionen cos und sin soll zwischen 0 und 2π liegen. Damit das Verfahren konvergiert, darf φ nach Gl. (2.14) den Wert 1.74 nicht überschreiten. Somit muss bei Auswahl des OP-Codes X“01“ das Argument (arg=Z_INT) im Automatenzustand Z1 gemäß nachfolgender Tabelle modifiziert sowie der 2-Bit-Signalvektor FALLUNTERSCHEIDUNG entsprechend gesetzt werden:

Argument FALLUNTER-SCHEIDUNG

modifiziertes Argument

cos-Zweig (x)

sin-Zweig (y)

arg 2π≤ 00 QZ_INT Z_INT ≤ 01 QZ_INT Z_INT 2π ≤ 10 QZ_INT Z_INT π ≤ 11 QZ_INT Z_INT 3 2π

25

zwischen 107d und 147d. Ausgehend von der Dezimalzahl 1 kann das Komma damit um 20 Stellen nach links bzw. rechts verschoben werden. Ist der Exponent ungerade, dann wird die Mantisse durch 2 geteilt, andernfalls bleibt sie unverändert. Der Teilung durch 2 entspricht ein Rechtsshift, bei welchem nach dem Vorzeichenbit von links eine ’0’ eingeschoben wird. 4.3 Erweiterung des Argument-Bereiches der exp-Funktion Nach Gl. (2.48) setzt sich die Exponentialfunktion additiv aus den hyperbolischen Funktionen cosh und sinh zusammen, für deren maximal zulässiges Argument Gl. (2.36) gilt. Für ein beliebiges Argument x ∈ℝ erhält man nach Gl. (2.52):

exp( ) exp( ln 2) 2 [cosh( ) sinh( )]mx z m z z= + ⋅ = ⋅ + Um die Exponentialfunktion berechnen zu können, müssen also zunächst die Größen z und m aus x bestimmt werden. Dies führt auf eine Division mit einem ganzzahligen Quotienten (m) und einem echt-gebrochenrationalen Rest (z): ln 2z x m= − ⋅ (4.1) Um die Vorgehensweise zu verdeutlichen, gehen wir zunächst von einer reduzierten Stellenanzahl mit einer 10-Bit-Mantisse (einschließlich Hidden-Bit) aus und betrachten allgemein folgendes Divisionsbeispiel: x = 9.171875000d = 1001.001011 = 1.001001011 * 2^3 ln2 = 0.693359375d = 0.101100011 Der Dezimalpunkt des Dividenden wird nun um eine Stelle nach links verschoben, wodurch sich seine Stellenzahl um ein Bit erhöht: x = 0.1001001011 * 2^4 Nun erfolgt die Division mit Rückstellen des Zwischenrestes (siehe z.B. [Beul08]): Stufe 0: 01001001011 : 0101100011 = 0.1101 * 2^4 = 1101 = 13d -0101100011 0 11110000101 => 0

Stufe 1: 10010010110 -0101100011 0 00111010000 => 1

Stufe 2: 01110100000 -0101100011 0 00011011010 => 1

Stufe 3: 00110110100 -0101100011 0 11011101110 => 0

Stufe 4: 01101101000 -0101100011 0 00010100010 => 1

� x = 1101 = 13d ; z = 0.001010001 = 0.158203125d

Die hier beschriebene Umsetzung des Divisionsverfahrens über 5 (stellenverschobene) Subtraktionen erlaubt einen maximalen Exponenten (ohne Bias) von E=3. Damit wird das maximal mögliche Argument für die Exponentialfunktion beschränkt auf:

26

3max 1111.1111111111 1.1111111111111 2 16x = = ⋅ 3, so führt dies zu falschen Ergebnissen! Der im Anhang angegebene VHDL-Quellcode fängt diese Wertebereichsüberschreitung nicht ab. Ist E 1

Stufe 2: 10011000100 -0101100011 0 00111111110 => 1

Stufe 3: 01111111100 -0101100011 0 00100110110 => 1

� x = 111 = 7d ; z = 0.010011011 = 0.302734375d

27

4.4 CORDIC-Iterationen Nach der Vorstufe werden, ausgehend von dem 32-Bit-Festkommaformat nach Bild 4.2, 24 Iterationen ausgeführt. Da für dieses Format das Zweierkomplement gewählt wurde, sind die arithmetischen Operationen “Addition“ und “Subtraktion“ besonders einfach. Die Vorgehensweise soll der Vollständigkeit halber für beide Operationen anhand von typischen Zahlenbeispielen demonstriert werden, auch wenn die Abfrage der hierbei auftretenden Bereichsüberschreitungen im VHDL-Code nicht implementiert ist. Addition: Nachfolgend aufgeführte Beispiele für die Addition von 4-stelligen Zweikomplementzahlen sind [Hoff93] entnommen. C OV C OV C OV a1) +2 0010 a2) +7 0111 b1) -1 1111 +2 0010 +7 0111 -1 1111 +4 00100 0 0 +14 01110 0 1 -2 11110 1 0

b2) -7 1001 c1) +2 0010 c2) +1 0001 -7 1001 -1 1111 -7 1001 -14 10010 1 1 +1 10001 1 0 -6 01010 0 0 Tabelle 4.4 zeigt die sich aus den Zahlenbeispielen ergebenden Fälle der Addition im Zweierkomplement.

Fall An-1 Bn-1 Sn Sn-1 Summe a1 0 0 0 0 S[n] = A + B a2 0 0 0 1 S[n+1] = A + B “Überlauf“ b1 1 1 1 1 S[n] = A + B – 2^n b2 1 1 1 0 S[n+1] = A + B “Überlauf“ c1 0 1 1 0 S[n] = A + B – 2^n c2 0 1 0 1 S[n] = A + B

Tabelle 4.4: Die verschiedenen Fälle der Addition im Zweikomplement In den Fällen a2 und b2 findet eine Bereichsüberschreitung (Überlauf, Overflow) statt, da das Ergebnis nicht mehr mit n=4 Stellen darstellbar ist. Die hinzukommende Summenstelle Sn kann dann als Vorzeichenstelle interpretiert werden. Die Stelle Sn entspricht dem Übertrag (carry out). Der Übertrag, der in den Fällen b1 und c1 auftritt, gehört nicht zur Darstellung des n -stelligen Ergebnisses und wird deshalb abgeschnitten (A+B-2n). Alle geschilderten Fälle können nach [Hoff93] in zwei Schritten zusammengefasst werden:

1. S[n]=A+B bilden und den Übertrag Sn merken. Die Überlaufbedingung berechnen. 2. Wenn kein Überlauf aufgetreten ist (OV=0), dann ist die Summe gleich S[n], ansonsten

S[n+1]. Subtraktion: Auch die Subtraktion von Zweikomplementzahlen soll anhand von Zahlenbeispielen für n=4 veranschaulicht werden. C OV C OV C OV a1) (+7) 0111 a2) (+5) 0101 b1) (+2) 0010 -(+3) 0011 -(+7) 0111 -(-5) 1011 +4 00100 0 0 -2 11110 1 0 +7 10111 1 0

b2) (+3) 0011 c1) (-5) 1011 c2) (-5) 1011 -(-6) 1010 -(+3) 0011 -(+6) 0110 +9 11001 1 1 -8 01000 0 0 -11 00101 0 1

28

C OV C OV d1) (-6) 1010 d2) (-3) 1101 -(-4) 1100 -(-7) 1001 -2 11110 1 0 +4 00100 0 0 Die sich ergebenden Fälle sind wieder in Tabelle 4.5 zusammengefasst.

Fall An-1 Bn-1 Dn Dn-1 Summe a1 0 0 0 0 D[n] = A - B a2 0 0 1 1 D[n] = A – B – 2^n b1 0 1 1 0 D[n] = A - B – 2^n b2 0 1 1 1 “Überlauf“ → Wert zu groß c1 1 0 0 1 D[n] = A - B c2 1 0 0 0 “Überlauf“ → Wert zu klein d1 1 1 1 1 D[n] = A - B – 2^n d2 1 1 0 0 D[n] = A - B

Tabelle 4.5: Die verschiedenen Fälle der Subtraktion im Zweikomplement Hier findet in den Fällen b2 und c2 eine Bereichsüberschreitung aufgrund einer zu geringen Stellenzahl statt. Invertiert man die hinzukommende Differenzstelle Dn, dann kann dieses Bit als Vorzeichenstelle interpretiert werden. Der Übertrag in den Fällen a2, b1 und d1 kann wieder abgeschnitten werden, da er nicht zum Ergebnis gehört (A-B-2n). Auch die Subtraktion im Zweikomplement kann in zwei Schritten zusammengefasst werden:

1. D[n]=A-B bilden und den Übertrag Dn merken. Die Überlaufbedingung berechnen. 2. Wenn kein Überlauf aufgetreten ist (OV=0), dann ist die Differenz gleich D[n], ansonsten wird

das Bit Dn invertiert und die Differenz ist D[n+1]. Wie bereits zu Beginn dieses Abschnitts erwähnt, fängt der im Anhang A10 aufgeführte VHDL-Quellcode die hier angegebenen Bereichsüberschreitungen sowohl bei Addition als auch bei Subtraktion nicht ab. Bei allen eingegebenen Testdaten ist ein solcher Fall (fehlerhaftes Ergebnis) bisher nicht aufgetreten! CORDIC-Architektur: Als Architektur kommt das Prinzip der in Bild 3.1 gezeigten bitparallelen iterativen CORDIC-Struktur zur Anwendung. Es werden im Prozess A_BEST im Zustand Z3 24 CORDIC-Iterationen ausgeführt. Nach der 24. Iteration wechselt der Automat in den Zustand Z4. Als Laufindex für den aktuellen CORDIC-Iterationszyklus dient der 5-Bit-Vektor QS[4:0]. Der Index i der Shiftfolgen wird im Quellcode über die Variable N ausgedrückt, die in jeder Iteration die Anzahl der Rechtsverschiebungen für die beiden bitparallelen Barrelshifter im x- und y-Zweig angibt. Sie ist vom Typ integer, da mit ihrem Wert innerhalb einer for-loop-Anweisung eine Umsortierung der Signalvektoren und damit die eigentliche Shiftoperation erfolgt. Zwischen QS und N besteht in Abhängigkeit von m folgender Zusammenhang:

QS (dezimal) 0 1 2 3 4 5 6 7 8 9 10 11 N (m=1) 0 1 2 3 4 5 6 7 8 9 10 11 N (m=-1) 1 2 3 4 4 5 6 7 8 9 10 11

QS (dezimal) 12 13 14 15 16 17 18 19 20 21 22 23

N (m=1) 12 13 14 15 16 17 18 19 20 21 22 23 N (m=-1) 12 13 13 14 15 16 17 18 19 20 21 22

Tabelle 4.6: Zusammenhang zwischen QS und N in Abhängigkeit von m Nachfolgender Codeausschnitt zeigt die Bestimmung von N aus QS gemäß Tabelle 4.6. when Z3 => if m = 1 then -- einfache Shiftfolge für zirkularen Modus N:=Bit5Int(QS); else -- komplexere Shiftfolge für hyperbolischen Modus if QS < "00100" then

29

N:=Bit5Int(QS) + 1; elsif QS < "01110" then N:=Bit5Int(QS); else N:=Bit5Int(QS) - 1; end if; end if; Die Ergebnisse der jeweiligen Shiftoperationen für den x- bzw. y-Zweig stehen in den Variablen VX_SHIFT[31:0] und VY_SHIFT[31:0]. Um die CORDIC-Berechnungen durchzuführen, ist für jeden Datenpfad eine zusätzliche 32-Bit-Variable vorgesehen, namentlich VX_INT, VY_INT und VZ_INT. Auf diese Weise kann durch Zugriff auf die Bitstelle VX_INT(31) bzw. VY_INT(31) das Vorzeichen oi für die nachfolgende CORDIC-Iteration ermittelt werden. Die Teilwinkel arctan(2-i) und artanh(2-i) sind jeweils als ROM-Speicher (ROM_arctan bzw. ROM_artanh) mit der Organisation 24*32 in der Quellcodedatei PACK_CORDIC.vhd abgelegt. Jedes 32-Bit-Datenwort setzt sich zusammen aus einem Vorzeichenbit (immer positiv, d.h. ’0’) sowie 3 Vor- und 28 Nachkommastellen. Die Umwandlung der Dezimalzahlen in die entsprechenden Dualzahlen erfolgt sowohl hier als auch bei dem Konstantenmultiplizierer im nächsten Abschnitt über das Programm dez_bin.exe der Website http://delphi.zsg-rottenburg.de/download.html. Im hyperbolischen Modus startet die Shiftfolge mit dem Wert 1, die Adressierung des Speichers ROM_artanh beginnt aber ab der Adresse 0. Deshalb wird hier von der Variablen N der Wert 1 subtrahiert. Die Iterationsgleichungen in Abhängigkeit von m und oi können nun direkt mittels Tabelle 2.4 angegeben werden: when Z3 => if m = 1 then -- zirkularer Modus if oi = 1 then VX_INT := QX_INT - VY_SHIFT; VY_INT := QY_INT + VX_SHIFT; VZ_INT := QZ_INT - ROM_arctan(N); else VX_INT := QX_INT + VY_SHIFT; -- oi ist -1 VY_INT := QY_INT - VX_SHIFT; VZ_INT := QZ_INT + ROM_arctan(N); end if; else -- hyperbolischer Modus if oi = 1 then VX_INT := QX_INT + VY_SHIFT; VY_INT := QY_INT + VX_SHIFT; VZ_INT := QZ_INT - ROM_artanh(N-1); else VX_INT := QX_INT - VY_SHIFT; -- oi ist -1 VY_INT := QY_INT - VX_SHIFT; VZ_INT := QZ_INT + ROM_artanh(N-1); end if; end if;

Abschließend findet die Vorzeichenbestimmung oi aus VZ_INT(31) bzw. VY_INT(31) für die nächste Iteration statt und die Variablenwerte werden in die Register QX_INT, QY_INT und QZ_INT übernommen. Der Abbruch des CORDIC-Verfahrens nach 24 Iterationen und damit der Wechsel in den Zustand Z4 erfolgt im Prozess F_BEST: when Z3 => if QS = "10111" then FOLGE_Z

30

4.5 Skalierung – Konstantenmultiplizierer Eine dem CORDIC-Verfahren innewohnende Eigenschaft ist die Skalierung der x- und y-Werte in Abhängigkeit von der verwendeten Shiftfolge, siehe hierzu die Gln. (2.11) und (2.35). Da hier ein Faktor konstant ist, kann die Hardware-Realisierung der Multiplikation vereinfacht werden. Bevor eine solche Konstantenmultiplikation ausgeführt wird, erfolgt die Umwandlung des nicht konstanten Faktors (QX_INT bzw. QY_INT; im Folgenden Dualzahl LD genannt) von der Zweikomplement- in die Vorzeichen-Betrag-Darstellung. Da die Konstanten zur Skalierung positiv sind, ist das Produktvorzeichen identisch mit dem Vorzeichen für LD und kann am Ende der Multiplikation einfach übernommen werden. Als einführendes Beispiel betrachten wir die Multiplikation einer 12-Bit-Dualzahl V[11:0] mit einer 8-Bit-Konstanten K[7:0]. Das Ergebnis ist eine 20-Bit-Dualzahl P[19:0]:

( )( )( )( )

11 10 9 8 2 1 011 10 9 8 2 1 0

11 10 9 8 11 10 9 8

7 6 5 47 6 5 4

3 2 1 03 2 1 0

2 2 2 2 2 2 2

2 2 2 2

2 2 2 2

2 2 2 2

P K v v v v v v v

K v v v v

K v v v v

K v v v v

= ⋅ ⋅ + ⋅ + ⋅ + ⋅ + + ⋅ + ⋅ + ⋅

= ⋅ ⋅ + ⋅ + ⋅ + ⋅

+ ⋅ ⋅ + ⋅ + ⋅ + ⋅

+ ⋅ ⋅ + ⋅ + ⋅ + ⋅

…

Die folgende Struktur aus [Jork04] zeigt die Bildung der Teilprodukte und ihre Zusammenfassung zum Gesamtprodukt: k7 k6 k5 k4 k3 k2 k1 k0 * v0 k7 k6 k5 k4 k3 k2 k1 k0 * v1 k7 k6 k5 k4 k3 k2 k1 k0 * v2 k7 k6 k5 k4 k3 k2 k1 k0 * v3 --------------------------------- u11u10 u9 u8 u7 u6 u5 u4 u3 u2 u1 u0

k7 k6 k5 k4 k3 k2 k1 k0 * v4 k7 k6 k5 k4 k3 k2 k1 k0 * v5 k7 k6 k5 k4 k3 k2 k1 k0 * v6 k7 k6 k5 k4 k3 k2 k1 k0 * v7 --------------------------------- r11r10 r9 r8 r7 r6 r5 r4 r3 r2 r1 r0

k7 k6 k5 k4 k3 k2 k1 k0 * v8 k7 k6 k5 k4 k3 k2 k1 k0 * v9 k7 k6 k5 k4 k3 k2 k1 k0 * v10 k7 k6 k5 k4 k3 k2 k1 k0 * v11 --------------------------------- s11s10 s9 s8 s7 s6 s5 s4 s3 s2 s1 s0

Die Zwischensummen von je vier Partialprodukten werden zu einem Modul zusammengefasst. Jede Teilsumme kann in Abhängigkeit von den vier Bitstellen der Variablen 16 verschiedene Werte mit einer Wortlänge von 12 Bit annehmen. Ein solches Modul lässt sich als ROM-Speicherblock mit der Größe 16*12 Bit bilden. Dieses Prinzip soll nun zur Skalierung der Ergebnisse des CORDIC-Verfahrens angewandt werden. Zunächst können wir einen allgemein gültigen Zusammenhang zwischen den Längen der Konstanten und der Dualzahl, der Anzahl der ROM-Speicher und deren jeweiliger Kapazität (Organisation) sowie der Anzahl der Adressbits angeben: Länge der Konstanten: KL

Länge der Dualzahl: DL

Anzahl der ROM-Speicher: ROMN

Anzahl der Adressbits: ROMA

Anzahl der Datenbits: ROMD

Kapazität des Speichers: 2 ROMAROM ROMK D= ∗

31

31 025 24

0

0

0

0

0

025

25

25

25

25

25

∑

Vorkommastellen Nachkommastellen

32-Bit-ROM-Datenwort31 025 24

0

0

0

0

0

025

25

25

25

25

25

∑

Vorkommastellen Nachkommastellen

32-Bit-ROM-Datenwort

D ROM ROML N A⇒ = ⋅ (4.2)

( ) 2 2ROM ROMA AROM ROM K ROMK D L A= ∗ = ∗ + (4.3) Die Dualzahl LD weist nach Bild 4.2 zunächst 31 Stellen (ohne Vorzeichenbit) auf. Um eine möglichst sinnvolle Einteilung zu erhalten, wird LD auf 30 Stellen reduziert, d.h. das niederwertigste Addition-Bit fällt weg. Mit 5 ROM-Speichern ergeben sich auf diese Weise nach Gl. (4.2) 6 Adressbits: 30 5 6D ROM ROML N A= = ⋅ = ⋅ Weiterhin wird als Datenwortbreite DROM eine Zweierpotenz angestrebt, was auf eine Konstantenlänge von LK=26 Bit führt, siehe Bild 4.4. Für die Darstellung eines Skalierungsfaktors von ≈1.2 ist eine Vorkommastelle ausreichend. Ein ROM-Speicher muss somit nach Gl. (4.3) folgende Kapazität aufweisen:

( ) 6 2 2 (26 6) 64 32ROMAROM K ROMK L A= ∗ + = ∗ + = ∗

Bild 4.4: Datenwortbreite des ROM-Speichers In Tabelle 4.7 ist ein kleiner Ausschnitt eines ROM-Speichers für die Skalierung mit K≈0.607252 (grau hinterlegte Felder) aufgeführt. Die vollständigen Tabellen für beide Skalierungsfaktoren (zirkular und hyperbolisch) sind im Package-File PACK_CORDIC.vhd im Anhang A10 zu finden. Bild 4.5 zeigt schließlich den sich ergebenden Konstantenmultiplizierer, bestehend aus 5 ROM-Speichern und 4 Addierern. Das Produkt (ohne Vorzeichen) umfasst LK+LD=56 Bits, wobei die höchstwertigen 4 Bits die Vorkommastellen sind. Die 6 niederwertigsten Bits des 1. ROMs sind bereits die 6 niederwertigsten Ergebnisbits. Jede der ersten drei Addierstufen liefert weitere 6 Bits zum Endergebnis. Die Summe der letzten Addierstufe stellt die 32 höchstwertigen Bits des Produktes dar. Eingang SUM SUM (dezimal) SUM (binär) 000000 0K 00.000000000000000 0000000.0000000000000000000000000 000001 1K 00.607252925634384 0000000.1001101101110100111011011 000010 2K 01.214505851268768 0000001.0011011011101001110110110 000011 3K 01.821758776903152 0000001.1101001001011110110010001 000100 4K 02.429011702537536 0000010.0110110111010011101101100 000101 5K 03.036264628171920 0000011.0000100101001000101000111 000110 6K 03.643517553806304 0000011.1010010010111101100100010 000111 7K 04.250770479440688 0000100.0100000000110010011111101

… … … … … … … …

111111 63K 38.256934314966192 0100110.0100000111000110011100101

Tabelle 4.7: Ausschnitt eines ROM-Speichers für 0.607252K ≈

32

ROM_1[31:0]ROM 64*32E[6:1]

ROM 64*32E[12:7]

ROM 64*32E[18:13]

ROM 64*32E[24:19]

ROM 64*32E[30:25]

ROM_2[31:0]

ROM_3[31:0]

ROM_4[31:0]

ROM_5[31:0]

ROM_1[31:6]

ADD_1[31:6]

“000000“

“000000“

ADD_2[31:6]

“000000“

“000000“

ADD_3[31:6]

ADD

ADD

ADD

ADD

E[31:0]

ADD_1[5:0]

ROM_1[5:0]

ADD_2[5:0]

ADD_4[31:0]

ADD_3[5:0]

A[55:0]

ROM_1[31:0]ROM 64*32E[6:1]

ROM 64*32E[12:7]

ROM 64*32E[18:13]

ROM 64*32E[24:19]

ROM 64*32E[30:25]

ROM_2[31:0]

ROM_3[31:0]

ROM_4[31:0]

ROM_5[31:0]

ROM_1[31:6]

ADD_1[31:6]

“000000“

“000000“

ADD_2[31:6]

“000000“

“000000“

ADD_3[31:6]

ADD

ADD

ADD

ADD

E[31:0]

ADD_1[5:0]

ROM_1[5:0]

ADD_2[5:0]

ADD_4[31:0]

ADD_3[5:0]

A[55:0]

ROM 64*32E[6:1]

ROM 64*32

A[31:0]E[12:7]

E[18:13]

E[24:19]

E[30:25]

MUX

“000000“

ADDADD[31:0]

MUX

REGREG[31:0]

REG[31:6]E[31:0]

X[31:6]

Y[31:0]

ROM 64*32E[6:1]

ROM 64*32

A[31:0]E[12:7]

E[18:13]

E[24:19]

E[30:25]

MUX

“000000“

ADDADD[31:0]

MUX

REGREG[31:0]

REG[31:6]E[31:0]

X[31:6]

Y[31:0]

Bild 4.5: Kombinatorische Realisierung der Konstantenmultiplikation

Der Aufwand an ROM-Speichern und Addierern lässt sich reduzieren, wenn die Operationen nicht kombinatorisch, sondern sequentiell ausgeführt werden. Der Konstantenmultiplizierer nach Bild 4.6 kommt neben zwei zusätzlichen Multiplexern mit nur 2 ROM-Speichern und einem Addierer aus. Weiterhin ist noch eine nicht skizzierte Ablaufsteuerung erforderlich. Diese Lösung erfordert 4 Taktzyklen für eine Multiplikation und ist im VHDL-Quellcode im Anhang A10 implementiert.

Bild 4.6: Sequentielle Realisierung der Konstantenmultiplikation Zur Simulation mit ModelSim XE werden die ROM-Speicher im Package-File als Array deklariert: type MMX64_32 is array(63 downto 0) of std_logic_vector(31 downto 0);

Im zirkularen Modus müssen für die Kosinus- und Sinusfunktion sowohl der x- als auch der y-Wert skaliert werden, d.h. in 4 Taktzyklen erfolgen parallel zwei Konstantenmultiplikationen. Somit sind für den Skalierungsfaktor nach Gl. (2.11) insgesamt vier ROM-Speicher notwendig: constant K1_ROM_X1,K1_ROM_X2,K1_ROM_Y1,K1_ROM_Y2: MMX64_32:= -- zirkularer Modus ("01001100100000111000110011100101", -- 111111 "01001011010011001010001100001010", -- 111110 ... ... "00000000000000000000000000000000"); -- 000000

33

Im hyperbolischen Modus muss für die Berechnung der Exponential- und Quadratwurzelfunktion lediglich der x-Wert skaliert werden, d.h. hier reichen für den Skalierungsfaktor nach Gl. (2.35) zwei ROM-Speicher aus: constant K2_ROM_X1,K2_ROM_X2: MMX64_32:= -- hyperbolischer Modus ("10011000001001010000011001110010", -- 111111 "10010101101110101100100101100100", -- 111110 ... ... "00000000000000000000000000000000"); -- 000000

In Bild 4.6 fällt auf, dass nicht alle 56 Ergebnisstellen verwendet werden, sondern die niederwertigsten 26 Stellen wegfallen. Die nicht skizzierte Ablaufsteuerung wird im Prozess A_BEST des Top-Level-Moduls durch die Automatenzustände Z4 bis Z7 übernommen. In Z4 wird zunächst das Zweierkomplement in die Vorzeichen-Betrag-Darstellung umgewandelt. Lediglich bei der Berechnung der exp-Funktion müssen vor dieser Umwandlung die Register QX_INT (cosh) und QY_INT (sinh) addiert werden. Um in Z4 direkt auf die Ergebnisse der Umwandlung zugreifen zu können, werden diese den 32-Bit-Variablen VX_INT bzw. VY_INT zugewiesen. Zusätzlich erfolgt am Prozessende die Speicherung in QX_INT bzw. QY_INT. Das in Bild 4.6 aufgeführte Register wird im VHDL-Entwurf um ein Bit auf 33 Bit zur Aufnahme des Vorzeichens erweitert und mit QX_MULT bzw. QY_MULT bezeichnet. Die serielle Konstantenmultiplikation sieht dann für den zirkularen Modus (OP_CODE = X“01“) prinzipiell wie folgt aus: when Z4 => -- OP_CODE = X“01“ -- Umwandlung des Zweierkomplements in die Vorzeichen-Betrag-Darstellung if QX_INT(31) = '1' then VX_INT(31 downto 0) := QX_INT(31) & (not(QX_INT(30 downto 0)) + 1); else VX_INT(31 downto 0) := QX_INT; end if; if QY_INT(31) = '1' then VY_INT(31 downto 0) := QY_INT(31) & (not(QY_INT(30 downto 0)) + 1); else VY_INT(31 downto 0) := QY_INT; end if; QX_INT

34

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

S Vorkommastellen Nachkommastellen

31 30 2425 0

S OV H 23-Bit-Mantisse AB

31 30 25 24 02 1

a)

b)

c) AOV

drei zusätzliche Overflow-Bits

0

32

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

S Vorkommastellen Nachkommastellen

31 30 2425 0

S OV H 23-Bit-Mantisse AB

31 30 25 24 02 1

a)

b)

c) AOV

drei zusätzliche Overflow-Bits

0

32

4.6 Zusätzliche Operationen für ln-Berechnung Da der Logarithmus über den Areatangens hyperbolicus bestimmt wird, ist hier keine Skalierung der CORDIC-Berechnung notwendig. Allerdings ist der artanh nach Gl. (2.40) noch mit dem Faktor 2 zu multiplizieren und hierzu nach Gl. (2.58) noch der Faktor E·ln2 zu addieren. All diese Operationen werden in den für die Konstantenmultiplikation vorgesehenen Zuständen Z4 bis Z7 ausgeführt. Das Produkt E·ln2 ist als ROM-Speicher (K3_ROM) der Organisation 64*32 im Package-File deklariert. In Z4 wird die Multiplikation mit 2 durch einen Linksshift realisiert. Zusätzlich wird der Referenzexponent QREF_EXP_XY in die Vorzeichen-Betrag-Darstellung umgewandelt. In Z5 werden die 6 niederwertigsten Bits des Referenzexponenten zur Adressierung von K3_ROM verwendet und das 32-Bit-Datenwort im 33-Bit-Register QE_LN_2[32:0] gespeichert. Die Breite dieses Registers orientiert sich an QX_MULT bzw. QY_MULT und spielt eine Rolle bei der Normalisierung der Mantisse innerhalb der Endstufe im Zustand Z8. In Z5 wird das höchstwertige Bit QE_LN_2(32) auf ’0’ gesetzt. Die übrigen Bits repräsentieren das adressierte ROM-Datenwort. Ist der Referenzexponent negativ, dann erfolgt noch eine entsprechende Komplementbildung, wie der nachfolgende Codeausschnitt zeigt: when Z5 => -- OP_CODE = X’04’ if QREF_EXP_XY(9) = '1' then QE_LN_2

35

when Z6 => -- OP_CODE = X’04’ QE_LN_2(31 downto 0) -- OP_CODE = X’04’ if QE_LN_2(31) = '1' then QE_LN_2 von links eine ’1’ nachziehen QREF_EXP_XY

36

S OV H 23-Bit-Mantisse Addition-Bits

32 31 28 27 05 4

1

1

X

1

1

X X X X X X X X X X X X X X X X X X X X X X X X X

1

1

0

0

1

1

0

0 0 0

XX

X

X

X

X

X

X

X

X

X

0

0 0

0

0 00

0 0 0

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

0

0 0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0

0 0 0

0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0

-3

-2

-1

0

1

2

22

23

S OV H 23-Bit-Mantisse Addition-Bits

32 31 28 27 05 4

1

1

X

1

1

X X X X X X X X X X X X X X X X X X X X X X X X X

1

1

0

0

1

1

0

0 0 0

XX

X

X

X

X

X

X

X

X

X

0

0 0

0

0 00

0 0 0

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

0

0 0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0

0 0 0

0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0

-3

-2

-1

0

1

2

22

23

QX_MULT[32:0]

COUNT_SH_OPCountX_mult[9:0]

VREF_EXP_X[9:0]

QY_MULT[32:0]

COUNT_SH_OPCountY_mult[9:0]

QE_LN_2[32:0]

COUNT_SH_OPCountZ_add[9:0]

VREF_EXP_Y[9:0] VREF_EXP_Z[9:0]

Korrektur des Referenzexponenten

SH

FT

32L

RM

anZ

out[2

3:0]

QE

_LN

_2[3

2:0]

SH

FT

32L

RM

anY

out[2

3:0]

QY

_MU

LT[3

2:0]

SH

FT

32L

RM

anX

out[2

3:0]

QX

_MU

LT[3

2:0]

Man

Xou

t[22

:0]

Man

You

t[22:

0]M

anZ

out[2

2:0]

VREF_EXP_X[7:0] VREF_EXP_Y[7:0] VREF_EXP_Z[7:0]

QX_MULT[32]

QY_MULT[32]

QE_LN_2[32]

QQUOTIENT[9:0]

OP_CODE[7:0]

FALLUNTERSCHEIDUNG[1:0]

X“03“: exp-Funktion

QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] löschen;

X“04“: l