#Clocks für float division

[WarChild]

Hallo,

Ich habe folgendes Problem:
Ich habe leichten Zeitmangel bei der ausführung meiner Berechnungen.
Ich berechne viele Sinus, Cosinus, Arussinus, Akuskosinus, und einige Wurzeln mittles der AVR math.h bibliothek. außerdem muss ich gelegentlich von Rad in Grad umrechnen. wenn ich im simulator durchklicke, dann scheint es mir so, als würden für die berechnung Radian/M_PI*180 3000 Clocks und für die eigentliche Winkelfunktion nochmal 4000 Clocks drauf gehen. Um jetzt zeit zu sparen, dachte ich mir, dass ich einfach die umrechnung umgehe, indem ich eigene funktionen schreibe, die vlt nicht ganz so exakt sind wie die der Bibliothek, aber dafür noch schneller. Mehr als 0.05 Grad kann ich eh nicht auflösen.

Kann es also realistisch sein, dass eine float division 3000 Takte frisst? Wenn ja, dann sollte ich bei eigenen Algorithmen eine float division vermeiden. Wenn die nur 50 Takte frisst dann lohnt es kaum darüber nachzudenken.


vielen Dank für Euer Feedback!

mfg WarChild

[Besserwessi]

Ein paar Daten zur geschwingdigkeit von GCC findet man z.B. hier :
"avr-libc-user-manual.pdf" unter Benchmarks.
Sollte sich im doc-ordner von WINAVR befinden.

Danach sollte die Division etwa 466 Zyklen brauchen. Die 4000 Zyklen für die Winkelfunktion sind grob OK.

Man sollte beim Compilieren mit GCC-AVR darauf achten die das man die OPtion -Lm mit drin hat, damit die schnellere Mathe bibliothek genutzt wird. Das reduziert auch die Programmlänge.

Die Rechnung Winkel / Pi*180 sollte eigentlich im wesenlichen zur Compilezeit berechent werden. Notfalls als Winkel*(180/Pi) schreiben. Dann sollte zur Laufzeit nur eine Multiplication übrig bleiben.

[Felix G]

Kann es also realistisch sein, dass eine float division 3000 Takte frisst? Wenn ja, dann sollte ich bei eigenen Algorithmen eine float division vermeiden.Float Divisionen sollte man eigentlich grundsätzlich immer vermeiden, vor allem auf einem 8Bit µC wie den AVRs.

Besser noch man vermeidet Floats allgemein, und berechnet so viel wie möglich mit Festkommaarithmetik.


Auch kann man mit ein paar Tricks bei spezielleren Problemen oftmals massenhaft Zyklen sparen. Nehmen wir mal an, man möchte einen Float mit 2^n multiplizieren, dann geht das am schnellsten indem man eine Funktion schreibt, die auf den Exponenten einfach n aufaddiert. Das ist z.B. auch hervorragend geeignet, wenn man einen Float in einen Int konvertieren und dabei gleich hochskalieren möchte (für mehr Genauigkeit).

Das ist nur ein Beispiel wie man schnellen Code schreiben kann, wenn man ein Problem erstmal genau betrachtet.

[WarChild]

-lm ist im makefile schon drin (große hoffnung kaputt), das umklammern kann ich nochmal ergänzen. Kann ich auch z.B. x*(1.0/(180*M_PI)) berechnen lassen?
eine float multiplikation sollte doch eigentlich nur zwei zyklen benötigen oder nicht? s.Datenblatt FMUL Anweisung. Jedoch dauert das bei mir ebenfalls wesenlich länger. 500 Takte.

Und ich habe noch eine frage zum AVR float format. Für eine signum(float x) funktion benötige ich die anzahl der bits, die eine float zahl einnimmt, da das erste bit das vorzeichen enthält.
Der normale vergleich der float zahl > 0 oder nicht benötigt 400 Takte.
Ich habe schon versucht einen simplen int8_t pointer auf die float Zahl zu setzen und dann nach größer oder kleiner null zu prüfen, aber das schien nicht zu funktionieren. Auf was zeigt die Adresse einer float zahl? auf eine struktur oder einfach nur einen zusammenhängenden 4Byte speicherbereich?

mfg Warchild

[Besserwessi]

Die Multiplication sollte gut 130 Zyklen dauern. Der ASM-befehl FMUL ist für Festkommazahlen, nicht für Fleißkomma. Der wird da nicht viel weiter helfen.

