Mikrocontroller.net - Benutzerbeiträge [de]

AVR Arithmetik/Sinus und Cosinus (Lineare Interpolation)

2012-04-22T21:43:12Z

95.208.157.192: Die -1 hatten mich erst auch verwirrt, daher meine Änderung. Dann hab ich's wieder rückgängig gemacht, weil eben Q-Format. Hiermit Versuch es etwas klarer zu formulieren

''von [[Benutzer:gjlayde|gjlayde]]''

<div style="float:right; margin:1em;">
{| {{Tabelle}}
|+ '''Lineare Interpolation: Ressourcenverbrauch'''
|- bgcolor="#ffddbb"
! Größe || maximaler Verbrauch
|-
| Laufzeit Sinus || 60 Ticks (incl. CALL+RET)
|-
| Laufzeit Cosinus || 74 Ticks (incl. CALL+RET)
|-
| RAM (statisch) || 0 Bytes
|-
| RAM (Stack) || 2–3 Bytes<ref name="pc">Abhängig von der Größe des PC</ref>
|-
| Flash || 300 Bytes<ref name="3.4.6">Übersetzt mit avr-gcc 3.4.6 für einen ATmega88</ref>
|-
| || 318 Bytes (avr-gcc 4.3.3 für ATmega88)
|}
</div>
Die AVR-Implementierung der linearen Interpolation von Sinus und Cosinus braucht mit 300 Bytes etwas mehr Platz als [[AVR Arithmetik/Sinus und Cosinus (CORDIC)|CORDIC]], arbeitet dafür aber deutlich schneller. Das Eingabeformat ist allerdings ein anderes als bei der obigen CORDIC-Implementierung. Die Werte streuen weniger als beim CORDIC, und sowohl Maximalfehler als auch Standardabweichungen sind kleiner.

Der Algorithmus ist in GNU-C implementiert, greift jedoch auf Inline-Assembler zurück. Dennoch dürfte er leichter zu portieren und zu verstehen sein als eine reine Assembler-Implementierung. Er verwendet Befehle wie MUL, setzt also einen AVR aus der ATmega-Familie voraus.

Der Algorithmus braucht für keine Berechnung (inclusive CALL und RET) länger als 60 Ticks<ref name="3.4.6"/> für den Sinus bzw. 74 Ticks<ref name="3.4.6"/> für den Cosinus. Für avr-gcc 4.x kommen u.U einige Ticks hinzu.

{{Clear}}

== Ein- und Ausgabe ==

Die Eingabe ist so skaliert, daß ein Grad aus 256 Inkrementen besteht. Gültige Eingaben liegen also in einem Bereich von 0° bis 255.99° bzw. 0·256...255·256+255. Für die Berechnung von sin(45°) ist der Aufruf beispielsweise
<c>
// si = sin (45°)
// co = cos (45°)
int16_t si = linsin (45*256);
int16_t co = lincos (45*256);
</c>

