Auf dieser Seite gibt's die Erläuterungen zur Sinusfunktion für Embedded-Anwendungen, die hier angegeben ist.

Zunächst der Queltext, den ich von Kommentaren befreit habe (die Erläuterungen decken dies bereits ab). Die Funktion besteht aus 8 Zeilen Code, hinzu kommt eine Wertetabelle und ein paar Präprozessordefinitionen, die lediglich die Lesbarkeit erhöhen sollen und schnelle Änderung der Tabellengröße ermöglichen (dazu muss man nur TABLE_BITS ändern - aber Vorsicht: Mit jedem Wert verdoppelt sich der Speicherplatz für die Tabelle). Für das Beispiel ist dieser Wert auf 5 eingestellt, damit erhält man eine Genauigleit von ca einem halben Promill. Auch die Variablengröße kann geändert werden, Z.B. indem int16_t durch int8_t ersetzt wird. In beiden Fällen muss natürlich die Tabelle auch neuberechnet werden.

Erläuterung

Das interessante passiert in der sin1()-Funktion, in der wir den Sinus mit Hilfe von einer Tabelle berechnen. Sie arbeitet mit Festkommaarithmetik, und erhält die Eingabewinkel normiert auf 1, nicht auf 2*PI wie die normale sin()-Funktion. Die Zahlen sind 16 Bit groß und das Komma sitzt an der 15ten Stelle (Q15). Die Ausgabe ist ebenfalls 16 Bit Q15, von -1 bis 1.

Mit der gewählten Eingabeskalierung wird die Auflösung des Winkels bestmöglich genutzt, und zusätzlich kommt sie vielen Anwendungen entgegen. Beispielsweise gibt es viele Winkelgeber, die 2^N Werte pro Umdrehung haben. Der Maximalwert der sin1()-Funktion ist 32767, und entspricht auch einer Umdrehung. Daher müssen die Bits der Winkelgeberwerte nur ein wenig geschoben werden bevor man den Sinus berechnen kann.

Tabellen

sin90 sagt bereits aus, dass sie Ergebnisse zwischen 0˚ und 90˚ enthält. Weil die Kurve symmetrisch ist können wir die Ergebnisse für alle anderen Eingabewinkel schnell berechnen. Damit sparen wir 75% Speicherbelegung. Die Tabelle ist 33 Werte groß, also "5 Bit + 1" (2^5=32). Der zusätzliche Eintrag ist nur für die Interpolation enthält den Wert 1 (bzw. "fast 1"). An den Stützstellen stimmt die Tabelle genau mit dem echten Sinus überein (bzw. so genau wie die Aufösung 16 Bit-Zahl ist). Wie Werte liegen also bei den Winkeln 0˚, 11.25˚, 22.5˚, usw. Die darin stehenden Zahlen sind keine Fließkommazahlen, sondern Ganzzahlen, bei denen wir uns das Komma ein eine feste Stelle denken, und zwar an die 15te (das 16te Bit ist schon das Vorzeichenbit). Genaueres über Festkommaarithmetik habe ich hier. Die Ausgabewerte der Tabelle sind also Q15 von -1 bis 1, und wir nehmen damit die maximale Genauigkeit aus einer 16 Bit-Zahl heraus. Praktisch heißt es, das alle Werte mal 32768 (d.h. 2^15) skaliert sind. Die Tabelle berechnen wir durch:

TABELLE[i] = (2^15) * sin(i * PI/2 / TABELLENGRÖßE)

... für alle i von 0 bis 31. Dann fügen wir noch den Zusatzwert 32767 hinzu.

Die Genauigkeit unserer Tabelleneinträge würde uns nicht viel nützen, wenn wir nur nachschlagen würden. Unsere Eingabeauflösung wäre nur 1/32 (außer wir brauchen nur genau diese Winkel). Deshalb interpolieren wir zwischen zwei Tabellenwerten, wenn wir "Zwischenwinkel" als Eingabe erhalten, also:

Ergebnis = TABELLE[a] + b * (TABELLE[a+1] - TABELLE[a])

Dabei a der Winkel, den wir in der Tabelle nachschlagen, und b der Zwischenwinkel, also so etwas wie die Nachkommastellen.

Nicht vergessen - die Tabelle geht nur bis 90˚.

Winkel jenseits von 0˚ .. 90˚

Schauen wir uns die Sinuskurve mal an, so sehen wir wo wir unser Programm optimieren können. Erstmal steigt sie von 0˚-90˚ an. Dise Werte haben wir gespeichert. Dann fällt sie von 90˚-180˚ genau so ab wie sie angesiegen ist. Wir erhalten also einen Wert zwischen 90˚-180˚ indem wir unseren Winkel einfach wieder zurückrechnen, d.h. 90˚ - Winkel und dann nachschlagen. Zwischen 180˚ und 360˚ wiederholt sich dies, nur negativ. D.h. dann müssen wir die Ausgabe mal -1 rechnen. Für Winkel, die nicht im Bereich von 0˚-360˚ liegen sind wir in einer anderen Periode. Alles was wir tun müssen ist so lange 360˚ hinzuaddieren oder abziehen bis der Winkel zwischen 0˚ und 360˚ liegt (oder Modulo rechnen). Das ist eigentlich alles. Um im Programm Winkel von 0˚ bis 360˚ statt 0˚-90˚ darzustellen brauchen müssen unsere Zahlen 4 mal so groß sein wie die Tabelle, also 2 bit größer. Zufälligerweise sagen diese beiden Bits einzeln betrachtet schon aus was wir tun müssen. Dazu mal eine kleine "Logiktabelle":

  BIT:          6   5   4  3  2  1  0
                N   F   L  L  L  L  L
               --- --- ---------------
  0˚            0   0   0  0  0  0  0

  0˚   - 90˚    0   0   ?  ?  ?  ?  ?
  90˚  - 180˚   0   1   ?  ?  ?  ?  ?
  180˚ - 270˚   1   0   ?  ?  ?  ?  ?
  270˚ - 360˚   1   1   ?  ?  ?  ?  ?

  359.99        1   1   1  1  1  1  1

