数字信号处理中的 SIMD

Posted 2023-02-23 芥末的无奈

系列文章目录

• 数字信号处理中的 SIMD
• Neon intrinsics 简明教程

---

文章目录

• 系列文章目录
• 0.
• 1. 前言
• 2. SIMD 是什么
• 3. SIMD 伪代码示例
• 4. SIMD 是如何被实现的？
• 5. SIMD 中哪些指令是可用的？
• • Load/Store
  • 寄存器上的操作
• 6. 如何获取 SIMD 指令集
• 7. MMX，SSE，AVX，NEON...
• 8. 选择哪种 SIMD ？
• 9. 为什么 SIMD 在 DSP 中很有用？
• 10. SIMD 的缺点
• 11. SIMD 简易示例
• 12. 总结
• 13. 参考

---

0.

本文翻译自 What is SIMD in Digital Signal Processing?，希望通过对这篇文章的翻译了解 SIMD，并为后续学习 SIMD 做铺垫。强烈建议有余力的同学可以阅读原文。

---

1. 前言

图像或声音的数字信号处理需要对大量的数据进行复杂的操作。例如，为了缩放（改变音量）一秒钟的音频数据，我们可能要进行44100次乘法运算。

如果我们实时地进行操作，且整个渲染过程在 10ms 内完成操作，事情就变得更加困难了。

值得庆幸的是，有一些编程工具可以让我们更高效地处理这些情况，其中之一就是本文的主题 — SIMD

2. SIMD 是什么

Single instruction, multiple data(SIMD)，即「单指令多数据」，它是一种特殊的处理器指令，一次对 1 个以上的变量进行操作。

用数学术语来说，我们可以说 SIMD 对变量的向量（数组）进行操作，就像正常代码对单个变量进行操作一样。

3. SIMD 伪代码示例

想象一下，我们有两个数组 $v_1$ 和 $v_2$，每个数组有 8 个 32 位的浮点数变量，假设 $v_{12}=v_1+v_2$。如何计算 $v_{12}$？在 C 代码中，你可以这么写：

for(auto i = 0; i < 8; ++i)
	v12[i] = v1[i] + v2[i]

上述代码中，有 8 次加法操作。如果我们可以简单地这么写：

v12 = vector8x32_add(v1, v2);

这样就只会产生一个加法操作？从理论上讲，这样的程序会快 8 倍。实际上，这与事实相差无几! 虽然执行速度的提高不一定与多线程的使用成正比，但SIMD的速度往往与数据向量中可容纳的变量数量成正比。

在AVX-512指令中，人们可以用一条指令对16个浮点数进行操作。你能想象你的代码的处理时间减少16倍吗？而如果我们使用16位整数，我们乐观地得到…没错，是32倍的提速。这就是 SIMID 的力量。

注意：使用SIMD的实现的性能必须在每个处理器上单独测量。此外，在第一次尝试时，线性提速是非常不可能的。但这是另一篇文章的主题。

4. SIMD 是如何被实现的？

SIMD 是又微处理器产商提供的。

一般来说，有几个处理器架构系列，其中主要是 x86 和 ARM。

每个架构系列都有一套的核心寄存器和相关指令，可以在该系列的任何处理器上运行。而新一代处理器都在上一代处理器的指令集上进行扩展。此外，每个处理器可以有一些列指令，这些指令是特定处理器型号特有的。

特定世代的指令集和特定处理器的指令通常利用额外的架构元素，如额外的寄存器或专用算术单元。

你可能会问自己：当每个处理器接受不同的指令时，怎么可能编写一个版本的源代码？

这就是为什么像 C 这样的高级编程语言的存在是一个福音。 🙂

编译器的作用是把你写的源代码，翻译成处理器专用的汇编语言。编译器对于哪些处理器包含哪些指令非常明智。他们也很善于确定调用这些指令中的哪些指令，并以何种顺序来使软件尽可能地高效。事实上，它们比至少95%的程序员（包括我）更擅长于此。

不幸的是，他们不知道我们想通过代码达到什么目的。例如，他们不知道这种大量的乘法和加法实际上是有限脉冲响应（FIR）过滤。因此，他们往往无法利用底层硬件的全部潜力。