Für die Ausgabe gilt das gleiche wie für [[AVR Arithmetik/Sinus und Cosinus (CORDIC)#Ein- und Ausgabe|CORDIC]]: Es sind vorzeichenbehaftete Werte im 1.15 Q-Format, wobei der kritische Wert 0x8000 (entspricht −1 im Q-Format) nicht produziert wird. Kritisch u.A. deshalb, weil er nicht negiert werden kann.

== Quellcode ==

=== linsin.h ===

<c>
#ifndef LINSIN_H
#define LINSIN_H

#include <stdint.h>

int16_t linsin (uint16_t);
int16_t lincos (uint16_t);
#endif // LINSIN_H
</c>

=== linsin.c ===

{{Scrollbox|40ex;|
<c>
#include <avr/pgmspace.h>

#include "linsin.h"

// Nachguck-Tabelle für eine Sinus-Interpolante.
// Die Stützstellen sind volle °-Werte von 0°...90°.
// Die Werte sind so gewählt, daß die Polygonzug-Interpolante
// den Maximalfehler minimiert; daher sind die Werte in der
// Tabelle nicht die Funktionswerte des Sinus an den
// entsprechenden Stellen.
static const int16_t linsin_tab[91] PROGMEM =
{
0x0000, 0x023c, 0x0478, 0x06b3, 0x08ee, 0x0b28, 0x0d61, 0x0f9a, 0x11d1, 0x1406,
0x163a, 0x186d, 0x1a9d, 0x1ccb, 0x1ef8, 0x2121, 0x2348, 0x256d, 0x278e, 0x29ad,
0x2bc8, 0x2ddf, 0x2ff4, 0x3204, 0x3410, 0x3619, 0x381d, 0x3a1d, 0x3c18, 0x3e0f,
0x4001, 0x41ed, 0x43d5, 0x45b7, 0x4794, 0x496c, 0x4b3d, 0x4d09, 0x4ecf, 0x508e,
0x5248, 0x53fb, 0x55a7, 0x574d, 0x58eb, 0x5a83, 0x5c14, 0x5d9e, 0x5f20, 0x609b,
0x620f, 0x637a, 0x64df, 0x663b, 0x678f, 0x68db, 0x6a1f, 0x6b5b, 0x6c8e, 0x6db9,
0x6edb, 0x6ff5, 0x7106, 0x720e, 0x730d, 0x7403, 0x74f0, 0x75d4, 0x76af, 0x7781,
0x7849, 0x7908, 0x79bd, 0x7a69, 0x7b0c, 0x7ba5, 0x7c34, 0x7cb9, 0x7d35, 0x7da7,
0x7e0f, 0x7e6e, 0x7ec2, 0x7f0d, 0x7f4e, 0x7f85, 0x7fb1, 0x7fd4, 0x7fed, 0x7ffc,
0x7fff
};

// Ein Word per post-Increment aus dem Flash lesen.
// Dies ist besser als pgm_read_word (p++)
#define pgm_read_word_inc(addr) \
(__extension__({ \
uint16_t __result; \
__asm__ \
( \
"lpm %A0, Z+" "\n\t" \
"lpm %B0, Z+" \
: "=r" (__result), "+z" (addr) \
); \
__result; \
}))

// Multiply-Accumulate
// return c + a*b
// c = 1.15 signed Q-Format
// a = 1.15 signed Q-Format
// b = 0.8 unsigned Q-Format
static inline int16_t fmac16_8 (int16_t c, int16_t a, uint8_t b)
{
asm ("; fmac16_8: %B[c]:%A[c] += %B[a]:%A[a] * %[b]" "\n\t"
"mulsu %B[a], %[b]" "\n\t" // ( (signed)ah * (unsigned)b ) << 1
"add %A[c], R0" "\n\t"
"adc %B[c], R1" "\n\t"
"mul %[b], %A[a]" "\n\t" // ( (unsigned) b * al ) << 1
"add %A[c], R1" "\n\t"
"clr __zero_reg__" "\n\t"
"adc %B[c], __zero_reg__"
: [c] "+r" (c)
: [a] "a" (a), [b] "a" (b)
);

return c;
}

// Sinus
// phi in [0°,256°) wobei 1° = 256
int16_t linsin (uint16_t phi)
{
uint8_t neg = 0;
uint16_t _180_grad = (uint16_t) 180*256;
// Wir wollen diese Konstante in einem Register
asm ("" : "+r" (_180_grad));

if (phi >= _180_grad)
{
// sin(x) = -sin(x+pi)
// phi in [0°,180°)
phi -= _180_grad;
neg = 1;
}

if (phi > 90*256)
// sin(x) = sin(pi-x)
// phi in [0°,90°]
phi = _180_grad - phi;

// Stützpunkte aus Tabelle lesen.
// Für phi = 90° ist das formal nicht 100% korrekt,
// da über das Tabellenende hinausgelesen wird.
// In diesem Falle ist aber phi.lo=0, so daß unten
// dann data1-data0 mit 0 multipliziert wird.
// Ausserdem gibt's auf AVR keine Segmentation Faults ;-)
const int16_t * p = & linsin_tab[phi >> 8];
uint16_t data0 = pgm_read_word_inc (p);
uint16_t data1 = pgm_read_word (p);

// sin = data0 + (data1-data0) * phi_lo
int16_t si = fmac16_8 (data0, data1-data0, phi);

return neg ? -si : si;
}

// Cosinus
// phi in [0°,256°) wobei 1° = 256
int16_t lincos (uint16_t phi)
{
uint16_t _90_grad = (uint16_t) 90*256;
// Wir wollen diese Konstante in einem Register
asm ("" : "+r" (_90_grad));

if (phi <= _90_grad)
// cos(x) = sin(pi/2 - x)
return linsin (_90_grad - phi);
else
// cos(x) = -sin (x - pi/2)
return -linsin (phi - _90_grad);
}
</c>
}}

=== Aufbau der Tabelle ===

Zur Erstellung der Tabelle gehen wir von äquidistanten Stützstellen im Abstand δ aus. Entwickeln wir die zu interpolierende Funktion in eine [http://de.wikipedia.org/wiki/Taylorreihe Taylor-Reihe] um den Mittelpunkt eines Teilintervalles der Breite δ, so erhalten wir für den absoluten Fehler
:<math>
\Delta = \frac{f''(\xi)}2 \cdot \left(\frac\delta 2\right)^2 + \mathcal O(f'''\cdot\delta^3)
</math>
Indem wir durch eine Sekante annähern anstatt durch – wie durch die Taylor-Reihe gegeben – eine Tangente, können wir die Größenordnung des Fehlers halbieren und erhalten:
:<math>
|\Delta| \approx \tfrac1{16}\cdot |f''|\cdot \delta^2
\;\leqslant\; \tfrac1{16} \cdot \delta^2
\;\stackrel!\leqslant\; 2^{-15}
</math>
Weil mit der Tabelle das Intervall [0..π/2] interpoliert werden soll, brauchen wir mindestens
:<math>
16\pi\sqrt{2} \approx 71
</math>
Teilintervalle. Wir machen also die naheliegende Wahl von 90 Teilintervallen bzw. 91 Stützstellen von 0°...90° und damit
:<math>
\delta = \tfrac1{90} \cdot \tfrac\pi2 \approx 0.01745 = 571.9\cdot 2^{-15}
</math>
Damit wäre ein Eingabeformat von 512=29 Inkrementen pro Grad zwar naheliegend (die Ableitung des Sinus/Cosinus ist höchstens 1), aus rechentechnischen Erwägungen unterteilen wir ein Grad jedoch nur in 256=28 Inkremente und kommen damit zum vorliegenden Eingabeformat.

== Genauigkeit ==

<div style="float:right; margin:1em;">
{| {{Tabelle|text-align:center;}}
|+ '''Lineare Interpolation: Genauigkeit'''
|- bgcolor="#ffddbb"
! Größe || Standardabweichung || Maximaler Fehler
|-
! sin
| 2.2·10-5 || 6.1·10-5
|-
! cos
| 2.3·10-5 || 6.1·10-5
|-
! Radius
| 2.9·10-5 || 6.1·10-5
|-
! φ
| 1.5·10-5 || 3.1·10-5
|}
</div>
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 2 Inkrementen bzw. 0.00006, was einer effektiven Auflösung von 14  Nachkomma-Bits bzw. insgesamt 15 Bits entspricht.

Im Gegensatz zum CORDIC sind die Ergebnisse der linearen Interpolation per Konstruktion ''monoton'': Im Intervall 0...π/2 ist die Sinus-Funktion z. B. monoton wachsend, und das gilt auch für die hier berechneten Werte. Bei CORDIC sind die Ergebnisse "verrauscht", so daß es auch bei etwas größerer Eingabe zu einem kleineren Wert kommen kann.

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.

<gallery perrow="2" caption="Standard- und Maximalabweichungen in Einheiten von  1/32767" widths="260" heights="130">
Bild:Linsin-dsin.png|Sinus
Bild:Linsin-dcos.png|Cosinus
Bild:Linsin-dr.png|Radius
Bild:Linsin-dphi.png|Winkel
</gallery>

== Siehe auch ==
* [[AVR-CORDIC]]
* [[AVR Arithmetik/Sinus und Cosinus (CORDIC)|Sinus und Cosinus mit CORDIC]]

== Anmerkungen ==
<references/>

[[Kategorie:AVR-Arithmetik|S]]

AVR Arithmetik/Sinus und Cosinus (Lineare Interpolation)

2012-04-22T15:55:52Z

95.208.157.192: Änderung 65864 von 95.208.157.192 (Diskussion) wurde rückgängig gemacht.

''von [[Benutzer:gjlayde|gjlayde]]''

<div style="float:right; margin:1em;">
{| {{Tabelle}}
|+ '''Lineare Interpolation: Ressourcenverbrauch'''
|- bgcolor="#ffddbb"
! Größe || maximaler Verbrauch
|-
| Laufzeit Sinus || 60 Ticks (incl. CALL+RET)
|-
| Laufzeit Cosinus || 74 Ticks (incl. CALL+RET)
|-
| RAM (statisch) || 0 Bytes
|-
| RAM (Stack) || 2–3 Bytes<ref name="pc">Abhängig von der Größe des PC</ref>
|-
| Flash || 300 Bytes<ref name="3.4.6">Übersetzt mit avr-gcc 3.4.6 für einen ATmega88</ref>
|-
| || 318 Bytes (avr-gcc 4.3.3 für ATmega88)
|}
</div>
Die AVR-Implementierung der linearen Interpolation von Sinus und Cosinus braucht mit 300 Bytes etwas mehr Platz als [[AVR Arithmetik/Sinus und Cosinus (CORDIC)|CORDIC]], arbeitet dafür aber deutlich schneller. Das Eingabeformat ist allerdings ein anderes als bei der obigen CORDIC-Implementierung. Die Werte streuen weniger als beim CORDIC, und sowohl Maximalfehler als auch Standardabweichungen sind kleiner.

Der Algorithmus ist in GNU-C implementiert, greift jedoch auf Inline-Assembler zurück. Dennoch dürfte er leichter zu portieren und zu verstehen sein als eine reine Assembler-Implementierung. Er verwendet Befehle wie MUL, setzt also einen AVR aus der ATmega-Familie voraus.

Der Algorithmus braucht für keine Berechnung (inclusive CALL und RET) länger als 60 Ticks<ref name="3.4.6"/> für den Sinus bzw. 74 Ticks<ref name="3.4.6"/> für den Cosinus. Für avr-gcc 4.x kommen u.U einige Ticks hinzu.

{{Clear}}

== Ein- und Ausgabe ==

Die Eingabe ist so skaliert, daß ein Grad aus 256 Inkrementen besteht. Gültige Eingaben liegen also in einem Bereich von 0° bis 255.99° bzw. 0·256...255·256+255. Für die Berechnung von sin(45°) ist der Aufruf beispielsweise
<c>
// si = sin (45°)
// co = cos (45°)
int16_t si = linsin (45*256);
int16_t co = lincos (45*256);
</c>

Für die Ausgabe gilt das gleiche wie für [[AVR Arithmetik/Sinus und Cosinus (CORDIC)#Ein- und Ausgabe|CORDIC]]: Es sind vorzeichenbehaftete Werte im 1.15 Q-Format, wobei der kritische Wert 0x8000 bzw. −1 nicht produziert wird.

== Quellcode ==

=== linsin.h ===

<c>
#ifndef LINSIN_H
#define LINSIN_H

#include <stdint.h>

int16_t linsin (uint16_t);
int16_t lincos (uint16_t);
#endif // LINSIN_H
</c>

=== linsin.c ===

{{Scrollbox|40ex;|
<c>
#include <avr/pgmspace.h>

#include "linsin.h"

// Nachguck-Tabelle für eine Sinus-Interpolante.
// Die Stützstellen sind volle °-Werte von 0°...90°.
// Die Werte sind so gewählt, daß die Polygonzug-Interpolante
// den Maximalfehler minimiert; daher sind die Werte in der
// Tabelle nicht die Funktionswerte des Sinus an den
// entsprechenden Stellen.
static const int16_t linsin_tab[91] PROGMEM =
{
0x0000, 0x023c, 0x0478, 0x06b3, 0x08ee, 0x0b28, 0x0d61, 0x0f9a, 0x11d1, 0x1406,
0x163a, 0x186d, 0x1a9d, 0x1ccb, 0x1ef8, 0x2121, 0x2348, 0x256d, 0x278e, 0x29ad,
0x2bc8, 0x2ddf, 0x2ff4, 0x3204, 0x3410, 0x3619, 0x381d, 0x3a1d, 0x3c18, 0x3e0f,
0x4001, 0x41ed, 0x43d5, 0x45b7, 0x4794, 0x496c, 0x4b3d, 0x4d09, 0x4ecf, 0x508e,
0x5248, 0x53fb, 0x55a7, 0x574d, 0x58eb, 0x5a83, 0x5c14, 0x5d9e, 0x5f20, 0x609b,
0x620f, 0x637a, 0x64df, 0x663b, 0x678f, 0x68db, 0x6a1f, 0x6b5b, 0x6c8e, 0x6db9,
0x6edb, 0x6ff5, 0x7106, 0x720e, 0x730d, 0x7403, 0x74f0, 0x75d4, 0x76af, 0x7781,
0x7849, 0x7908, 0x79bd, 0x7a69, 0x7b0c, 0x7ba5, 0x7c34, 0x7cb9, 0x7d35, 0x7da7,
0x7e0f, 0x7e6e, 0x7ec2, 0x7f0d, 0x7f4e, 0x7f85, 0x7fb1, 0x7fd4, 0x7fed, 0x7ffc,
0x7fff
};

// Ein Word per post-Increment aus dem Flash lesen.
// Dies ist besser als pgm_read_word (p++)
#define pgm_read_word_inc(addr) \
(__extension__({ \
uint16_t __result; \
__asm__ \
( \
"lpm %A0, Z+" "\n\t" \
"lpm %B0, Z+" \
: "=r" (__result), "+z" (addr) \
); \
__result; \
}))

// Multiply-Accumulate
// return c + a*b
// c = 1.15 signed Q-Format
// a = 1.15 signed Q-Format
// b = 0.8 unsigned Q-Format
static inline int16_t fmac16_8 (int16_t c, int16_t a, uint8_t b)
{
asm ("; fmac16_8: %B[c]:%A[c] += %B[a]:%A[a] * %[b]" "\n\t"
"mulsu %B[a], %[b]" "\n\t" // ( (signed)ah * (unsigned)b ) << 1
"add %A[c], R0" "\n\t"
"adc %B[c], R1" "\n\t"
"mul %[b], %A[a]" "\n\t" // ( (unsigned) b * al ) << 1
"add %A[c], R1" "\n\t"
"clr __zero_reg__" "\n\t"
"adc %B[c], __zero_reg__"
: [c] "+r" (c)
: [a] "a" (a), [b] "a" (b)
);

return c;
}

// Sinus
// phi in [0°,256°) wobei 1° = 256
int16_t linsin (uint16_t phi)
{
uint8_t neg = 0;
uint16_t _180_grad = (uint16_t) 180*256;
// Wir wollen diese Konstante in einem Register
asm ("" : "+r" (_180_grad));

if (phi >= _180_grad)
{
// sin(x) = -sin(x+pi)
// phi in [0°,180°)
phi -= _180_grad;
neg = 1;
}

if (phi > 90*256)
// sin(x) = sin(pi-x)
// phi in [0°,90°]
phi = _180_grad - phi;

// Stützpunkte aus Tabelle lesen.
// Für phi = 90° ist das formal nicht 100% korrekt,
// da über das Tabellenende hinausgelesen wird.
// In diesem Falle ist aber phi.lo=0, so daß unten
// dann data1-data0 mit 0 multipliziert wird.
// Ausserdem gibt's auf AVR keine Segmentation Faults ;-)
const int16_t * p = & linsin_tab[phi >> 8];
uint16_t data0 = pgm_read_word_inc (p);
uint16_t data1 = pgm_read_word (p);

// sin = data0 + (data1-data0) * phi_lo
int16_t si = fmac16_8 (data0, data1-data0, phi);

return neg ? -si : si;
}

// Cosinus
// phi in [0°,256°) wobei 1° = 256
int16_t lincos (uint16_t phi)
{
uint16_t _90_grad = (uint16_t) 90*256;
// Wir wollen diese Konstante in einem Register
asm ("" : "+r" (_90_grad));

if (phi <= _90_grad)
// cos(x) = sin(pi/2 - x)
return linsin (_90_grad - phi);
else
// cos(x) = -sin (x - pi/2)
return -linsin (phi - _90_grad);
}
</c>
}}

=== Aufbau der Tabelle ===

Zur Erstellung der Tabelle gehen wir von äquidistanten Stützstellen im Abstand δ aus. Entwickeln wir die zu interpolierende Funktion in eine [http://de.wikipedia.org/wiki/Taylorreihe Taylor-Reihe] um den Mittelpunkt eines Teilintervalles der Breite δ, so erhalten wir für den absoluten Fehler
:<math>
\Delta = \frac{f''(\xi)}2 \cdot \left(\frac\delta 2\right)^2 + \mathcal O(f'''\cdot\delta^3)
</math>
Indem wir durch eine Sekante annähern anstatt durch – wie durch die Taylor-Reihe gegeben – eine Tangente, können wir die Größenordnung des Fehlers halbieren und erhalten:
:<math>
|\Delta| \approx \tfrac1{16}\cdot |f''|\cdot \delta^2
\;\leqslant\; \tfrac1{16} \cdot \delta^2
\;\stackrel!\leqslant\; 2^{-15}
</math>
Weil mit der Tabelle das Intervall [0..π/2] interpoliert werden soll, brauchen wir mindestens
:<math>
16\pi\sqrt{2} \approx 71
</math>
Teilintervalle. Wir machen also die naheliegende Wahl von 90 Teilintervallen bzw. 91 Stützstellen von 0°...90° und damit
:<math>
\delta = \tfrac1{90} \cdot \tfrac\pi2 \approx 0.01745 = 571.9\cdot 2^{-15}
</math>
Damit wäre ein Eingabeformat von 512=29 Inkrementen pro Grad zwar naheliegend (die Ableitung des Sinus/Cosinus ist höchstens 1), aus rechentechnischen Erwägungen unterteilen wir ein Grad jedoch nur in 256=28 Inkremente und kommen damit zum vorliegenden Eingabeformat.

== Genauigkeit ==

<div style="float:right; margin:1em;">
{| {{Tabelle|text-align:center;}}
|+ '''Lineare Interpolation: Genauigkeit'''
|- bgcolor="#ffddbb"
! Größe || Standardabweichung || Maximaler Fehler
|-
! sin
| 2.2·10-5 || 6.1·10-5
|-
! cos
| 2.3·10-5 || 6.1·10-5
|-
! Radius
| 2.9·10-5 || 6.1·10-5
|-
! φ
| 1.5·10-5 || 3.1·10-5
|}
</div>
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 2 Inkrementen bzw. 0.00006, was einer effektiven Auflösung von 14  Nachkomma-Bits bzw. insgesamt 15 Bits entspricht.

Im Gegensatz zum CORDIC sind die Ergebnisse der linearen Interpolation per Konstruktion ''monoton'': Im Intervall 0...π/2 ist die Sinus-Funktion z. B. monoton wachsend, und das gilt auch für die hier berechneten Werte. Bei CORDIC sind die Ergebnisse "verrauscht", so daß es auch bei etwas größerer Eingabe zu einem kleineren Wert kommen kann.

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.

<gallery perrow="2" caption="Standard- und Maximalabweichungen in Einheiten von  1/32767" widths="260" heights="130">
Bild:Linsin-dsin.png|Sinus
Bild:Linsin-dcos.png|Cosinus
Bild:Linsin-dr.png|Radius
Bild:Linsin-dphi.png|Winkel
</gallery>

== Siehe auch ==
* [[AVR-CORDIC]]
* [[AVR Arithmetik/Sinus und Cosinus (CORDIC)|Sinus und Cosinus mit CORDIC]]

== Anmerkungen ==
<references/>

[[Kategorie:AVR-Arithmetik|S]]

AVR Arithmetik/Sinus und Cosinus (Lineare Interpolation)

2012-04-22T15:50:27Z

95.208.157.192: /* Ein- und Ausgabe */

''von [[Benutzer:gjlayde|gjlayde]]''

<div style="float:right; margin:1em;">
{| {{Tabelle}}
|+ '''Lineare Interpolation: Ressourcenverbrauch'''
|- bgcolor="#ffddbb"
! Größe || maximaler Verbrauch
|-
| Laufzeit Sinus || 60 Ticks (incl. CALL+RET)
|-
| Laufzeit Cosinus || 74 Ticks (incl. CALL+RET)
|-
| RAM (statisch) || 0 Bytes
|-
| RAM (Stack) || 2–3 Bytes<ref name="pc">Abhängig von der Größe des PC</ref>
|-
| Flash || 300 Bytes<ref name="3.4.6">Übersetzt mit avr-gcc 3.4.6 für einen ATmega88</ref>
|-
| || 318 Bytes (avr-gcc 4.3.3 für ATmega88)
|}
</div>
Die AVR-Implementierung der linearen Interpolation von Sinus und Cosinus braucht mit 300 Bytes etwas mehr Platz als [[AVR Arithmetik/Sinus und Cosinus (CORDIC)|CORDIC]], arbeitet dafür aber deutlich schneller. Das Eingabeformat ist allerdings ein anderes als bei der obigen CORDIC-Implementierung. Die Werte streuen weniger als beim CORDIC, und sowohl Maximalfehler als auch Standardabweichungen sind kleiner.

Der Algorithmus ist in GNU-C implementiert, greift jedoch auf Inline-Assembler zurück. Dennoch dürfte er leichter zu portieren und zu verstehen sein als eine reine Assembler-Implementierung. Er verwendet Befehle wie MUL, setzt also einen AVR aus der ATmega-Familie voraus.

Der Algorithmus braucht für keine Berechnung (inclusive CALL und RET) länger als 60 Ticks<ref name="3.4.6"/> für den Sinus bzw. 74 Ticks<ref name="3.4.6"/> für den Cosinus. Für avr-gcc 4.x kommen u.U einige Ticks hinzu.

{{Clear}}

== Ein- und Ausgabe ==

Die Eingabe ist so skaliert, daß ein Grad aus 256 Inkrementen besteht. Gültige Eingaben liegen also in einem Bereich von 0° bis 255.99° bzw. 0·256...255·256+255. Für die Berechnung von sin(45°) ist der Aufruf beispielsweise
<c>
// si = sin (45°)
// co = cos (45°)
int16_t si = linsin (45*256);
int16_t co = lincos (45*256);
</c>

Für die Ausgabe gilt das gleiche wie für [[AVR Arithmetik/Sinus und Cosinus (CORDIC)#Ein- und Ausgabe|CORDIC]]: Es sind vorzeichenbehaftete Werte im 1.15 Q-Format, wobei der kritische Wert 0x8000 bzw. −32768 nicht produziert wird.

== Quellcode ==

=== linsin.h ===

<c>
#ifndef LINSIN_H
#define LINSIN_H

#include <stdint.h>

int16_t linsin (uint16_t);
int16_t lincos (uint16_t);
#endif // LINSIN_H
</c>

=== linsin.c ===

{{Scrollbox|40ex;|
<c>
#include <avr/pgmspace.h>

#include "linsin.h"

// Nachguck-Tabelle für eine Sinus-Interpolante.
// Die Stützstellen sind volle °-Werte von 0°...90°.
// Die Werte sind so gewählt, daß die Polygonzug-Interpolante
// den Maximalfehler minimiert; daher sind die Werte in der
// Tabelle nicht die Funktionswerte des Sinus an den
// entsprechenden Stellen.
static const int16_t linsin_tab[91] PROGMEM =
{
0x0000, 0x023c, 0x0478, 0x06b3, 0x08ee, 0x0b28, 0x0d61, 0x0f9a, 0x11d1, 0x1406,
0x163a, 0x186d, 0x1a9d, 0x1ccb, 0x1ef8, 0x2121, 0x2348, 0x256d, 0x278e, 0x29ad,
0x2bc8, 0x2ddf, 0x2ff4, 0x3204, 0x3410, 0x3619, 0x381d, 0x3a1d, 0x3c18, 0x3e0f,
0x4001, 0x41ed, 0x43d5, 0x45b7, 0x4794, 0x496c, 0x4b3d, 0x4d09, 0x4ecf, 0x508e,
0x5248, 0x53fb, 0x55a7, 0x574d, 0x58eb, 0x5a83, 0x5c14, 0x5d9e, 0x5f20, 0x609b,
0x620f, 0x637a, 0x64df, 0x663b, 0x678f, 0x68db, 0x6a1f, 0x6b5b, 0x6c8e, 0x6db9,
0x6edb, 0x6ff5, 0x7106, 0x720e, 0x730d, 0x7403, 0x74f0, 0x75d4, 0x76af, 0x7781,
0x7849, 0x7908, 0x79bd, 0x7a69, 0x7b0c, 0x7ba5, 0x7c34, 0x7cb9, 0x7d35, 0x7da7,
0x7e0f, 0x7e6e, 0x7ec2, 0x7f0d, 0x7f4e, 0x7f85, 0x7fb1, 0x7fd4, 0x7fed, 0x7ffc,
0x7fff
};

// Ein Word per post-Increment aus dem Flash lesen.
// Dies ist besser als pgm_read_word (p++)
#define pgm_read_word_inc(addr) \
(__extension__({ \
uint16_t __result; \
__asm__ \
( \
"lpm %A0, Z+" "\n\t" \
"lpm %B0, Z+" \
: "=r" (__result), "+z" (addr) \
); \
__result; \
}))

// Multiply-Accumulate
// return c + a*b
// c = 1.15 signed Q-Format
// a = 1.15 signed Q-Format
// b = 0.8 unsigned Q-Format
static inline int16_t fmac16_8 (int16_t c, int16_t a, uint8_t b)
{
asm ("; fmac16_8: %B[c]:%A[c] += %B[a]:%A[a] * %[b]" "\n\t"
"mulsu %B[a], %[b]" "\n\t" // ( (signed)ah * (unsigned)b ) << 1
"add %A[c], R0" "\n\t"
"adc %B[c], R1" "\n\t"
"mul %[b], %A[a]" "\n\t" // ( (unsigned) b * al ) << 1
"add %A[c], R1" "\n\t"
"clr __zero_reg__" "\n\t"
"adc %B[c], __zero_reg__"
: [c] "+r" (c)
: [a] "a" (a), [b] "a" (b)
);

return c;
}

// Sinus
// phi in [0°,256°) wobei 1° = 256
int16_t linsin (uint16_t phi)
{
uint8_t neg = 0;
uint16_t _180_grad = (uint16_t) 180*256;
// Wir wollen diese Konstante in einem Register
asm ("" : "+r" (_180_grad));

if (phi >= _180_grad)
{
// sin(x) = -sin(x+pi)
// phi in [0°,180°)
phi -= _180_grad;
neg = 1;
}

if (phi > 90*256)
// sin(x) = sin(pi-x)
// phi in [0°,90°]
phi = _180_grad - phi;

// Stützpunkte aus Tabelle lesen.
// Für phi = 90° ist das formal nicht 100% korrekt,
// da über das Tabellenende hinausgelesen wird.
// In diesem Falle ist aber phi.lo=0, so daß unten
// dann data1-data0 mit 0 multipliziert wird.
// Ausserdem gibt's auf AVR keine Segmentation Faults ;-)
const int16_t * p = & linsin_tab[phi >> 8];
uint16_t data0 = pgm_read_word_inc (p);
uint16_t data1 = pgm_read_word (p);

// sin = data0 + (data1-data0) * phi_lo
int16_t si = fmac16_8 (data0, data1-data0, phi);

return neg ? -si : si;
}

// Cosinus
// phi in [0°,256°) wobei 1° = 256
int16_t lincos (uint16_t phi)
{
uint16_t _90_grad = (uint16_t) 90*256;
// Wir wollen diese Konstante in einem Register
asm ("" : "+r" (_90_grad));

if (phi <= _90_grad)
// cos(x) = sin(pi/2 - x)
return linsin (_90_grad - phi);
else
// cos(x) = -sin (x - pi/2)
return -linsin (phi - _90_grad);
}
</c>
}}

=== Aufbau der Tabelle ===

Zur Erstellung der Tabelle gehen wir von äquidistanten Stützstellen im Abstand δ aus. Entwickeln wir die zu interpolierende Funktion in eine [http://de.wikipedia.org/wiki/Taylorreihe Taylor-Reihe] um den Mittelpunkt eines Teilintervalles der Breite δ, so erhalten wir für den absoluten Fehler
:<math>
\Delta = \frac{f''(\xi)}2 \cdot \left(\frac\delta 2\right)^2 + \mathcal O(f'''\cdot\delta^3)
</math>
Indem wir durch eine Sekante annähern anstatt durch – wie durch die Taylor-Reihe gegeben – eine Tangente, können wir die Größenordnung des Fehlers halbieren und erhalten:
:<math>
|\Delta| \approx \tfrac1{16}\cdot |f''|\cdot \delta^2
\;\leqslant\; \tfrac1{16} \cdot \delta^2
\;\stackrel!\leqslant\; 2^{-15}
</math>
Weil mit der Tabelle das Intervall [0..π/2] interpoliert werden soll, brauchen wir mindestens
:<math>
16\pi\sqrt{2} \approx 71
</math>
Teilintervalle. Wir machen also die naheliegende Wahl von 90 Teilintervallen bzw. 91 Stützstellen von 0°...90° und damit
:<math>
\delta = \tfrac1{90} \cdot \tfrac\pi2 \approx 0.01745 = 571.9\cdot 2^{-15}
</math>
Damit wäre ein Eingabeformat von 512=29 Inkrementen pro Grad zwar naheliegend (die Ableitung des Sinus/Cosinus ist höchstens 1), aus rechentechnischen Erwägungen unterteilen wir ein Grad jedoch nur in 256=28 Inkremente und kommen damit zum vorliegenden Eingabeformat.

== Genauigkeit ==

<div style="float:right; margin:1em;">
{| {{Tabelle|text-align:center;}}
|+ '''Lineare Interpolation: Genauigkeit'''
|- bgcolor="#ffddbb"
! Größe || Standardabweichung || Maximaler Fehler
|-
! sin
| 2.2·10-5 || 6.1·10-5
|-
! cos
| 2.3·10-5 || 6.1·10-5
|-
! Radius
| 2.9·10-5 || 6.1·10-5
|-
! φ
| 1.5·10-5 || 3.1·10-5
|}
</div>
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 2 Inkrementen bzw. 0.00006, was einer effektiven Auflösung von 14  Nachkomma-Bits bzw. insgesamt 15 Bits entspricht.

Im Gegensatz zum CORDIC sind die Ergebnisse der linearen Interpolation per Konstruktion ''monoton'': Im Intervall 0...π/2 ist die Sinus-Funktion z. B. monoton wachsend, und das gilt auch für die hier berechneten Werte. Bei CORDIC sind die Ergebnisse "verrauscht", so daß es auch bei etwas größerer Eingabe zu einem kleineren Wert kommen kann.

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.

<gallery perrow="2" caption="Standard- und Maximalabweichungen in Einheiten von  1/32767" widths="260" heights="130">
Bild:Linsin-dsin.png|Sinus
Bild:Linsin-dcos.png|Cosinus
Bild:Linsin-dr.png|Radius
Bild:Linsin-dphi.png|Winkel
</gallery>

== Siehe auch ==
* [[AVR-CORDIC]]
* [[AVR Arithmetik/Sinus und Cosinus (CORDIC)|Sinus und Cosinus mit CORDIC]]

== Anmerkungen ==
<references/>

[[Kategorie:AVR-Arithmetik|S]]