Es könnte sein, das der Vergleich x>0 immer noch als Test von (x-0) gemacht wird. Eine Subtraktion ist dabei nur unwesentlich schneller als die Multiplikation. Das liegt daran, da man zur addition die beiden Zahlen erst mal auf den gleichen exponenten bringen muß. Beim verschieben der Bits ist des AVR nicht besonders effektiv.

Bei den Routinen für Fließkomma soll es irgenwann bei etwa GCC Version 4 eine deutliche Verbesserung gegeben haben. Es könnte sich lohnen eine neuere VErsion von WINAVR zu nutzen.



Das Format der Floats sollte etwa so aussehen:
3 bytes für die mantisse, wobei das führende Bit fehlt, denn es ist immer eine 1.
1 byte für Vorzeichen, Exponent und Speziallwerte wie Null, Nan.

Fürs separieren des Vorzeichens gibt es in math.h sogar eine extra Funktion:
int signbit (double __x) __ATTR_CONST__;
Sonst sollte das Bit 7 im 4 ten Byte für das Vorzeichen stehen. Also nicht das erste Bit.

Die addresse ein float Zahl sollte auf das erste der 4 bytes zeigen. Man könnte das auch als Strucktur sehen, oder einfach nur als 4 bytes.

[Felix G]

Auf was zeigt die Adresse einer float zahl? auf eine struktur oder einfach nur einen zusammenhängenden 4Byte speicherbereich?Im Normalfall auf einen 4Byte Speicherbereich, der eine IEEE 754 Single Precision Gleitkommazahl enthält.

Wie genau dieser Datentyp aufgebaut ist kannst du z.B. unter http://de.wikipedia.org/wiki/Gleitkommazahl nachlesen. Da steht unter anderem auch, daß es für ein und dieselbe Zahl mehrere mögliche Darstellungen geben kann (aber nicht muss), und was der Compiler alles beachten muss um zum Standard kompatibel zu sein. Spätestens wenn man sich das mal durchgelesen hat wird einem klar, wieso Float-Operationen so verdammt lange dauern.

Andererseits helfen einem diese Informationen aber eben auch, um für bestimmte Spezialfälle sehr schnelle Funktionen schreiben zu können (wie z.B. die von mir schon beschriebene Multiplikation mit 2^n).

[WarChild]

Ja, danke für eure tips.

Ich habe jetzt meinen Code bezüglich der Laufzeit optimiert und konnte folgende Ergebnisse erzielen:
Clocks
vorher / nachher Funktion
2500/325 Berechnung der benötigten Pulslängen (8 mal schneller)
34000/382 Berechnung der Translation (100 mal schneller)
175000/152000 Berechnung der Rotation (18% schneler)
400000/375000 Berechnung der Gelenkwinkel IK (8% schneller)

Alle anderen Funktionen benötigten bereits zuvor nur wenige hundert zyklen. Aber es bleibt die Problematik, dass die Rotation und die IK zu langsam sind. Ich werde die rotation, wenn meine Umrechnungen und meine eigene Sinus/Cosinus Funktion funktionieren wahrscheinlich ähnlich wie die Translation beschleunigen können, aber es bleibt das Problem der IK. Gibt es ein 32Bit Integer? ich habe das problem, dass meine zahlen, die im bereich -240 bis 240 liegen bei festkommadarstellung mit 16 Bit (7nachkomma Bits) beim quadrieren überlaufen, und meine funktionen nicht mehr funktionieren wenn ich vorm quadrieren selbst erst herabskaliere,dann können zahlen nahe der null einfach verschwinden.

mfg WarChild

[Besserwessi]

Natürlich gibt es 32 bit Integer Zahlen: long und unsigned long.

[Felix G]

Wenn du in C Integer Variablen mit definierter Länge brauchst (long kann bei der einen Architektur 32 Bit haben, bei einer anderen aber auch gern mal 16 oder 64), dann verwende die stdint.h oder inttypes.h (letztere bietet zusätzlich zu den neuen Typen noch ein paar mehr oder weniger nützliche Makros).

da drin sind Integer-Typen definiert wie uint8_t (= unsigned char) oder int32_t (= long int), bei denen man ganz genau weiss wie lang sie sind.

[WarChild]

