避免浮点运算

Posted 2023-02-16

【中文标题】避免浮点运算

【英文标题】：Avoiding floating point arithmetic

【发布时间】：2009-09-18 15:42:09

【问题描述】：

我为 iPhone 编写了一个小型软件合成器。 为了进一步调整性能，我使用 Shark 测量了我的应用程序，发现我在 float/SInt16 转换中浪费了很多时间。 因此，我重写了一些部分以通过预先计算返回“即用型”SInt16 样本的查找表来绕过转换。到目前为止，这工作正常。 目前，我正在尝试重写一些过滤器和我的 ADSR 包络实现以仅使用整数算术，但我可以使用一些提示如何在没有浮点数的情况下执行乘法/除法。 我的目标是iPhone canonical format：

LPCM 16 位整数样本

在不使用浮点数的情况下将振幅应用于最终样本的好方法是什么？

编辑： 到目前为止，我唯一想到的是，我可以通过右移当前样本来除以 2 的幂。

inBuffer[frame] = wavetable[i % cycleLengthInSamples] >> 4;

但我想不出任何优雅的方式来创建一个平滑的 ADSR 信封。

编辑2： 感谢您的所有精彩回答！ 我目前的做法：

带上我所有的 ADSR 包络值 进入正 SInt16 范围 乘以波表中的当前值（将中间值存储为 SInt32） 将结果右移 16

这似乎有效:)

【参考方案1】：

固定点很好，因为在这种情况下您使用的是 16 位。最简单的方法是根据您需要的精度乘以 10 的幂。如果您可以使用 32 位整数作为中间值，您应该能够获得不错的精度。最后，您可以转换回 16 位 int，根据需要进行舍入或截断。

编辑： 您想向左移动，以使值更大。将移位结果存储在更精确的类型中（32 位或 64 位，具体取决于您的需要）。 如果您使用带符号的类型，简单的转换将不起作用

请注意您是否将两个定点数相乘或相除。乘以 (a*n) * (bn) 最终得到 abn^2 而不是 abn。除法是 (an) / (bn) 是 (a/b) 而不是 ((an)/b)。这就是为什么我建议使用 10 的幂，如果您不熟悉定点，可以很容易地发现您的错误。

当您完成计算后，您移回右侧以返回 16 位 int。 如果你想花哨，你也可以在移位前进行四舍五入。

如果您真的对实现高效的定点感兴趣，我建议您阅读一下。 http://www.digitalsignallabs.com/fp.pdf

【讨论】：

我就是这么想的。但是你不会乘以 2 的幂来得到你的固定点吗？这只涉及移位。

当你使用有符号整数时，移位会变得很时髦。 2 的幂更有效，但更难调试。如果是第一次使用定点，我建议使用 10 的幂。

感谢您的回答！我还没有找到一种优雅的方式来实现无浮动 ADSR 信封。我只是尝试右移样本，这会导致除以 2 的任意幂，因此会降低幅度 - 但我不知道如何用它创建平滑的包络。

【参考方案2】：

this SO question 的答案在实施方面非常全面。这里的解释比我在那里看到的要多一点：

一种方法是强制所有数字进入一个范围，例如 [-1.0,1.0)。然后将这些数字映射到范围 [-2^15,(2^15)-1] 中。例如，

Half = round(0.5*32768); //16384
Third = round((1.0/3.0)*32768); //10923

这两个数相乘得到

Temp = Half*Third; //178962432
Result = Temp/32768; //5461 = round(1.0/6.0)*32768

在最后一行中除以 32768 是 Patros 所做的关于需要额外缩放步骤的乘法运算。如果您明确编写 2^N 缩放，这将更有意义：

x1 = x1Float*(2^15);
x2 = x2Float*(2^15);
Temp = x1Float*x2Float*(2^15)*(2^15);
Result = Temp/(2^15); //get back to 2^N scaling

