Gleitkomma-Linearinterpolation

Question 1

Um eine lineare Interpolation zwischen zwei Variablen durchzuführen a und b einen Bruchteil gegeben fverwende ich derzeit diesen Code:

float lerp(float a, float b, float f) 
{
    return (a * (1.0 - f)) + (b * f);
}

Ich denke, es gibt wahrscheinlich einen effizienteren Weg, dies zu tun. Ich verwende einen Mikrocontroller ohne FPU, daher werden Gleitkommaoperationen in Software ausgeführt. Sie sind ziemlich schnell, aber es sind immer noch etwa 100 Zyklen zum Addieren oder Multiplizieren.

Irgendwelche Vorschläge?

Hinweis: Aus Gründen der Übersichtlichkeit in der Gleichung im obigen Code können wir die Angabe weglassen 1.0 als explizites Fließkommaliteral.

Question 2

Wie Jason C in den Kommentaren betont, ist die von Ihnen gepostete Version aufgrund ihrer überlegenen Präzision in der Nähe der Randfälle höchstwahrscheinlich die beste Wahl:

float lerp(float a, float b, float f)
{
    return a * (1.0 - f) + (b * f);
}

Wenn wir die Genauigkeit für eine Weile außer Acht lassen, können wir den Ausdruck wie folgt vereinfachen:

a(1 − f) × (b − a)
= a − af + bf

= a + f(b − a)

Das heißt, wir könnten es so schreiben:

float lerp(float a, float b, float f)
{
    return a + f * (b - a);
}

In dieser Version haben wir eine Multiplikation entfernt, aber etwas Genauigkeit verloren.

Question 3

Unter der Annahme, dass Gleitkommamathematik verfügbar ist, ist der Algorithmus des OP gut und der Alternative immer überlegen a + f * (b - a) aufgrund von Präzisionsverlust, wenn a und b unterscheiden sich deutlich in der Größenordnung.

Zum Beispiel:

// OP's algorithm
float lint1 (float a, float b, float f) {
    return (a * (1.0f - f)) + (b * f);
}

// Algebraically simplified algorithm
float lint2 (float a, float b, float f) {
    return a + f * (b - a);
}

In diesem Beispiel gehen wir von 32-Bit-Gleitkommazahlen aus lint1(1.0e20, 1.0, 1.0) wird korrekt 1.0 zurückgeben, wohingegen lint2 gibt fälschlicherweise 0.0 zurück.

Der Großteil des Genauigkeitsverlusts liegt in den Additions- und Subtraktionsoperatoren, wenn sich die Operanden in der Größe erheblich unterscheiden. Im obigen Fall sind die Übeltäter die Subtraktion in b - aund der Zusatz in a + f * (b - a). Der Algorithmus des OP leidet darunter nicht, da die Komponenten vor der Addition vollständig multipliziert werden.

Für die a=1e20, b=1 Fall, hier ist ein Beispiel für unterschiedliche Ergebnisse. Testprogramm:

#include <stdio.h>
#include <math.h>

float lint1 (float a, float b, float f) {
    return (a * (1.0f - f)) + (b * f);
}

float lint2 (float a, float b, float f) {
    return a + f * (b - a);
}

int main () {
    const float a = 1.0e20;
    const float b = 1.0;
    int n;
    for (n = 0; n <= 1024; ++ n) {
        float f = (float)n / 1024.0f;
        float p1 = lint1(a, b, f);
        float p2 = lint2(a, b, f);
        if (p1 != p2) {
            printf("%i %.6f %f %f %.6e\n", n, f, p1, p2, p2 - p1);
        }
    }
    return 0;
}

Ausgabe, leicht angepasst für die Formatierung:

    f            lint1               lint2             lint2-lint1
0.828125  17187500894208393216  17187499794696765440  -1.099512e+12
0.890625  10937500768952909824  10937499669441282048  -1.099512e+12
0.914062   8593750447104196608   8593749897348382720  -5.497558e+11
0.945312   5468750384476454912   5468749834720641024  -5.497558e+11
0.957031   4296875223552098304   4296874948674191360  -2.748779e+11
0.972656   2734375192238227456   2734374917360320512  -2.748779e+11
0.978516   2148437611776049152   2148437474337095680  -1.374390e+11
0.986328   1367187596119113728   1367187458680160256  -1.374390e+11
0.989258   1074218805888024576   1074218737168547840  -6.871948e+10
0.993164    683593798059556864    683593729340080128  -6.871948e+10
1.000000                     1                     0  -1.000000e+00

Question 4

Wenn Sie sich auf einem Mikrocontroller ohne FPU befinden, wird Gleitkomma sehr teuer. Könnte für eine Gleitkommaoperation leicht zwanzigmal langsamer sein. Die schnellste Lösung besteht darin, die gesamte Berechnung mit Ganzzahlen durchzuführen.

Die Anzahl der Stellen nach dem festen Binärkomma (http://blog.credland.net/2013/09/binary-fixed-point-explanation.html?q=fixed+binary+point) ist: XY_TABLE_FRAC_BITS.

Hier ist eine Funktion, die ich verwende:

inline uint16_t unsignedInterpolate(uint16_t a, uint16_t b, uint16_t position) {
    uint32_t r1;
    uint16_t r2;

    /* 
     * Only one multiply, and one divide/shift right.  Shame about having to
     * cast to long int and back again.
     */

    r1 = (uint32_t) position * (b-a);
    r2 = (r1 >> XY_TABLE_FRAC_BITS) + a;
    return r2;    
}

Mit der Funktion inline sollte es ca. 10-20 Zyklen.

Wenn Sie einen 32-Bit-Mikrocontroller haben, können Sie größere Ganzzahlen verwenden und größere Zahlen oder mehr Genauigkeit erhalten, ohne die Leistung zu beeinträchtigen. Diese Funktion wurde auf einem 16-Bit-System verwendet.

Question 5

Wenn Sie für einen Mikrocontroller ohne Gleitkommaoperationen codieren, ist es besser, überhaupt keine Gleitkommazahlen zu verwenden und zu verwenden Festpunktarithmetik stattdessen.

Question 6

Seit C++20 können Sie verwenden std::lerp()was wahrscheinlich die bestmögliche Implementierung für Ihr Ziel ist.

Question 7

Es sei darauf hingewiesen, dass die standardmäßigen linearen Interpolationsformeln f1

Diese Funktion hat für mich auf Prozessoren funktioniert, die IEEE754-Gleitkomma unterstützen, wenn ich brauche, dass sich die Ergebnisse gut verhalten und die Endpunkte genau treffen (ich verwende sie mit doppelter Genauigkeit, aber Float sollte auch funktionieren):

double lerp(double a, double b, double t) 
{
    if (t <= 0.5)
        return a+(b-a)*t;
    else
        return b-(b-a)*(1.0-t);
}

Question 8

Wenn Sie möchten, dass das Endergebnis eine Ganzzahl ist, ist es möglicherweise schneller, auch Ganzzahlen für die Eingabe zu verwenden.

int lerp_int(int a, int b, float f)
{
    //float diff = (float)(b-a);
    //float frac = f*diff;
    //return a + (int)frac;
    return a + (int)(f * (float)(b-a));
}

Dadurch werden zwei Casts und ein Float multipliziert. Wenn ein Cast auf Ihrer Plattform schneller ist als ein Float-Addieren/Subtrahieren und wenn eine Integer-Antwort für Sie nützlich ist, könnte dies eine vernünftige Alternative sein.