浮点数组的更快 abs-max答案

【问题标题】：Faster abs-max of float array浮点数组的更快 abs-max
【发布时间】：2010-10-23 11:43:08
【问题描述】：

我需要实时绘制音频峰值表。每秒最少 44100 个样本乘以最少 40 个流。每个缓冲区介于 64 到 1024 个样本之间。我需要从每个缓冲区中获取最大绝对值。（然后将它们通过一种低通滤波器并以大约 20 毫秒的间隔进行绘制。）

for(int i = 0; i < numSamples; i++)
{ 
      absMaxOfBuffer = MAX( fabs( buffer[i] ), absMaxOfBuffer);
}

我现在就是这样做的。我想做得更快。缓冲区的浮点数在 -1 到 1 范围内，因此是晶圆厂。

问题，有没有一些棘手的 comp-sci quicksort-esque 方法可以更快地完成这项工作？

如果做不到这一点，浮点数的无分支 ABS 和 MAX 函数是否存在？

编辑：主要平台是 Linux/gcc，但计划有一个 windows 端口（可能与 mingw）。

编辑，第二个：
我接受了一个人，因为关于实际算法结构的一点是问题的核心。
我将尝试将循环展开到四个，将符号位归零，然后使用 SSE（maxps 指令）获得最大值，看看是否不会剥皮。感谢您的建议，我投票赞成你们中的一些人，作为亚军。 :)

【问题讨论】：

什么编译器？什么平台？如果你使用 gcc，你可能会对 __builtin_expect 之类的东西感兴趣。

标签： c++ performance optimization signal-processing

【解决方案1】：

晶圆厂和比较对于 IEEE 浮点数都非常快（例如，原则上单整数运算速度很快）。

如果编译器没有内联这两个操作，那么要么戳它直到它内联，要么找到适合您架构的实现并自己内联。

您也许可以从 positive IEEE 浮点数与具有相同位模式的整数的顺序相同的事实中得到一些信息。也就是说，

f > g   iff   *(int*)&f > *(int*)&g

因此，一旦您完成了 fabs'ed，我认为 int 的无分支最大值也适用于 float（当然假设它们的大小相同）。有一个解释为什么在这里工作：http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm。但是你的编译器已经知道这一切，你的 CPU 也知道，所以它可能没有任何区别。

没有比复杂性更快的方法了。您的算法已经是 O(n)，您无法击败它并仍然查看每个样本。

我猜你的处理器的 SIMD（即 Intel 上的 SSE2）中可能有一些东西会有所帮助，因为每个时钟周期处理的数据比你的代码多。但我不知道是什么。如果有，那么它很可能会快几倍。

您可能可以在多核 CPU 上进行并行化，尤其是因为您要处理 40 个独立的流。这将是最好的几个因素更快。 “只需”启动适当数量的额外线程，将它们之间的工作分开，并使用最轻量的原语来指示它们何时全部完成（可能是线程屏障）。我不太清楚您是在绘制所有 40 个流的最大值，还是分别绘制每个流的最大值，所以除了确保将结果传递到下一阶段之外，您实际上可能不需要同步工作线程没有数据损坏。

可能值得看一下反汇编，看看编译器展开了多少循环。尝试再展开它，看看是否有什么不同。

另一件需要考虑的事情是您获得了多少缓存未命中，以及是否可以通过向缓存提供一些线索来减少该数量，以便它可以提前加载正确的页面。但是我没有这方面的经验，我不会抱太大希望。 __builtin_prefetch 是 gcc 上的魔法咒语，我猜第一个实验类似于“在进入该块的循环之前预取下一个块的开头”。

您目前处于所需速度的百分之几？还是“尽可能快”的情况？

【讨论】：

“已经太多了”。只要这个函数已经被多个线程使用。否则，如果一个线程正在完成所有工作，那么无论存在的原始线程数量如何，您都没有充分利用 CPU。你可以有 1000 个线程，他们都在等待这个线程然后你没有足够的线程 ;-)
正确，我不知道更好的算法，假设您必须查看每个样本。您也许可以尝试只取每个缓冲区前半部分的最大值，看看您出错的频率和程度，然后决定您是否关心。当然，如果波形被截断，您知道它不会变得更大，并且可以忽略该缓冲区的其余部分。但是要检测到您需要在循环中进行额外的比较，这会在未剪辑的情况下减慢速度。而且我假设浮动音频数据不一定很明显构成“剪辑”的幅度。

【解决方案2】：

http://www.scribd.com/doc/2348628/The-Aggregate-Magic-Algorithms 上记录了一个无分支晶圆厂

还请注意，最新版本的 GCC 将使用 MMX 指令为您内联无分支 fabs。还有fmin 和fmax，但GCC 不会内联它们（你会得到call fmin）。

【讨论】：

【解决方案3】：

要尝试的事情：

fabs() 可能不是内联函数。
特殊的组装说明可能会有所帮助。在 Intel 上，SSE 有一条指令可以一次计算最多四个浮点数。
否则，IEEE 754 规范规定如果 a 和 b 非负浮动，则 a
您真的需要每个样本的最大值吗？也许最大值可能不止一次出现，从而可能导致不检查每个输入。
您能以流方式计算最大值吗？

【讨论】：

我试图不去假设任何事情。
或者只是构建而不首先组装（使用gcc我认为这是-S？）。

【解决方案4】：

您可能想查看Eigen。

它是一个 C++ 模板库，使用 SSE（2 及更高版本）和 AltiVec 指令集，优雅地回退到非矢量化代码。

快。（参见基准）。
表达式模板允许在适当的时候智能地删除临时变量并启用惰性求值 - Eigen 会自动处理此问题并在大多数情况下也处理别名。
对 SSE（2 和更高版本）和 AltiVec 指令集执行显式矢量化，并优雅地回退到非矢量化代码。表达式模板允许对整个表达式全局执行这些优化。
使用固定大小的对象，可以避免动态内存分配，并且在有意义时展开循环。
对于大型矩阵，需要特别注意缓存友好性。

【讨论】：

+1 缓存行友好性。还有另一个模板库：nt2.sourceforge.net

【解决方案5】：

用另一个问题回答一个问题并不完全是在回答，但是嘿...我也不是 C++ 开发人员。

既然你是用 C++ 开发这个并且你正在做 dsp，你不能连接到 matlab 或 octave（它是开源的）来为你做数学并得到结果吗？

在这些软件中已经实现了一些功能（如 conv、fft、ifft、fir 和绘图实用程序功能，如 freqz、stem、graph、plot、mesh 等）。我在 Photoshop 中查看了一个名为 MATLAB 的大文件夹……其中包含一些从应用程序获取调用的 .m 文件，将其发送到动态 matlab，然后将结果返回给应用程序。

希望对你有帮助。

【讨论】：

【解决方案6】：

我看到的简单优化：

将循环转换为指针算术版本（假设您的编译器看不到这个）
IEEE 754 标准将其表示定义为符号幅度。因此，将最高有效位设置为 0 与调用 fabs() 相同。也许这个函数已经使用了这个技巧。

【讨论】：

【解决方案7】：

出于您的目的，您可以将其平方而不是取绝对值；正如数学上的 |a|

【讨论】：