这就是算术。对于实现，请注意两个 16 位整数的乘法需要 32 位结果，因此 Temp 应该是 32 位。此外，32768 不能用 16 位变量表示，因此请注意编译器将生成 32 位立即数。正如您已经指出的那样，您可以转换为乘/除以 2 的幂，这样您就可以编写

N = 15;
SInt16 x1 = round(x1Float * (1 << N));
SInt16 x2 = round(x2Float * (1 << N));
SInt32 Temp = x1*x2;
Result = (SInt16)(Temp >> N);
FloatResult = ((double)Result)/(1 << N);

但是假设 [-1,1) 不是正确的范围？如果您希望将数字限制为 [-4.0,4.0)，则可以使用 N = 13。然后您有 1 个符号位，二进制点之前有 2 位，之后有 13 位。这些分别称为 1.15 和 3.13 定点小数类型。你用分数的精度换取净空。

只要注意饱和度，添加和减去小数类型就可以正常工作。正如帕特罗斯所说，对于除法，缩放实际上抵消了。所以你必须这样做

Quotient = (x1/x2) << N;

或者，为了保持精度

Quotient = (SInt16)(((SInt32)x1 << N)/x2); //x1 << N needs wide storage

乘以和除以整数正常工作。例如，要除以 6，您可以简单地写

Quotient = x1/6; //equivalent to x1Float*(2^15)/6, stays scaled

在除以2的幂的情况下，

Quotient = x1 >> 3; //divides by 8, can't do x1 << -3 as Patros pointed out

但是，添加和减去整数并不天真。您必须首先查看整数是否适合您的 x.y 类型，制作等效的小数类型，然后继续。

我希望这对这个想法有所帮助，请查看其他问题中的代码以获得干净的实现。

【参考方案3】：

看看这个描述快速乘法算法的页面。

http://www.newton.dep.anl.gov/askasci/math99/math99199.htm

【参考方案4】：

一般来说，假设您将使用带符号的 16.16 定点表示。这样一个 32 位整数将有一个有符号的 16 位整数部分和一个 16 位小数部分。那我不知道iPhone开发用的是什么语言（也许是Objective-C？），但是这个例子是C语言的：

#include <stdint.h>

typedef fixed16q16_t int32_t ;
#define FIXED16Q16_SCALE 1 << 16 ;

fixed16q16_t mult16q16( fixed16q16_t a, fixed16q16_t b )

    return (a * b) / FIXED16Q16_SCALE ;


fixed16q16_t div16q16( fixed16q16_t a, fixed16q16_t b )

    return (a * FIXED16Q16_SCALE) / b ;

请注意，以上是一个简单的实现，并且不提供算术溢出保护。例如，在 div16q16() 中，我在除法之前进行多次以保持精度，但根据操作数，操作可能会溢出。您可以使用 64 位中间体来克服这个问题。此外，除法总是向下舍入，因为它使用整数除法。这提供了最佳性能，但可能会影响迭代计算的精度。修复很简单，但会增加开销。

请注意，当乘以或除以恒定的 2 次幂时，大多数编译器会发现微不足道的优化并使用移位。然而，C 并没有定义负符号整数右移的行为，所以为了安全和可移植性，我把它留给编译器来解决。 YMV 在您使用的任何语言上。

在 OO 语言中，fixed16q16_t 自然会成为具有运算符重载的类的候选对象，因此您可以像使用普通算术类型一样使用它。

您可能会发现在类型之间进行转换很有用：

double fixed16q16_to_double( fixed16q16_t fix )

    return (double)fix / FIXED16Q16_SCALE ;


int fixed16q16_to_int( fixed16q16_t fix )

    // Note this rounds to nearest rather than truncates
    return ((fix + FIXED16Q16_SCALE/2)) / FIXED16Q16_SCALE ;


fixed16q16_t int_to_fixed16q16( int i )

    return i * FIXED16Q16_SCALE ;


fixed16q16_t double_to_fixed16q16( double d )

    return (int)(d * FIXED16Q16_SCALE) ;

这是基础，可以变得更复杂并添加三角函数和其他数学函数。

固定的加法和减法适用于内置的 + 和 - 运算符及其变体。