mikrocontroller.net

Forum: Mikrocontroller und Digitale Elektronik Code für Logarithmus gesucht


Autor: Jörg H. (idc-dragon)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich suche einen AVR-tauglichen (Festkomma-)Algorithmus, der den 
Logarithmus einer 32-Bit Zahl bildet. Zu welcher Basis ist egal, das ist 
ja nur ein Faktor im Ergebnis. Vermutlich ist die Basis 2 am 
geeignetsten für Bitschiebereien. Meine Eingangswerte liegen im Bereich 
von ca. 1000 bis 100.000.000.
Es muß nicht so fürchterlich genau sein, aber nur das höchste gesetzte 
Bit auszählen reicht mir nicht.

Hat da jemand vielleicht was Nettes parat? Vielen Dank!

Autor: Der Dude (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Potenzreihenentwicklung. Frage Wikipedia.

Autor: sechseinssechs (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Der Logarithmus ist relativ einfach. Es beginnt damit zu zaehlen auf 
welcher Position das MSB ist. Fuer die folgenden bits nimmt man ein 
Tabelle. Eine Potenzreihe ist unpassend.

Autor: sechseinssechs (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert

Autor: Stefan (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> Es muß nicht so fürchterlich genau sein, aber nur das höchste gesetzte
> Bit auszählen reicht mir nicht.

Ein genaueres Resultat kann man aber mit Festkomma-Zahlen gar nicht 
speichern. Das höchste Bit zu suchen ist also der einfachste, schnellste 
und genaueste Algorithmus der möglich ist mit Fix-Point-Arithmetik.

Beweisskizze:
x sei ein unsigned 32-bit Integer (ganz normal, kein Komma). Der 
Wertebereich liegt zwischen 0 und 2^32-1. log2(x) wird sich entsprechend 
zwischen -unendlich und 32 bewegen. Soll der log2 wiederum in dem selben 
Datentyp abgelegt werden, muss man auf ganze Zahlen runden (z.B. hier 
z.B. abrunden). Somit sind nur 32 Werte als Resultat möglich. Wenn du 
die wieder zurück-exponenzierst wirst du feststellen das dies genau jene 
32 Fälle sind, die du per MSB Bit-Zählen bekommst.

Nun zur Fix-Point-Arithmetik. Sei b=Komma-Position dann kann man eine 
32-bit-Fix-Punkt Zahl y Schreiben als y = x / 2^b. Wobei x ein 
32-bit-Int ist. Es folgt: log2(y) = log2(x / 2^b) = log2(x) - log2(2^b) 
= log2(x) - b. Sprich der fix-point-log unterscheidet sich nur duch eine 
ganzzahlige Subtraktion vom integer-log, eine höhere Genauigkeit ist 
also nicht möglich.

Wenn dir die Genauigkeit nicht ausreicht, müssen input und output 
Datentyp des Logarithmus eine unterschiedliche Komma-Position haben.

Autor: Stefan (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hmm Moment, irgenwas stimmt da doch ganz und gar nicht. Ich glaub was 
ich da geschrieben hab ist falsch, zumindest der fix-point Teil. Ich 
schau mir das Morgen nochmals an (brauch wohl etwas Schlaf) :-)

Autor: yalu (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@Stefan:

Endlich mal einer, der seine Aussagen ordentlich begründet. Aber ich
glaube, in dem Beweis steckt ein kleiner Fehler:

> log2(y) = log2(x / 2^b) = log2(x) - log2(2^b) = log2(x) - b

= (log2(x) - b) * 2^b / 2^b

(log2(x) - b) * 2^b (gerundet) ist also der Wert, der als Ergebnis in
die 32-Bit-Variable geschrieben wird und repräsentiert eine Zahl mit b
Nachkommastellen. Anders als in deiner ersten Beweishälfte ist log2(x)
jetzt natürlich keine Ganzzahl mehr.

Folglich ist es kein Problem, aus einer Festkommazahl mit b Nachkomma-
stellen den Logarithmus mit ebenfalls b Nachkommastellen zu berechnen.
Oder denke ich da falsch?

> Wenn dir die Genauigkeit nicht ausreicht, müssen input und output
> Datentyp des Logarithmus eine unterschiedliche Komma-Position haben.

Was sicher auch kein Problem wäre.

Mein Favorit wäre übrigens der erste Vorschlag von sechseinssechs
(Schieben und Tabelle).

Autor: Jörg H. (idc-dragon)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ich habe mich an den Vorschlag Schieben+Tabelle gehalten. Per Tabelle 
nähere ich aus 16 Geradenstücken an. Meine Annäherung liegt wegen der 
Krümmung unter der Originalkurve, für kleineren durchschnittlichen 
Fehler könnte man die etwas höher schieben, ich weiß aber nicht wie man 
diesen Ausgleich berechnen täte. Die Genauigkeit ist schon bemerkenswert 
hoch, nahe 10e-5.

Hier ist mein erster Code:
uint32_t log2fix(uint32_t n)
{
    static const uint16_t lut[17] = { // log2(1.0)...log2(2.0) in 1/16 interval
        0, 5732, 11136, 16248, 21098, 25711, 30109, 34312, 38336, 42196, 45904, 49472, 52911, 56229, 59434, 62534, 65535
    };
    uint8_t ri = 31; // result integer part
    uint16_t rf; // result fractional part

    uint8_t j; // table index
    uint16_t nf; // input value, fractional part

    while (!(n & 0x80000000)) // count bits to MSB
    {
        n <<= 1;
        ri--;
    }

    n <<= 1; // shift out the mantissa '1', now only fractional part left

    j = (uint8_t)(n >> (32-4)); // use 4 bit for table entry selection
    nf = (uint16_t)(n >> (16-4)); // use the bits below for linear interpolation

    rf = lut[j] + ((nf * (lut[j+1] - lut[j])) >> 16);
    
    return ((uint32_t)ri << 16) + rf; // return as 8.16 fixpoint
}
Die Hi/Lo Teile kann man mit Unions wohl noch optimieren.

Autor: Stefan (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Merci yalu
Genau da lag mein Fehler. Offensichtlich kann man mit Fix-Point 
Arithmetik den Logarithmus genauer darstellen. Das ist mir beim erneuten 
durchlesen dann auch aufgefallen, nur um mit dem Finger darauf zu zeigen 
war ich dann doch zu müde.
Sorry wenn ich für Verwirrung gesorgt habe :)

Schliesse mich 616 an. Die Zahl ins 1..2 Intervall schieben und dann 
eine der üblichen Approximationsmethoden verwenden. Der Logartihmus ist 
in dem Bereich ja ziemlich "zahm". Für eine erste grobe Approximation 
könnte man einfach den Wert übernehmen (Approx durch Gerade mit Steigung 
1).

Autor: Jörg H. (idc-dragon)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Mit Unions zum partiellen Zugriff auf die Teile eines 32-bit Wertes 
kriege ich es effizenter, im Hinblick auf 8 Bit Controller ohne 
Barrel-Shifter.
Außer der MSB-Suche sind keine Shifts mehr nötig, wenn ich damit passend 
für eine Bytegrenze aufhöre. Mit diesem Code kann ich die oberen 3 Bit 
des Eingangswertes nicht nutzen, aber so große Werte treten bei mir 
nicht auf. (Bei Verdopplung der Tabellengröße würde man ein Bit davon 
zurückgewinnen)
#pragma pack(1) // for MSVC
// gcc needs __attribute__((packed))

// this struct allows efficient access to the parts of a 32 bit value, without shifting/masking
typedef union
{   // union for little endian CPU, make sure this gets packed
    struct 
    {
        uint8_t ll8;
        uint8_t lh8;
        uint8_t hl8;
        uint8_t hh8;
    };
    struct
    {
        uint16_t l16;
        uint16_t h16;
    };
    struct
    {
        uint8_t l8;
        uint16_t m16; // this needs unaligned read capability
        uint8_t h8;
    };
    uint32_t val32;
} variant_t;


// logarithm to base 2, returns a 16.16 fixed point value
// restrictions: input value must not be 0 and not exceed 0x1FFFFFFF (because of an implementation shortcut)
uint32_t log2fix(uint32_t n)
{
    static const uint16_t lut[17] = { // log2(1.0)...log2(2.0) in 1/16 interval
        0, 5732, 11136, 16248, 21098, 25711, 30109, 34312, 38336, 42196, 45904, 49472, 52911, 56229, 59434, 62534, 65535
    };
    variant_t v; // alias for the input value
    variant_t r; // result
    uint8_t j; // table index

    r.hh8 = 0; // upper integer part is always zero
    r.hl8 = 28; // integer part is a bit countdown

    v.val32 = n;

    while (!(v.hh8 & 0xF0)) // count bits until MSB reaches place 28
    {
        v.val32 <<= 1;
        r.hl8--;
    }

    j = v.h8 & 0x0F; // use 4 bit in places 27...24 for table entry selection
    
    // use the bits below (in v.m16) for linear interpolation
    r.l16 = lut[j] + ((v.m16 * (lut[j+1] - lut[j])) >> 16);

    return r.val32; // return as 16.16 fixpoint
}

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.