这就是你需要你参与的地方：一个程序员，他知道如何用处理器特定的指令来编码 DSP 算法。

5. SIMD 中哪些指令是可用的？

让我们再次声明。SIMD指令是简单的特殊处理器指令。这意味着它们可以在汇编代码中找到，也就是直接在寄存器和内存上操作的代码。

Load/Store

SIMD通常使用专用寄存器。这意味着，有三种指令必须是肯定可以使用的：

• 将数据从内存移动到专用寄存器（Load）
• 将数据从特殊寄存器中转移到内存（Store）
• 设置特殊寄存器的值（Set）

例如，x86 处理器的 AVX 指令集有 vmovups 指令，它将 8 个 32位浮点数（ps 表示"单精度"）的向量（v）加载（mov）到 AVX寄存器中。指令中的 u 代表 “未对齐”，这意味着我们加载的内存位置不必在32字节边界上对齐。

寄存器上的操作

一旦数据进入专用寄存器，我们就可以对它们进行一些操作：

• 在一个寄存器上进行算术运算（例如，floor, ceiling、平方根）。
• 在 2 个寄存器上进行算术运算（例如，加、减、乘、除）。
• 在 2 个寄存器上进行逻辑运算（例如，AND，OR）。
• 转换寄存器中的数据类型，例如，从整数到浮点数。
• 操作寄存器中的变量位置（例如，置换）。
• 对寄存器中的变量进行不同的操作，取决于其索引是偶数还是奇数。
• 对寄存器中的数值进行比较。
• 执行DSP操作，如乘法和加法。

正如你所看到的，可用操作的数量是巨大的。

为了在任何给定的硬件上有效地利用 SIMD 的力量，我们可以使用处理器制造商提供的文档。

例如，AVX指令集的所有可用指令的列表可以在英特尔的网站上找到。

6. 如何获取 SIMD 指令集

我们知道 SIMD 指令是什么，也知道了它能做什么，作为软件开发者，我们如何才能访问它们呢？我见过 4 中不同的方法可以将 SIMD 指令嵌入到你的代码中。

1. 自动向量化。某些情况下，编译器本身可以聪明到将没有向量化的代码进行自动向量化。
2. 使用汇编指令。可以使用汇编来写整个软件，或者在 C/C++ 中使用 asm 块
3. 内置函数。处理器制造商通常提供 C 函数来执行引擎盖下的专用处理器指令（或其组合）。程序员必须简单地包括相关的头文件，并用专用的编译器选项来编译他们的程序。
   例如，英特尔已经发布了其架构上所有可用的 SIMD 指令的列表，称为英特尔指令指南。ARM为基于ARM的体系结构发布了一个类似的列表。
   特殊编译标志的例子是 MSVC 上的 /arch:AVX 以及 gcc 和clang上的 -mavx（以便能够使用AVX指令）。问题是，运行该软件的处理器必须有这些指令。这是开发人员的一个额外责任。
4. 专门的库。有一些软件库在编写的代码和它所运行的硬件之间提供了一个抽象层。如果我们想到矢量加法，我们可以猜测每个支持 SIMD 的处理器可能都有一些矢量加法指令。它的调用方式可能不同，但功能是一样的。这种库的一个例子是开放计算语言（OpenCL）。一个更贴近音频开发者的例子是 JUCE 框架。JUCE 提供了诸如 SIMDRegister 这样的抽象，它是对平台特定的扩展寄存器类型的封装。这种方法可能是最舒服的方法，但它要求你使用（通常是付费的）第三方软件。

7. MMX，SSE，AVX，NEON…

你可能遇到过许多代表 SIMD 指令集系列的缩写。下面我列出了对数字信号处理最重要的指令以及简短的描述。

简称	全名	架构	可用寄存器	备注
Intel® MMX™		x86	8 个 64-bit 的寄存器	只支持整数操作
Intel® SSE	Streaming SIMD Extensions	x86	8 个 128-bit 的寄存器	在 MMX 上的扩展，支持浮点操作
Intel® AVX	Advanced Vector Extensions	x86	8 个 256-bit 的寄存器（AVX-512 则是 512-bit）	SSE 的扩展
Arm Neon		ARM	32 个 128-bit 的寄存器	android 6.0 以上设备全部支持；苹果设备也支持：iPhone、iPad、一些型号的 Mac；可用在 x86 架构上使用 ARM_NEON_2_x86_SSE

