AVR Arithmetik/Sinus und Cosinus (CORDIC)
von gjlayde
| Größe | maximaler Verbrauch |
|---|---|
| Laufzeit | 560 Ticks (incl. CALL+RET) |
| RAM (statisch) | 0 Bytes |
| RAM (Stack) | 2–3 Bytes[1] |
| Flash | 278 Bytes |
CORDIC ist ein iterativer Algorithmus, mit dessen Hilfe sich Sinus und Cosinus effizient implementieren lassen. Dieser Artikel beschreibt eine CORDIC-Implementierung für 16-Bit Werte auf AVR.
Auf den eigentlichen Algorithmus wird in diesem Artikel nicht eingegangen. Für weiterführende Informationen siehe die Weblinks.
Die vorliegende Implementierung hat eine Optimierung in Bezug auf die Ausführungszeit zum Ziel.
Einsatzgebiet
Der Algorithmus findet dort Einsatz, wo eine Lookup-Tabelle zu groß oder zu ungenau ist. Mit 278 Bytes verbraucht CORDIC in etwa so viel Platz wie eine 256 Byte große Lookup-Tabelle mit 8-Bit Werten.
Der Algorithmus ist zwar in Assembler implementiert, kann aber problemlos zusammen mit avr-gcc verwendet werden. Er verwendet Befehle wie FMUL, setzt also einen AVR aus der ATmega-Familie voraus.
Der Algorithmus braucht für keine Berechnung (inclusive CALL und RET) länger als 560 Ticks.
Ein- und Ausgabe
Die Eingabe ist so skaliert, daß ein Grad aus 360 Inkrementen besteht. Gültige Eingaben liegen also in einem Bereich von ca. –91° bis 91° bzw. -91·360...91·360. Für die Berechnung von sin(45°) ist der Aufruf somit
// cos_sin.sin = sin (45°)
// cos_sin.cos = cos (45°)
cos_sin_t cos_sin = cordic (45*360);
Die Struktur cos_sin_t wird im Header cordic.h definiert und übergibt Sinus und Cosinus in einer Rückgabestruktur mit den 16-Bit Komponenten .sin und .cos.
Die Ausgabewerte liegen im 1.15 signed Q-Format vor, d.h. als Festkomma-Werte mit 15 binären Nachkommastellen und ohne Vorkommastelle. Das MSB gibt das Vorzeichen des Ergebnisses an. Ein Inkrement entspricht 1/32767. Dieser krumme Wert ist notwendig, um auch den Wert 1 ohne Überlauf darstellen zu können. Man kann die Ergebniswerte jedoch auch einfach als mit 1/32768 skaliert betrachten, was den Fehler bei der Berechnung dann um 1 Inkrement (also 1/32768) vergößert.
Die Ausgabe liegt – dargestellt als Hex-Zahl – im Bereich von 0x8001...0x7fff d.h. −0x7fff...0x7fff. Der problematische Wert 0x8000 wird nicht produziert.
Quellcode
cordic.h
Der Header regelt den Datenaustausch zwischen Assembler-Routine und C-Programm, indem er die Struktur cos_sin_t definiert und einen Prototypen für die cordic-Funktion deklariert.
#ifndef CORDIC_H
#define CORDIC_H
#include <stdint.h>
#include <math.h>
typedef struct
{
int16_t cos;
int16_t sin;
} cos_sin_t;
// Input in (-91°,91°) bzw. -91*360...91*360
// 1° = 360
#define CORDIC_SCALE_IN(X) ((180.*360./M_PI)*(X)+0.5)
// Output in (-1,1) als 1.15 signed Q-Format
#define CORDIC_SCALE_OUT(X) (32767.*(X))
// Gain = \prod_i \sqrt{1+4^{-i}}
// 1 / Gain
#define CORDIC_GAIN 0.60725294
extern cos_sin_t cordic (int16_t phi);
#endif // CORDIC_H
cordic.c
Dieses C-Modul dient lediglich dazu, ein Array mit den im Assembler benötigten Konstanten anzulegen. Das geht in C komfortabler als in Assembler.
Der CORDIC-Algorithmus bewirkt abhängig von der Anzahl der ausgeführten Iterationsschritte einen Zuwachs der Werte. Daher wird cos zu Beginn des Algorithmus nicht mit 1 initialisiert, sondern mit 1/GAIN wobei
- [math]\displaystyle{ \text{Gain} = \prod_i \sqrt{1+4^{-i}} }[/math]
Als kleine Optimierung könnte 1/GAIN direkt per LDI geladen werden anstatt ihn aus dem Array zu lesen.
#include <avr/pgmspace.h>
#include "cordic.h"
const int16_t cordic_atan_tab[1+16] PROGMEM =
{
// 1 / GAIN
CORDIC_SCALE_OUT (CORDIC_GAIN),
CORDIC_SCALE_IN (0.78539816), // atan (2^{-0})
CORDIC_SCALE_IN (0.46364761), // atan (2^{-1})
CORDIC_SCALE_IN (0.24497866), // atan (2^{-2})
CORDIC_SCALE_IN (0.12435499), // atan (2^{-3})
CORDIC_SCALE_IN (0.06241881), // atan (2^{-4})
CORDIC_SCALE_IN (0.03123983), // atan (2^{-5})
CORDIC_SCALE_IN (0.01562373), // atan (2^{-6})
CORDIC_SCALE_IN (0.00781234), // atan (2^{-7})
CORDIC_SCALE_IN (0.00390623), // atan (2^{-8})
CORDIC_SCALE_IN (0.00195312), // atan (2^{-9})
CORDIC_SCALE_IN (0.00097656), // atan (2^{-10})
CORDIC_SCALE_IN (0.00048828), // atan (2^{-11})
CORDIC_SCALE_IN (0.00024414), // atan (2^{-12})
CORDIC_SCALE_IN (0.00012207), // atan (2^{-13})
CORDIC_SCALE_IN (0.00006104), // atan (2^{-14})
CORDIC_SCALE_IN (0.00003052) // atan (2^{-15})
};
cordic.inc
Diesese Assembler Include-Datei definiert einige mehrfach verwendete Assembler-Makros. Es werden Assembler-Makros verwendet, da diese mächtiger sind als die gewohnten C-Makros.
;; a = 16 Bit signed
;; b = 8 Bit unsigned
;; c = 16 Bit signed
;; zero = 0
;; c += a*b / 2**8
.macro MADD16_8 c a b zero
fmulsu \a+1, \b
add \c, R0
adc \c+1, R1
fmul \a, \b
sbc \c+1, \zero
add \c, R1
adc \c+1, \zero
.endm
;; a = 16 Bit signed
;; b = 8 Bit unsigned
;; c = 16 Bit signed
;; zero = 0
;; c -= a*b / 2**8
.macro MSUB16_8 c a b zero
fmulsu \a+1, \b
sub \c, R0
sbc \c+1, R1
fmul \a, \b
adc \c+1, \zero
sub \c, R1
sbci \c+1, 0
.endm
;; a = 16 Bit signed
;; b = 8 Bit unsigned
;; c = 16 Bit signed
;; zero = 0
;; c += a*b / 2**16
.macro MADD16_16 c a b zero
fmulsu \a+1, \b
sbc \c+1, \zero
rol R0
adc \c, R1
adc \c+1, \zero
fmul \a, \b
adc \c, \zero
adc \c+1, \zero
.endm
;; a = 16 Bit signed
;; b = 8 Bit unsigned
;; c = 16 Bit signed
;; zero = 0
;; c -= a*b / 2**16
.macro MSUB16_16 c a b zero
fmulsu \a+1, \b
adc \c+1, \zero
rol R0
sbc \c, R1
sbci \c+1, 0
fmul \a, \b
sbci \c, 0
sbci \c+1, 0
.endm
.macro PHI_ADD
;; phi += atan (2^{-i})
lpm _tmp, z+
lpm _tmp+1, z+
add _phi, _tmp
adc _phi+1, _tmp+1
.endm
.macro PHI_SUB
;; phi -= atan (2^{-i})
lpm _tmp, z+
lpm _tmp+1, z+
sub _phi, _tmp
sbc _phi+1, _tmp+1
.endmcordic-asm.S
Der eigentliche CORDIC-Algorithmus. Er ist avr-gcc ABI-konform, kann also problemlos von avr-gcc aus verwendet werden, ohne mit dessen Registerverwendung zu kollidieren.
Der Algorithmus führt 16 Iterationsschritte von n=0...15 aus, wobei zur Erhöhung der Geschwindigkeit drei Fälle unterschieden werden:
- n = 0
- n = 1...7
- n = 8...15
Die im CORDIC üblichen Shifts um variable Offsets sind auf AVR ungünstig bzw. zweitraubend. Daher wird anstatt der Shifts mit Konstanten 2-n multipliziert.
;; CORDIC
;; 278 Bytes
;; Ticks < 560 (inclusive call+ret)
;; Georg-Johann Lay
.include "cordic.inc"
__zero_reg__ = 1
.text
.global cordic
.type cordic, @function
;; CORDIC-Algorithmus für cos() und sin()
;; Register-Verwendung ist avr-gcc ABI-konform
_phi = 26 ;; Winkel: 16 Bit signed
_cos = 22 ;; cos: 16 Bit signed R16..R22
_sin = 24 ;; sin: 16 Bit signed
_pow2 = 20 ;; 2^{-i}: 8 Bit unsigned R16..R23
_zero = 21 ;; immer 0: 8 Bit
_tmp = 18 ;; 16 Bit R16..R22
; R31:R30 ;; 16 Bit cordic_atan_tab[]
/*
cos = 1/1.6468 = 1 / CORDIC_GAIN
sin = 0
for i=0..15
tmp = sin
if (phi > 0)
sin += cos * 2^(-i)
cos -= tmp * 2^(-i)
phi -= arctan(2^(-i))
else
sin -= cos * 2^(-i)
cos += tmp * 2^(-i)
phi += arctan(2^(-i))
endif
next i
*/
cordic:
;; end prolog
movw _phi, R24
clr _zero
;; Z = cordic_atan_tab[]
ldi r30, lo8 (cordic_atan_tab)
ldi r31, hi8 (cordic_atan_tab)
;; cos = 1 / CORDIC_GAIN
lpm _cos, z+
lpm _cos+1, z+
ldi _pow2, 0x80
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Shift = 0
movw _sin, _cos
tst _phi+1
;; phi >= 0: sin = 0+cos
brpl .L_phi_pos_1_7
;; phi < 0: sin = 0-cos
com _sin+1
neg _sin
sbci _sin+1, -1
rjmp .L_phi_neg_1_7
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Shifts 1..7
.L_loop_1_7:
movw _tmp, _sin;
tst _phi+1
brmi 0f
;; phi > 0
MADD16_8 _sin _cos _pow2 _zero
MSUB16_8 _cos _tmp _pow2 _zero
.L_phi_pos_1_7:
PHI_SUB
lsr _pow2
brne .L_loop_1_7
breq .L_done_1_7
;; phi < 0
0: MSUB16_8 _sin _cos _pow2 _zero
MADD16_8 _cos _tmp _pow2 _zero
.L_phi_neg_1_7:
;; phi += atan (2^{-i})
PHI_ADD
lsr _pow2
brne .L_loop_1_7
.L_done_1_7:
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Shifts 8..15
ldi _pow2, 0x80
.L_loop_8_15:
movw _tmp, _sin;
tst _phi+1
brmi 0f
;; phi > 0
MADD16_16 _sin _cos _pow2 _zero
MSUB16_16 _cos _tmp _pow2 _zero
PHI_SUB
lsr _pow2
brne .L_loop_8_15
breq .L_done_8_15
;; phi < 0
0: MSUB16_16 _sin _cos _pow2 _zero
MADD16_16 _cos _tmp _pow2 _zero
PHI_ADD
lsr _pow2
brne .L_loop_8_15
.L_done_8_15:
;; avr-gcc ABI
clr __zero_reg__
;; Mittleren Fehler minimieren
;; cos += 1
subi _cos, lo8 (-1)
sbci _cos+1, hi8 (-1)
brvc 0f
;; Saturiere bei Überlauf
ldi _cos, lo8 (0x7fff)
ldi _cos+1, hi8 (0x7fff)
;, Epilog
0: ret
.size cordic, .-cordicGenauigkeit
| Größe | Standardabweichung | Maximaler Fehler |
|---|---|---|
| sin | 6.7·10-5 | 2.13·10-4 |
| cos | 6.1·10-5 | 2.13·10-4 |
| Radius | 5.8·10-5 | 2.13·10-4 |
| φ | 6.8·10-5 | 2.44·10-4 |
Die Güte der Ergebnisse kann aus der nebenstehenden Tabelle entnommen werden.
Die Ergebnisse sind recht gut normalverteilt um das korrekte Ergebnis. Der maximale absolute Fehler liegt bei 7 Inkrementen bzw. 0.0002, was einer effektiven Auflösung von 12 Bits entspricht. Mehr ist bei 16 Iterationsschritten auf 16-Bit Werte auch nicht zu erwarten, da in jedem Schritt Fehlerfortpflanzung und Rundung stattfinden.
Die untenstehenden Diagramme zeigen die Verteilungen der absoluten Fehler. In x-Richtung ist der Fehler aufgetragen, wieder in Einheiten von 1/32767. In y-Richtung ist die Anzahl der Werte aufgetragen, die mit dem entsprechenden absoluten Fehler behaftet sind.
- Standard- und Maximalabweichungen in Einheiten von 1/32767
Sinus
Cosinus
Radius
Winkel
Siehe auch
Anmerkungen
- ↑ abhängig von Byte-Größe des PC