Den beiden Bits habe ich bereits Kürzel gegeben. Man sieht, dass das bit6 nur zwischen 180˚-360˚ gesetzt ist, wir können es also abfragen um zu erkennen ob der Sinus negativ ist (deshalb N). Das bit5 ist gesetzt wenn wir nicht vorwärts, sondern die Tabelle rückwärts lesen müssen, deshalb F wie "flip".

Interpolation

Jetzt haben wir schon 7 Bits mit denen wir eine ganze Kurve nachschlagen können. Unsere Eingabe ist allerdings Q15, und interpoliert haben wir auch noch nicht. Glücklicherweise haben wir noch ein paar Bits frei. Da alles mal 2^15 skaliert sieht eine Q15 Zahl so aus:

    BIT     15 |  14 |  13 |  12 |  11  |  10  | ... |   1   |   0
Q15 WERT   -/+ | 1/2 | 1/4 | 1/8 | 2^-4 | 2^-5 | ... | 2^-14 | 2^-15

  0˚ = 0    0     0     0     0     0       0    ...
 90˚ = 1/4  0     0     1     0     0       0    ...
180˚ = 1/2  0     1     0     0     0       0    ...
270˚ = 3/4  0     1     1     0     0       0    ...

Aus dieser Tabelle soll folgendes hervorgehen: Unsere bits für's Nachschlagen (inklusive den beiden Bits um bis 360˚ zu kommen) liegen nicht "unten" in der Zahl, sondern kommen direkt nach dem Komma-Bit. D.h. sie liegen eigentlich so in der Zahl:

 BIT     F   E   D   C   B   A   9   8   7   6   5   4   3   2   1   0
 -----------------------------------------------------------------------
 BIT     S | N | F | L   L   L   L   L | I   I   I   I   I   I   I   I
 ----------|-------|---|---|---|---|---|---|---|---|---|---|---|---|----
           |       |                   |
         S | N   F |     Lookup        |        Interpolieren
         i | e   l |                   |
         g | g   i |    (0-90deg)      |
         n | a   p |                   |
           | t     |                   |
           | e     |                   |

Und damit sehen wir auch mit welchen Werten wir interpolieren, es sind die I-Bits, deren Wertigkeit kleiner als die der Tabelle ist. Abzüglich dem Vorzeichenbit, den beiden "360 Grad"-Bits und der Tabelle bleiben noch 8 Bit übrig, die wir zur Interpolation verwenden.

Implementierung

Gießen wir das in eine Funktion wie unten in sin1.c angezeigt.

1. Makros

2. Unsere Eingabe ist Q15, d.h. es können auch negative Winkel "reinkommen". Die handeln wir damit schnell ab, dass wir 360˚ dazuaddieren. Allerdings wäre das 0x8000. Das kann ein Problem auf manchen Prozessoren sein, weil ein signed int16 nur bis 0x7fff geht. Also tun wir das sicherheitshalber in zwei Schritten. Wenn die Platform es kann wird der Compiler das wegoptimieren.

3. Jetzt schieben wir uns den Eingabewinkel zum Nachschlagen zurecht, also um die Anzahl an Interpolationsbits nach rechts.

4. Dann schauen wir nach, ob wir die Tabelle rückwärts lesen müssen. Ist das so, so rechnen wir 90˚ - Winkel. Hier können wir tricksen, den 90˚ sind 32, in unserer Tabelle, und 32 ist eine Dualzahl. Wenn ich die Bits alle umkippe, also das Zweierkomplement bilde, habe ich 32-Winkel schon gerechnet (Variable v0). Wir müssen das auch für die Interpolation tun, also mit den unteren Bits des Eingabewinkels. Das Ergebnis landet in v1.

5. Weil wir alle Bits in v0 gekippt haben müssen wir die Zahl zurechtstutzen, also alle Bits, die nicht zu unserer Tabelle gehören, maskieren (=0 setzen).

6. Jetzt schlagen wir mit v1 nach und nehmen auch den darauf folgenden Tabellenwert. Die Differenz ist quasi das Maximum des Interpolationswertes. Den Rechnen wir * v1 / 256, wobei /256 durch achtfaches Rechtsschieben gerechnet wird. Zu der Interpolation addieren wir noch den Tabellenwert selbst.

7. Zum Schluss prüfen wir noch, ob wir das Ergebnis negieren müssen weil der Winkel größer als 180 ist.

8. return