AVR Arithmetik/Sinus und Cosinus (CORDIC)

Wechseln zu: Navigation, Suche

von gjlayde

CORDIC: Ressourcenverbrauch
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

\text{Gain} = \prod_i \sqrt{1+4^{-i}}

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
.endm

cordic-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, .-cordic

Genauigkeit

CORDIC: Genauigkeit
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.

Siehe auch

Anmerkungen

  1. abhängig von Byte-Größe des PC

Weblinks