【问题标题】:Floating point linear interpolation浮点线性插值
【发布时间】:2011-05-20 04:56:25
【问题描述】:

要在两个变量 ab 之间进行线性插值,给定一个分数 f,我目前正在使用此代码:

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

我认为可能有一种更有效的方法。我使用的是没有 FPU 的微控制器,因此浮点运算是在软件中完成的。它们的速度相当快,但相加或相乘仍然需要 100 个周期。

有什么建议吗?

n.b.为了清楚上面代码中的等式,我们可以省略将 1.0 指定为显式浮点文字。

【问题讨论】:

    标签: c algorithm embedded interpolation linear-interpolation


    【解决方案1】:

    忽略精度差异,该表达式等效于

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

    这是 2 次加法/减法和 1 次乘法,而不是 2 次加法/减法和 2 次乘法。

    【讨论】:

    • 当 a 和 b 的指数显着不同时,由于精度损失,这不是等效算法。 OP 的算法总是更好的选择。例如,此答案中的算法 lerp(-16.0e30, 16.0, 1.0) 将返回 0,而不是 OP 算法产生的正确结果 16。精度损失发生在加法运算符中,当a 明显大于f * (b - a) 时,以及在(b - a) 中的减法运算符中。
    • 原始算法在性能方面也没有太大损失:FP乘法比FP加法快得多,如果f保证在0和1之间,对@987654327进行某些优化@ 是可能的。
    • @Sneftel:您能详细说明一下1 - f 的优化吗?我恰好处于那种情况并且很好奇:D
    • @coredump 抱歉 2 年前没有注意到您的评论(呵呵...)。 OP 仍然会更精确,特别是如果f * (b - a) 在该算法中的幅度与a 有显着差异,那么加法就会分崩离析。这是您遇到麻烦的加法/减法。也就是说,如果f 相对于1.0f 太大,即使是OP 也会失败,因为对于非常大的f1.0f - f 可能等同于-f。因此,如果您使用 f 的巨大值,您需要认真考虑一下数学。问题是你遇到了1.0 + 1.0e800 == 1.0e800之类的东西。
    • 把浮点数想象成定点尾数和指数(它比这更复杂,但以这种方式查看它们足够发现许多 麻烦区域)。因此,如果您超出尾数的精度,您将开始丢失信息。在概念上类似于我们不能用十进制表示 1,230,000 只有两个有效数字(1.2 * 10^6 是我们能得到的最接近的),所以如果你做 1,200,000 + 30,000 但你只有两个有效数字你的处置,你失去了那 30,000。
    【解决方案2】:

    假设浮点数学可用,则 OP 的算法是一个很好的算法,并且由于在 ab 的幅度显着不同时精度损失,OP 的算法总是优于替代 a + f * (b - a)

    例如:

    // 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);
    }
    

    在该示例中,假设 32 位浮点数 lint1(1.0e20, 1.0, 1.0) 将正确返回 1.0,而 lint2 将错误地返回 0.0。

    当操作数的大小差异很大时,大部分精度损失出现在加法和减法运算符中。在上述情况下,罪魁祸首是b - a 中的减法和a + f * (b - a) 中的加法。 OP 的算法不会受此影响,因为在相加之前组件已完全相乘。


    对于 a=1e20, b=1 的情况,以下是不同结果的示例。测试程序:

    #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;
    }
    

    输出,格式稍作调整:

    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

    【讨论】:

    • 有趣的是,OP 的版本并不总是更好。我以为它被这个例子咬了:lerp(0.45, 0.45, 0.81965185546875)。它显然应该给出 0.45,但至少对于双精度,我得到 0.45000000000000007,而当 a==b 时,a + (b-a)*f 版本显然给出了 a。我希望看到一种算法具有以下属性:lerp(a, b, f) 如果f==0,则返回a,如果f==1,则b,并且对于ab] 保持在@987654339 的范围内@ 在 [0,1] 中。
    • 首先,您需要案例if a == b -&gt; return a。但是,精确的 0.45 不可能以双精度或浮点精度表示,因为它不是 2 的精确幂。在您的示例中,所有参数 a, b, f 在函数调用内部时都存储为双精度 - 返回 a 永远不会返回正好 0.45。 (当然,对于像 C 这样的显式类型语言)
    • 这看起来是更好的选择。有趣的是,标准库 lerp 似乎与 algebraically simplified version 一起使用。想法?
    • @Don Well;事实是相关的,因为它是本观察的关键;被忽视的是,它与 lerp 实现的联系是一条红鲱鱼:是的 lerp(a, a, anything) 应该返回 a,但 0.45 无法表示,因此 在该函数的域之外,所以谈论它是没有意义的。另请注意,两个版本的 lerp 都不会精确地产生 0.45。即使return 0.45 也不会返回 0.45。不过,使用此类语言的程序员通常不会在谈话中提及这一点,因为它通常是含蓄且无趣的。
    • @LorahAttkins 而 C++ 标准将 std::lerp 指定为计算 $a+t(b-a)$,这用作函数计算内容的数学定义。该标准还对std::lerp 的实现施加了更多限制:它必须是单调的,对于$t\in\{0,1\}$ 和$a = b$ 必须是精确的。这意味着lint1lint2 都不是std::lerp 的有效实现。因此,没有人会使用std::lerp,因为它太长太慢了。
    【解决方案3】:

    如果您使用的是没有 FPU 的微控制器,那么浮点将非常昂贵。浮点运算很容易慢 20 倍。最快的解决方案是使用整数进行所有数学运算。

    固定二进制点(http://blog.credland.net/2013/09/binary-fixed-point-explanation.html?q=fixed+binary+point)后的位数为:XY_TABLE_FRAC_BITS。

    这是我使用的一个函数:

    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;    
    }
    

    使用内联函数,它应该是大约。 10-20 个周期。

    如果您有一个 32 位微控制器,您将能够使用更大的整数并获得更大的数字或更高的精度,而不会影响性能。此函数用于 16 位系统。

    【讨论】:

    • 我阅读了该网站,但对于应该在什么位置仍然有点困惑。这是 0 到 0xFFFF 的值吗?还是 0 到 0xFFFE?还有什么是 XY_TABLE_FRAC_BITS? 8 个?
    • @jjxtra: XY_TABLE_FRAC_BITS 只是(可怜的)命名整数常量,其值指定假定的二进制点在使用的定点整数值中的位置(因为它不会“浮动”就像在浮点数中一样)。
    【解决方案4】:

    如果您正在为没有浮点运算的微控制器进行编码,那么最好不要使用浮点数,而是使用fixed-point arithmetic

    【讨论】:

    • 我打算迁移到定点,但是浮点已经很快了。
    【解决方案5】:

    值得注意的是,标准线性插值公式f1(t)=a+t(ba)、f2(t)=b-(ba)(1-t)、f3(t)=a (1-t)+bt 不保证在使用浮点运算时表现良好。 即,如果 a != b,则不保证 f1(1.0) == b 或 f2(0.0) == a,而对于 a == b,不保证 f3(t) 等于 a , 当 0

    这个函数在支持 IEEE754 浮点的处理器上对我有用,当我需要结果表现良好并准确地达到端点时(我以双精度使用它,但浮点也应该工作):

    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);
    }
    

    【讨论】:

    • 在 c++20 中他们添加了 std::lerp,它保证了单调的行为。
    • 这似乎是我见过的最好的解决方案。我想看看证明它是单调的。 (似乎是,因为我找不到反例,但我不明白为什么。)
    • @DonHatch 根据您的要求更改了措辞。谢谢!
    • @DonHatch 我暂时从答案中删除了“单调”,因为我没有证据。
    • 哦,但单调性是最好的部分! :-) 显然 f1 和 f2 两部分是单调的,仍有待证明它在切换点 t=0.5 处是单调的。我认为是(仅从我寻找反例失败的事实来看),只是还没有证明。也许这对于其他一些更有理论意识的网站(例如 cs.stackechange.com)来说是一个很好的问题。请注意,那里有一个相关问题:cs.stackexchange.com/questions/59625/…
    【解决方案6】:

    从 C++20 开始,您可以使用 std::lerp(),这可能是您目标的最佳实现。

    【讨论】:

    • std::lerp 在我看来完全不应该被使用。您实际上很少需要插值外推,以及大量的分支行为,数值不稳定的内部实现之上。我对std::lerp 的实现方式有很多不同意见,很难推荐。
    • @jeremyong 你能举一个std::lerp 表现不佳的例子吗?它的合约在几个重要方面看起来确实不错:它是单调的,lerp(a,b,0)==a, lerp(a,b,1)==b(这两个事实意味着它保持在范围内 [ a,b] 对于 [0,1] 中的 t),lerp(a,a,t)==a。所以通常的抱怨似乎都被覆盖了。
    【解决方案7】:

    如果您希望最终结果为整数,则将整数也用于输入可能会更快。

    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));
    }
    

    这会进行两次强制转换和一次浮点数相乘。如果在您的平台上强制转换比浮点加/减更快,并且整数答案对您有用,那么这可能是一个合理的选择。

    【讨论】:

    • 对于f * (b - a),类型提升将允许将(b - a) 提升为float,因为ffloat 类型。所以,在(float)(b - a) 中对(float) 的显式转换充其量只是说明性的,但实际上没有必要,不是吗?
    • @Scheff - 是的,你是对的,浮点转换纯粹是为了引起人们注意编译器无论如何都会插入的东西。
    猜你喜欢
    • 1970-01-01
    • 2017-02-22
    • 1970-01-01
    • 2013-07-27
    • 2015-09-25
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多