Also Festpunktarithmetik ist schon so eine sache.
Auf dem PC funktionieren diese berechnungen, aber im µC nicht. Ich vermute da ist irgendwo ein Überlauf, den ich nicht sehe.
Ich habe alle Typen entsprechend umgewandelt.
hier meine Festpunkt sinus aproximation:


const int16_t W=180; //Winkel eines Halbkreises
const int16_t W2=360; //Winkel eines Vollkreises
const int16_t Wd2=90; //Winkel eines Viertelkreises

int16_t sinus(int16_t x) //Eingabe mit Skalierung 128 = 7 Festkommabits Ausgabe 14 festkommabits
{

if((x>>7)>W) x=-x+(W<<7); //x wird auf das Intervall [-W;W] normiert (aufgrund der Achsensymmetrie um x=0)
if((x>>7)<-W) x=-x-(W<<7);

//Parabel
const int16_t B = 182; //2^-13 //(4/W); //linearer Formfaktor der Parabel
const int16_t C = -259;//2^-21 //(-4/(W*W)); //quadratischer Formfaktor der Parabel

long y=((B*x)>>6)+((((C*x)>>11)*x*signi(x))>>10); //2^-14 //Funktionswert der einfachen Parabel

//Parabel Korrektur
const int16_t Q = 99; //2^-7 //0.775; //Linearfaktor der einfachen parabel
const int16_t P = 29; //2^-7 //0.225; //Linearfaktor der quadrierten korektur Parabel

y=((Q*y)>>7)+((((P*y)>>7)*y*signl(y))>>14); //2^-14 //Endergebnis nach gewichteter Summenbildung

return y;
}


mfg Warchild

[WarChild]

Hier noch ergänzend meine Signum funktionen.


inline int8_t signf(float number)
{
if (signbit(number))
return -1;
else
return 1;
}
inline int8_t signi(int16_t number)
{
if(number>>15)
return -1;
else
return 1;
}
inline int8_t signl(long number)
{
if(number>>31)
return -1;
else
return 1;
}


nicht, dass dort der fehler herkommt

[Felix G]

Erstmal hätte ich ein paar kleine Hinweise zur Performance:


if((x>>7)>W) x=-x+(W<<7); //x wird auf das Intervall [-W;W] normiert (aufgrund der Achsensymmetrie um x=0)
if((x>>7)<-W) x=-x-(W<<7);Die Shifterei kannst du dir sparen, wenn du die Konstanten für 90°, 180° und 360° gleich in skalierter Form definierst. Das spart immerhin allein bei diesen 2 Zeilen schon 28 Shift-Operationen (der AVR muss ja für ein <<7 tatsächlich 7x shiften)

Auch die Signum Funtionen lassen sich beschleunigen, indem du statt if(number >> 31) folgendes machst: if(number & (1 << 31)). Dafür braucht der AVR keine Shifts, da 1 << 31 als Konstante schon vom Compiler berechnet wird.


Was den Überlauf.. bzw. die Überläufe betrifft sind mir mehrere Stellen aufgefallen:

1. x hat den Typ int16_t und enthält den Winkel in Grad * 128, korrekt?
Dann möchte ich mal anmerken, daß ein Winkel von 360° einem x von 46080 entspricht. Bei einem int16_t ist allerdings schon bei (2^15)-1 = 32767 Schluss.

2. B und x sind beide int16_t und werden multipliziert, da ist es sicherer zuerst beide auf uint32_t zu casten, damit der Compiler nicht auf die Idee kommt das Ergebnis in 16 Bit quetschen zu wollen.

Es sind noch mehr Probleme deser Art im Code zu finden, daher mein Tipp:

Für jeden einzelnen Berechnungsschritt nochmal den worst-case Fall prüfen. Prinzipiell gilt: wenn irgendein Zwischenergebnis größer werden könnte als der Typ der beiden Operatoren, sollten beide vorher zur Sicherheit auf den nächst größeren Typ gecastet werden.

[WarChild]

ok vielen dank,
zu 1. meine funktion soll nur für die Rotation meiner Koordinaten verwendet werden. Den Torso kann ich nur um max 60° drehen um Z sogar nur um 15°, sonst kippt er um, daher war ich von maximal 180° ausgegangen, und dafür genügt ein int16_t mit 7 Festkommabits. Dadurch habe ich das leiche format wie bei meinen Drehwinkeln. Also bevor die funktion einen Überlauf hat, sind bereits meine Drehwinkel übergelaufen.

