Forum: Compiler & IDEs sinus berechnen


von Gast (Gast)


Lesenswert?

hi leute,

kurzer umriss:
ich suche nach einer alternative zur gnu bibliothek... ich brauche eine 
implementierung, die schneller rechnet... und auf die 4te 
nachkommastelle noch genau ist...

für eine tabelle ist nicht genug platz im speicher...

...daher meine frage, kennt jemand einen schnellen algorhitmus oder gar 
eine bibliothek die ich verwenden könnte??

(GNU-C)

vielen dank mal :)

von Tom E. (tkon)


Lesenswert?

hallo,

also ich denk mal wenn'st den vollen Sinus (0-45 Grad wegen Symmetrie) 
brauchst, wirst kaum eine Chance haben etwas merklich schnelleres zu 
finden da es dafür etablierte Standardverfahren gibt die praktisch 
überall eingesetzt werden und sehr schnell konvergieren.

Wie sieht es mit Speicherplatz für eine Tabelle aus? Vielleicht wär es 
möglich die zu komprimieren wenn es nur ein wenig an Platz fehlt

gr

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

Gast wrote:

> implementierung, die schneller rechnet... und auf die 4te
> nachkommastelle noch genau ist...

Binär oder Dezimal? **g**

> für eine tabelle ist nicht genug platz im speicher...

Algorithen bekommst du auch nicht für lau.

> ...daher meine frage, kennt jemand einen schnellen algorhitmus oder gar
> eine bibliothek die ich verwenden könnte??

Wo du nie dran vorbei kommst ist die Normierung auf einen kleinen 
Argumentbereich, also erstmal mod 2pi, mod pi, mod pi/2 oder mod pi/4 
rechnen.

Für den sin selber gibt's 3 Ideen:

1) Taylor um 0
2) Taylor um pi/4 oder pi/8 (dann auch den cos um pi/8)
3) Approximation durch ein Polynom, z.b. 4. Grades.

Ansatz ist

Das Polynom wird so gewählt, daß
was also eine Optimierungsaufgabe für die 5 Koeffizienten von p 
darstellt. Da es sich um ein Polynom handelt, ist es unangenehm, für 
Argumente die betragsmäßig größer als 1 sind, auswerten zu müssen. Soll 
der Bereuch daher auf 0..pi/4 beschränkt werden, kommt noch eine analoge 
Aufgabe für den Cosinus hinzu bzw. man zwelegt das Intervall der Länge 
pi/2 in 2 der Länge pi/4 und nimmt 2 Polynome, eines für jedes 
Teilintervall.

Als Maß bietet sich an
oder

Das ganze dürfte auf ne Taylor-Entwicklung rauslaufen, deren 
Entwicklungspunkt etwa in der Intervallmitte ist.

Allerdings ist die Funktion per Konstruktion wegen der Minimalität der 
Norm optimal geeignet für den Zweck.

Für die Arithmetik:

Wie wär's mit Q-Format in Assembler? Also erst mal checken, welches 
Q-Format passen würde. 1.15 ist zwar recht easy zu implementieren, 
allerdings bleiben die (Zwischen)Werte nicht im Invervall (-1,1). Genehm 
wäre ein 3.29 Q-Format, das immerhin (-4,4) darstellen kann. Bei 3.13 
reichen die Nachkommastellen wohl nicht für die Genauigkeit aus. 
Speziell für den Sinus hatte ich ein 16-er Q-Format mal implementiert, 
aber wo...?

von Gast (Gast)


Lesenswert?

da kommen wir doch schon näher :)

Q-format ist gut... weiss nicht wie sinnvoll, aber dachte an 1.15

aber implementiert die libc des gcc nicht eh schon die 
taylor-approximation?

ich dachte an CORDIC in Q

aber ich raff den CORDIC nich, bin zu doof :(

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

Gast wrote:

> aber implementiert die libc des gcc nicht eh schon die
> taylor-approximation?

Doch. Aber deren Entwicklungspunkt ist 0, wie gesagt.

Stell dir einfach vor, du willst den sin durch ein Polynom 1. Grades 
(also eine Gerade) annähern. Bei Taylor mit Entwicklungspunkt 0 ist das 
p(x)= x. Wenn du die Gerade so legen kannst, daß zB der über alle Werte 
genommene Fehler (in der gewählten Norm) minimal ist, bekommst du sowas 
wie ne Bestgerade an sin in [0, pi/2]. Das approximiert besser und ist 
aber auch nur ein Polynom ersten Grades.

Daran sieht man (entgegen meiner obigen Vermutung) daß diese 
Approximation keine Taylor-Entwicklung ist, weil sie keine Tangente an 
sin liefert.

Oben ist das Verfahren eben für ein Polynom höheren Grades gemacht. Der 
sin hat den Vorteil, daß die geraden Terme wegfallen. Die Approximation 
hat den Vorteil, daß sie sich nicht auf den Punkt 0 fixiert und mit 
roher Gewalt ;-) so lange rechnet, bis der Fehler/Restterm klein genug 
ist, sondern der Fehler per Konstruktion klein ist. Ob das im Endeffekt 
besser/schneller ist, muss man sich anschauen.