8. 选择哪种 SIMD ？

由于有很多指令集可供选择，你可能会发问：如何选择合适的指令集？答案总是：看情况。

首先要考虑你的软件在在什么平台上运行。例如如果是你为 Android 设备编写软件，你可能会使用 Neon。但并不总是如此，有一些 Android 驱动的设备也运行在 X86 架构上。

如果你只针对x86架构，并想尽可能地利用SIMD，你应该为所有的x86 SIMD指令集提供实现，并给软件提供后备可能性；例如，如果AVX-512不存在，也许可以使用SSE。请记住，有SSE、SSE2、SSE3、SSE4…还有一些小版本。

如果你不知道你的软件将在哪个平台上运行，你就会面临实现每个可用指令集支持的可能性。

这就是为什么编译器是如此的幸运：它们可以为我们做这些繁重的工作。而且在这个过程中犯的错误要少得多。 🙂

9. 为什么 SIMD 在 DSP 中很有用？

SIMD在数字信号处理应用中特别有优势。为什么呢？

1. DSP 算法通常是以向量来定义的。此外，科学家们习惯于在Matlab或Python中对向量进行操作。借助SIMD的力量，我们也可以在C/C++代码中有效地使用向量。
2. DSP算法经常在不同的数据上执行相同的任务。SIMD指令是为了 一次对多个变量进行相同的操作。
3. 信号处理最常在块中进行。音频样本、图像数据或电影数据块的长度通常等于 SIMD 寄存器大小的倍数。这使得它们很容易被向量化。
4. SIMD指令集通常包含 DSP 特定的功能。例如，Neon指令包含一个乘加运算。这意味着我们可以用一条指令执行点乘。

10. SIMD 的缺点

不幸的是，SIMD 并不全是美好的。下面是SIMD在DSP方面的一些缺点（但不限于此）。

1. 程序员的噩梦：支持所有的指令集。如果你建立的不仅仅是跨平台的应用程序，而且是应该在不同的处理器架构上类似地工作的应用程序，你可能会遇到确定底层架构、其特征和处理每一种可能情况的问题。想想看：你不仅要使用特定的AVX、SSE或Neon指令编写代码。你还需要到处实现编译时甚至运行时的检查，以确定使用哪些例程。这在算法代码的基础上增加了很多额外的代码，包括宏的技巧。
   不要认为如果你只写Android的代码，你的代码就会一直在ARM架构上运行。现在，有很多设备在x86系列的处理器上运行Android。这就是为什么，最好使用一个能告知你当前处理器能力的库，如谷歌的cpu_features库。
2. 运行时的可用性检查。如上所述，你有时需要在运行时确定使用哪个指令集。这增加了额外的代码和执行时间的开销。
3. 不对齐的数据。扩展指令集在对齐的数据上效果最好。不幸的是，你的数据本身通常不会被对齐。在一个特定的边界上对齐向量的必要性在你的算法上又增加了一层代码和复杂性。
4. 边缘情况，只有一个采样点的输入。如果你想处理的信号数据不是在大小等于SIMD寄存器大小的倍数的块中出现的呢？那么你就需要 "手动 "完成算法的标量版本。这意味着更多的复杂性。
5. SIMD 网上资源很少。用 SIMD 处理信号在网络上或 YouTube 上不是一个很好的解释的话题。人们需要求助于专业书籍和研究论文来了解基础知识。SIMD 是一个必须在其应用环境中讨论的话题。把图像处理或3D图形的教程转移到音频处理上并不容易。这些原因使得入门门槛相当高。
6. 可读性低。SIMD代码中充满了像 vrecpeq_f32 或 _mm256_testnzc_ps 这样的函数，当你看代码或与同事交谈时，不容易阅读和理解，也不容易发音。另一方面，一旦你对内在函数有了一些了解，这些名字就会成为非常好的记忆法。

11. SIMD 简易示例

为了完成这篇文章，我们将用C++语言编写一个小的例子，它使用了内置函数。