Aber mit 2. hast du recht. Kann sehr gut sein, dass auf meinem 32Bit rechner die überläufe nicht auffielen, da er ehh mit "long" für AVRs arbeitet.

mfg WarChild

[WarChild]

So: Felix, du hattest recht
jetzt habe ich die minimale Anzahl an Float operationen erreicht, es bleiben nur je sechs mal zwei Wurzeln, vier Winkelfunktionen und vier Divisionen bei der einen und sechs Winkelfunktionen bei der anderen Funktion, die ich nur durch eigene Festkommafunktionen beschleunigen könnte.

vorher / nachher Funktion
2500/325 Berechnung der benötigten Pulslängen (8 mal schneller)
34000/382 Berechnung der Translation (100 mal schneller)
175000/51000 Berechnung der Rotation (3.4 mal schneler)
400000/276000 Berechnung der Gelenkwinkel IK (45% schneller)
jetzt liege ich bei eine gesamten Rechndauer mit Interrups bei 22.5ms,
also fast echtzeit. Eigentlich muss mich nurnoch die Beschleunigung von der IK Interessieren, da diese mit 17,5ms den größten Ausschlag macht.

vielen Dank an Alle, ohne euch wäre ich vorher warscheinlich schon verzweifelt!!! =D>
edit: meine Festpunkt Sinus/Cosinus Funktion mit max 0.0001 Abweichung funktioniert jetzt. Also:
175000/11100 Berechnung der Rotation (15 mal schneler) \:D/
Fehlt nurnoch eine Festpunkt Wurzel und Arkus-funktionen

mfg WarChild

[WarChild]

So meine Quadratwurzel ist fertig.
Hab mir ein Iteratives VErfahren ausgedacht, dass nach dem Prinzip der Annäherung durch Intervallhalbierung funktioniert.
dadurch konnte ich folgende verbesserung erziehlen:
400000/232000 Berechnung der Gelenkwinkel IK (72% schneller)

hat jemand von euch eine Idee wie man die Arkus Winkelfunktionen Approximieren kann? In Festkomma Arithmetik...


uint16_t sqroot(unsigned long number) // O(n)=log_2(n) Max Abweichung meistens 0 oder 1, selten 2 bzw. 0.1% für große zahlen
{
if(number==0) return 0; //triviale Fälle abfangen
if(number==1) return 1;

unsigned long r2;
uint16_t root=0,delta;
if(number>16777216UL) //falls die Zahl größer als eine 24 Bit Zahl ist
delta=(number>>16);
else if(number>64536) //falls die Zahl größer als eine 16 Bit Zahl ist
delta=(number>>8); //günstigeren Anfangspunkt Wählen --> führt zu leichter Unschärfe 0.1%
else
delta=(number>>1); //sonst beginn bei der Zahl selbst

while (delta) //sollte innerhalb von 16 Iterationen eine Lösung liefern
{
r2=(unsigned long)root*root;
if(r2>number)
root-=delta;
else
root+=delta;
delta=(delta>>1);
}
return (uint16_t)root;
}


mfg WarChild

[Besserwessi]

Im Prinzip sollten die Potenzreihen auch in Festkommaarithmetik gehen. Die Reihe vom Arctan ist aber bekannt dafür, dass sie nicht besonders schnell konvertiert.

Ein guter Link dazu :
http://de.wikipedia.org/wiki/Arctan#N.C3.A4herungsweise_Berechnung

[WarChild]

vielen Dank für den Tipp
Der Arkustangens ist schon so gut wie erledigt.
Das Problem bei den letzten beiden (Arkus-Co-Sinus) ist die abweichung in den Randbereichen. Selbst Taylorpolynome 7. Grades haben in der nähe um die -1 und 1 Abweichungen von bis zu 4 Grad. Mal sehen, vlt. nehm aber rotzdem ein Taylorpolynom, da meine mechanischen Gelenke abgesehen von der schulter nur in den Intervallen [5°;175°] und [10°;180°] arbeiten können. So wird eine Überstrapazierung der Getriebe durch das Erreichen des Endandschlages von vorn herein vermieden.

[Besserwessi]

Beim ArcSin / cos macht es eventuell Sinn, den Wertebereich aufzuteilen. Bei Werten von etwa -0.5 ... 0.5 mit dem Arcussinus, von 0.5 ... 1 mit dem Arcus cosinus.
Der Arccos um 1 ist etwas problematisch, schon wegen der großen Steigung. Eventuell wäre hier eine andere Näherung, z.B. über die Wurzel sinnvoller. Für Werte nahe 1 sollte sqrt(2-2*x)) ein brauchbare Näherung für Arccos(x) sein.