Zudem implementiert die libc in float und nicht in Q-Format, was 
aufwändig ist. Des weiteren muss sie sich nicht im NaNs etc kümmern und 
die floats aus- und einpacken und ist nicht in asm geschrieben.

von Gast (Gast)


Lesenswert?

hmmmm....

in welche richtung empfielt sich nun zu gehen?

Q-format steht fest, aber soll ich in Q nun taylorreihen verfolgen oder 
soll ich mich weiter in CORDIC einlesen (hier soll der geschwindigkeits 
vorteil in der ausschliesslichen verwendung von additionen und shifts 
liegen)...

ich bewege mich nur im bereich 0..PI/2

danke soweit mal für die infos

von dvh (Gast)


Lesenswert?

Wenn man Stützstellen bei 30 und 60 Grad (zwischen 0 und 90) nimmt, dann 
ist dies eine Naeherung  auf besser als 1.0E-6:
wenn z.B. x=75.0*pi/180.0 dann ist 60° zu verwenden
mit dx=x-x60 (x60=60.0*pi/180.0) gilt dann:
sin(x)=sin(x60)*cos(dx)+cos(x60)*sin(dx)
mit sin(x60)=sqrt(3)/2 und cos(x60)=0.5.
Für die kleinen Werte dx gilt in der 1.0E-6 Naeherung
sin(dx)=dx*(1-1/6*dx^2*(1-1/20*dx^2)) und
cos(dx)=1-1/2*dx^2*(1-1/12*dx^2)
(Taylorentwicklung)
ich bekomme zB. bei 75°
s60*(1.0-dx**2*0.5*(1.0-dx**2/12.0))+c60*( 
dx*(1.0-1.0/6.0*dx**2*(1.0-0.05*dx**2)))-sin(x)=3.9514708150001354e-07

von Alex (Gast)


Lesenswert?

Hier mal die Implementierung aus einer von mir genutzten fixed point 
(16.16 Format) Klasse.
1
    Fxp32 Fxp32::sineInternal(Fxp32 Value)
2
    {
3
        // Polynomial length is chosen such that theoretical error < 2*10^(-4)
4
        // Matlab code:
5
        // x  = 0:.001:pi/2;
6
        // y  = sin(x);
7
        // xn = x ./ pi; %normalize x from (0, pi/2) to (0, .5)
8
        // p = polyfit(xn, y, 4);
9
        // p = p * 2^16; % convert to 16.16 format
10
        const UInt8 PolyLength = 5;
11
        const Int32 PolyCoeff[PolyLength] = { 0x0000000FL,  //  15.137052889215985
12
                                              0x0003209AL,  //   2,049535691178602e+05
13
                                              0x000034CAL,  //   1.351355002615320e+04
14
                                              0xFFF9AE1CL,  //  -4.141804497187611e+05                                             
15
                                              0x0002CA0FL   //   1.827989287236855e+05
16
                                            };
17
18
        Fxp32 Temp;
19
        Fxp32 Result;
20
        Fxp32 Exp;
21
        Int8  i;
22
        bool IsNeg = false;
23
24
        // Use identity sin(x) = -sin(-x) and make make Value positive, remember sign
25
        if (Value < Fxp32(0))
26
        {
27
            Value = -Value;
28
            IsNeg = !IsNeg;
29
        }
30
31
        // Use identity sin(x) = -sin(-x) to map 3rd and 4th quadrant to 1st and 2nd
32
        if (Value.m_Value & ConversionFactor)
33
        {
34
            IsNeg = !IsNeg;
35
        }
36
37
        // Normalize to range (0, 1]
38
        Value.m_Value &= 0x0000FFFFL;
39
40
        // Now m_Value maps to the first positive half wave
41
        // Check if 1st or 2nd quadrant, if 2nd use identity sin(pi/2+x) = sin(pi/2-x)
42
        if (Value.m_Value & HalfIncrement)
43
        {
44
            // This means value >= 0.5, i.e. 2nd quadrant => map to 1st
45
            Value.m_Value = (ConversionFactor - Value.m_Value);
46
        }
47
48
        // Polynom evaluation in this way:
49
        // Result = X*(X*(A*X+B)+C)+D
50
        // (requires less multiplications)
51
        // Constant term
52
        Result.m_Value = PolyCoeff[PolyLength-1];
53
        Result        *= Value;
54
        Temp.m_Value   = PolyCoeff[PolyLength-2];
55
        Result        += Temp;
56
57
        // Use polynom to approximate sine
58
        for (i=PolyLength-3; i>=0; i--)
59
        {
60
            Temp.m_Value = PolyCoeff[i];
61
            Result      *= Value;
62
            Result      += Temp;
63
        }
64
65
        if (IsNeg)
66
        {
67
            Result = -Result;
68
        }
69
70
        return Result;
71
    }