具体来说，我们将使用英特尔处理器上的AVX指令集。

这个小程序的目标是计算两个浮点数向量的内积。AVX寄存器有256位，所以每个寄存器应该正好可以容纳8个32位浮点数。

#include <vector>
#include <array>
#include <cassert>
#include <random>
#include <algorithm>
#include <chrono>
#include <iostream>

#include <immintrin.h>

using Vector = std::vector<float>;

Vector scalarAdd(const Vector& a, const Vector& b) 
    assert(a.size() == b.size());

    Vector result(a.size());

    for (auto i = 0u; i < a.size(); ++i) 
        result[i] = a[i] + b[i];
    

    return result;


Vector simdAdd(const Vector& a, const Vector& b) 
    assert(a.size() == b.size());

    Vector result(a.size());

    constexpr auto FLOATS_IN_AVX_REGISTER = 8u;

    const auto vectorizableSamples = (a.size() / FLOATS_IN_AVX_REGISTER) 
                                     * FLOATS_IN_AVX_REGISTER;

    auto i = 0u;
    for (; i < vectorizableSamples; i += FLOATS_IN_AVX_REGISTER) 
        // load unaligned data to SIMD registers
        auto aRegister = _mm256_loadu_ps(a.data() + i);
        auto bRegister = _mm256_loadu_ps(b.data() + i);

        // perform the addition
        auto intermediateSum = _mm256_add_ps(aRegister, bRegister);

        // store data back in the data vector
        _mm256_storeu_ps(result.data() + i, intermediateSum);
    
    // process the remaining (unvectorized) samples
    for (; i < a.size(); ++i) 
        result[i] = a[i] + b[i];
    

    return result;


Vector randomVector(Vector::size_type size) 
    Vector v(size);

    auto randomEngine = std::default_random_engine();
    std::uniform_real_distribution<float> uniformDistribution(-1.f, 1.f);
    auto generator = [&]()  return uniformDistribution(randomEngine); ;
    std::generate(v.begin(), v.end(), generator);

    return v;


int main() 
    constexpr auto TEST_VECTOR_SIZE = 1000001;

    const auto a = randomVector(TEST_VECTOR_SIZE);
    const auto b = randomVector(TEST_VECTOR_SIZE);

    assert(scalarAdd(a, b) == simdAdd(a, b));

    constexpr auto TEST_RUN_COUNT = 1000;

    using namespace std::chrono;
    milliseconds totalScalarTime;
    for (auto i = 0u; i < TEST_RUN_COUNT; ++i) 
        auto start = high_resolution_clock::now();
        auto result = scalarAdd(a, b);
        auto end = high_resolution_clock::now();
        totalScalarTime += duration_cast<milliseconds>(end - start);
    

    std::cout << "Average scalarAdd() execution time: " 
              << totalScalarTime.count() / static_cast<float>(TEST_RUN_COUNT) 
              << " ms." << std::endl;

    milliseconds totalSimdTime;
    for (auto i = 0u; i < TEST_RUN_COUNT; ++i) 
        auto start = high_resolution_clock::now();
        auto result = simdAdd(a, b);
        auto end = high_resolution_clock::now();
        totalSimdTime += duration_cast<milliseconds>(end - start);
    

    std::cout << "Average simdAdd() execution time: " 
              << totalSimdTime.count() / static_cast<float>(TEST_RUN_COUNT) 
              << " ms." << std::endl;

输出：

.\\InnerProductSIMD.exe 
Average scalarAdd() execution time: 11.695 ms.
Average simdAdd() execution time: 4.113 ms.

12. 总结

在这篇文章中，我们讨论了 SIMD 指令在数字信号处理中的用处。

SIMD指令让我们使用专用的处理器寄存器同时对一个以上的变量进行操作。

不同的处理器架构和型号有不同的SIMD指令可用。

主要的收获应该是。SIMD指令可以使你的DSP代码大大加快，但代价是

1. 代码的复杂性。
2. 可移植性。
3. 处理器方面的专家知识。

13. 参考

• Intel® SIMD technology comparison
• Intel® Intrinsics Reference
• Arm Neon Intrinsics Reference
• Advanced Vector Extensions (AVX) on Wikipedia