[WarChild]

Danke Beserwessi:
Dank der Festkomma ArcTan(x) und ArcTan2(y,x) funktionen konnte ich wieder eine erhebliche Beschleunigung erziehlen.
400000/149000 Berechnung der Gelenkwinkel IK (2.7 mal schneller)
Aufgrund dieses Leistungssprungs schafft der µC alle nötigen berechnungen in nur 10ms zu erledigen. Also noch reichlich Zeit für die PWMs vlt. eine dynamische Fußnachsetz Funktion und I²C kommunikation.
Allerdings reizt es mich die floats gänzlich aus meinem Bot zu verbannen.
Vlt klapts ja noch mit den letzten beiden Funktionen, auch wenn mich die Überlaufkontrolle mitlerweile ankotzt.


//**************************************************// 974 Takte*
int16_t arctan(long number) //Eingabe mit Skalierung 1024 = 10 Festkommabits Ausgabe 7 Festkommabits
{
int8_t v;
unsigned long x;
if (number==0) return 0;
if (number>0) {x= number;v=1;} // Betrag und Vorzeichen von number
else {x=-number;v=-1;}
if (x> 8388608UL) return (v* Wd2); //für |x|>2^13 gilt +-90°
if (x<1024) //für x<1 gilt allg.: x/(1+0.28x)
return v*(((209539UL*x)>>3)/(((x*x)>>10)+3657)); //zähler/nenner= 2^-17/2^-10=2^-7
else if(x<60000UL) // sonst gilt allg.: PI/2-x/(x²+0.28)
return v*(Wd2-((58671UL*(x>>3))/(((x*x)>>10)+287))); //zähler/nenner= 2^-17/2^-10=2^-7
else
{
x=(x>>5); //anpassung zum Schutz vor Überläufen
return v*(Wd2-(((58671UL*x)>>8)/(((x*x)+287)>>10))); //zähler/nenner= 2^-7/2^0=2^-7
}
}
//**************************************************//

//**************************************************// 1729 Takte*
int16_t arctan2(int16_t y,int16_t x) //Eingaben Ganzzahlig, Ausgabe 7 Festkommabits
{
if(x==0) //Vorzeichenorientierte Fallunterscheidung
{
if (y==0) return 0;
else return signi(y)*Wd2;
}
if(x>0) return arctan(((long)y<<10)/x);
else return arctan(((long)y<<10)/x)+signi(y)*W;
}
//**************************************************//



mfg WarChild

[WarChild]

@Besserwessi
Ich hab mir mal die kurve der Arcos Funktion angesehen und der Tipp mit der Wurzel ist gar nicht schlecht. Ist ja auch logisch, wenn die Funktion einer Parabel Ähnelt, so muss die Umkehrfunktion einer Quadratwurzel ähneln.
Hier also meine Formel zur Aproximation: -1<x<1
für x==0 -->90 //Optionale abkürzungen, nicht notwendig
für x==1 -->0
für x==-1 -->180
für x>=0 --> 0.822*90*sqrt(1-x)-0.178*90*(1-x)
für x<0 --> 180-(0.822*90*sqrt(1+x)-0178*90*(1+x))

Ich denke mit meiner Ganzzahligen Wurzel, als Festkomma Interpretiert sollte das recht schnell funktionieren.
Diese Funktion hat eine maximale Absolute Abweichung bei 0.9 von 0.845°, die zweitgrößte Abweichung liegt bei 0.2 mit 0.6° aber an den Kritischen stellen stimmt Sie absolut.
Ich habe die Linearfaktoren der Wurzelfunktion und der Korrekturgeraden auf das Minimale Integral der absoluten Abweichung justiert.

naja und der Asin hat halt nur etwas andere Y-Verschiebung und eine Spiegelung. asin(x)=-acos(x)+90

edit: Funktion ist implementiert:
400000/78000 Berechnung der Gelenkwinkel (5.2 mal schneller)
jetzt schaffe ich die gesamte IK und Transformation in nur 5.8ms. Bäm
Keine Floats sind mehr vorhanden! Alles zusammen hat eine Verzehnfachung der Rechengeschwindigkeit bewirkt.

mfg WarChild