72
73
    Fxp32 Fxp32::sine(Fxp32 Value)
74
    {
75
        Fxp32 Temp;
76
77
        // Normalize Value by dividing by PI (map interval +/-Pi to +/-1)
78
        Temp.m_Value = PiInverse;
79
        Value       *= Temp;
80
81
        return sineInternal(Value);
82
    }

von Werner Grimmel (Gast)


Lesenswert?

Ohne alles obige im Detail gelesen zu haben:

Google mal nach "Bezierkurven" oder "Beziersplines". Das ist nichts 
anderes als ein Polynom (Bernsteinpolynome etc...). Du betrachtest den 
Sinus als Graph in einem kartesischen Koordinatensystem und rechnest 
also die x/y Koordinaten aus für alle Werte bis Pi/2. Der Rest ist ja 
symmetrisch.

Du musst für einen Bezierspline einen Startpunkt und einen Endpunkt 
festlegen
(also z.B. (0,0) und (Pi/2,1)). Jetzt benötigst du Stützstellen. Die 
Stützstellen liegen hier aber NICHT auf der Sinuskurve, sondern sind 
Punkte im Koordinatensystem die den Verlauf des Graphen beeinflussen.

Du nimmst einen Algorithmus zur Berechnung des Splines auf Basis von 
Additionen und Shifts und fertig. Keine explizite Multiplikation nötig.

Ist sehr sehr einfach, ganz wenig Code und recht schnell da keine 
Multiplikation nötig. Genauigkeit theoretisch beliebig erweiterbar durch 
Umsaklieren des Koordinatensystems oder Einführung zusätzlicher 
Stützstellen.

google Stichwörter:

- Satz von de Casteljau
- Fadenkonstruktion Parabel
- Bernsteinpolynome
- Bezierspline

--> Mal eines der zahlreichen Java-Applets zu den Splines im Internet 
ansehen! Es ist fantastisch was man mit diesen Kurven alles machen 
kann!! Und superleicht!!

von Gast (Gast)


Lesenswert?

das mit dem splining ist interessant, aber wenn ich dann meinen punkt 
auf dem sinusähnliche gebilde habe, bin ich noch immer nicht bei dem 
gesuchten winkel... (asin())

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

Gast wrote:
> das mit dem splining ist interessant, aber

Auch nur Polynome, die es auszuwerten gilt.

> wenn ich dann meinen punkt
> auf dem sinusähnliche gebilde habe, bin ich noch immer nicht bei dem
> gesuchten winkel... (asin())

Wenn du sin willst, dann berechne sin.
Wenn du asin willst, dann berechne asin.

von Gast (Gast)


Lesenswert?

>Wenn du sin willst, dann berechne sin.
>Wenn du asin willst, dann berechne asin.

der weg ist mein ziel ;)

ich danke allen für eure infos, ist alles sehr wertvoll für mich... 
danke :)

von Hans-jürgen H. (hjherbert) Benutzerseite


Angehängte Dateien:

Lesenswert?

Ein Polynom 4. Ordnung sollte reichen fuer Genauigkeit 200 ppm

Und wenns schneller gehen soll siehst Du da
Beitrag Beitrag "float arithmetic"

(Jene Programme verzichen auf Behandlung von Sonderfällen)

von Dennis (Gast)


Lesenswert?

Blöd gefragt:

Externen Sinusgenerator auf den ADC anschliessen und Meßwert 
verarbeiten...

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

Ist das schneller? Wie lange muss man warten, wenn man sin(1000) will?

von Dennis (Gast)


Lesenswert?

sin(1000) = sin(280), höher als 360 macht es eh kein Sinn...

Ob es schneller ist? Kommt auf die Frequenz und Anwendung an.

Vorteil: SW in 10 Min geschrieben, ext. Sinusgenerator kann evtl. auch 
getriggert werden (Anwendungsabhängig). Wie im vorigen Beitrag gesagt: 
blöde Idee, aber es kann evtl. auch funktionieren. Wurde uns im Studium 
nicht jeden Tag eingebläut, auch mal "unübliche" Wege zu gehen, wenn sie 
den einfacher/besser/störfester/wasauchimmer sind?

Gruss:

Dennis

von Gast (Gast)


Lesenswert?

danke für die anregung....

ich suche etwas extrem schnelles

Je Sekunde muss ich 20000 mal je einen sin() einen cos(mit anderem 
winkel) und einen arccos ausfürhen...

(un jede menge anderes)

denke ich muss dann doch 4kb flash opfern für eine tabelle...

danke euch allen soweit :)

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.