系列文章目录

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

文章目录

• 系列文章目录
• 写在前面
• 前言
• FIR 滤波器
• 优化方案
• 何时应该使用时域的滤波器？
• 如何加快FIR滤波器的实现速度？
• 初步假设
• 有限长度信号
• 反转的滤波器系数
• 实际的卷积公式
• 卷积过程可视化
• 线性卷积最简单的实现
• 循环向量化
• Inner loop vectorization (VIL)
• VIL NEON 实现
• Outer Loop Vectorization (VOL)
• VOL NEON Implementation
• Outer and Inner Loop vectorization (VOIL)
• 为什么 VOIL 更优？
• 数据对齐
• VOIL NEON Implementation
• 总结
• 参考

写在前面

本文多数内容翻译自 Efficient FIR Filter Implementation with SIMD。原文在 SIMD 代码实现中使用到了 AVX，本文将使用 NEON 实现，关于 NEON 如何使用，请参考 Neon intrinsics 简明教程。

前言

如何让你的 FIR 滤波器在时域中更快的运行？

FIR 滤波器是数字信号处理中的基石。它在将混响应用于音频信号时尤其重要，例如在虚拟现实音频或数字音频工作站的VST插件中。它还被广泛用于移动电话（甚至是前智能手机！）和嵌入式设备的声音应用。

如何让 FIR 运行地更快呢？

FIR 滤波器

FIR 滤波器由它的有限长脉冲响应信号 h[n]h[n]h[n] 定义，FIR 滤波器的输出 y[n]y[n]y[n] 是输入信号 x[n]x[n]x[n] 与脉冲响应的卷积。我们把它写成：
y[n]=x[n]∗h[n]=∑k=−∞∞x[k]h[n−k]y[n]=x[n] * h[n]=\sum_{k=-\infty}^{\infty} x[k] h[n-k] y[n]=x[n]∗h[n]=k=−∞∑∞​x[k]h[n−k]

关于卷积的知识你可以参考我之前发的文章 【音频处理】Fast Convolution 快速卷积算法简介 或者 Convolution: The secret behind filtering

优化方案

如果你想让任何软件尽可能快地执行，有两种方法可以实现。

1. 最佳的算法
2. 高效的执行

同样的原则也适用于 DSP 代码。要想在代码中拥有一个快速的有限脉冲响应（FIR）滤波器，你可以这样做：

1. 使用一种时间复杂度更低的算法，例如通过傅里叶变换到频域后进行快速卷积，或者
2. 或者，利用硬件和软件资源的优势来有效地实现时域卷积。这通常意味着使用 SIMD 指令对你的代码进行矢量化。

何时应该使用时域的滤波器？

你或许会这么问自己：

既然我们有一个快速卷积算法，为什么我们还需要一个时域的实现呢？

快速卷积的算法复杂度是 O(Nlog⁡N)O(N\log N)O(NlogN)，其中 N 可能输入信号或滤波器的长度。而线性卷积的算法复杂度是 O(N2)O(N^2)O(N2)

首先，这意味着对于足够小的 NNN，线性卷积算法将比快速卷积算法更快。

其次，快速卷积在复数域上操作，而线性卷积总是使用实数。这实际上意味着，快速卷积需要处理的数据是线性卷积的两倍。

如果时域卷积可能比快速卷积更快这一事实让你感到困惑，请想想排序算法。一般来说，快排被认为是最快的排序算法。但是如果你看一下排序的实现，比如std::sort，你会发现它最初只使用quicksort。在将排序的容器划分为足够小的部分后，会使用插入排序对其中的元素进行排序。

在解释了这一点之后，我们可以研究如何在时域中有效地实现FIR滤波。

如何加快FIR滤波器的实现速度？

简而言之：通过 SIMD。

加快滤波器执行速度的最好方法是使用 SIMD 指令一次处理多个样本。为了实现这一点，我们需要重写线性卷积算法，使我们的代码在向量上运行。

这个过程被称为循环向量化。

循环向量化通常由编译器完成，但这种自动向量化的程度通常对实时 DSP 来说是不够的。相反，我们需要准确地指示编译器要做什么。

SIMD 指令在对对齐的数据进行操作时能达到最佳性能。因此，数据对齐是我们应该考虑的另一个因素。

总之，一个高效的FIR滤波器的实现要同时使用两种策略。

1. 循环矢量化
2. 数据对齐

现在我们将详细讨论这两种策略。

初步假设

线性卷积公式如下：
x[n]∗h[n]=∑k=−∞∞x[k]h[n−k]=y[n],n∈Z.(1)x[n] * h[n]=\sum_{k=-\infty}^{\infty} x[k] h[n-k]=y[n], \quad n \in \mathbb{Z} . \tag{1} x[n]∗h[n]=k=−∞∑∞​x[k]h[n−k]=y[n],n∈Z.(1)

正如你可能猜到的那样，一个无限长信号在实际操作中是不切实际的。此外，在代码中反转 h[n]h[n]h[n] 信号中的时间是相当有问题的。因此，我们将做一些假设，然而，这不会改变我们讨论的一般性质。

有限长度信号

我们将假设我们的信号是有限的。当然，h[n]h[n]h[n]是这样的，但x[n]x[n]x[n]不一定是如此。

我们用 NxN_xNx​ 表示 x[n]x[n]x[n] 的长度，用 NhN_hNh​ 表示 h[n]h[n]h[n] 的长度。

且我们假设 Nx>NhN_x \gt N_hNx​>Nh​。

对于下标 n<0n < 0n<0 或者 n>=Nxn >=N_xn>=Nx​，信号都为 0。

反转的滤波器系数

在实际的实时音频场景中，如虚拟现实、计算机游戏或数字音频工作站，我们知道 h[n]h[n]h[n] 但不知道 x[n]x[n]x[n]

因此，我们可以对 h[n]h[n]h[n] 进行反转，我们定义长度为 NhN_hNh​ 的 信号 ccc 为：
c[n]=h[Nh−n−1],n=0,…,Nh−1(2)c[n]=h\left[N_{h}-n-1\right], \quad n=0, \ldots, N_{h}-1 \tag{2} c[n]=h[Nh​−n−1],n=0,…,Nh​−1(2)

实际的卷积公式

引入上述两个假设后，我们可以将公式（1）改写为：
y[n]=(x[n]∗h[n])[n]=∑k=0Nh−1x[n+k]c[k],n=0,…,Nx−1.(3)\begin{array}{c} y[n]=(x[n] * h[n])[n] \\ =\sum_{k=0}^{N_{h}-1} x[n+k] c[k], \quad n=0, \ldots, N_{x}-1 . \end{array} \tag{3} y[n]=(x[n]∗h[n])[n]=∑k=0Nh​−1​x[n+k]c[k],n=0,…,Nx​−1.​(3)

这个公式与相关公式非常相似，但请记住，它仍然是卷积，尽管写法不同。

在这个新的公式中，一个卷积输出只是向量 ccc 和 xxx 的内积。

注意公式（3）的卷积使用的是 “same” 模式，如果我们往 xxx 中追加Nh−1N_h-1Nh​−1个0的话，那就变成了 “full” 模式。

正如你将看到的，full 模式将大大简化我们的讨论。

卷积过程可视化

下图说明了卷积是如何计算的。

橙色的框框表明哪些元素进行相乘，然后将乘法的结果累加得到 y[0]y[0]y[0]

有了上面的假设和卷积格式，我们可以写出它的实现。

线性卷积最简单的实现

在我们用 SIMD 提高 FIR 滤波器的速度之前，我们需要从一个基线开始：一个非 SIMD 的实现。

这可以通过以下方式实现：

struct FilterInput {
// assume that these fields are correctly initializedconst float* x;  // input signal with (N_h-1) zeros appendedsize_t inputLength;   // N_xconst float* c;  // reversed filter coefficientssize_t filterLength;  // N_hfloat* y;  // output (filtered) signal; // pointer to preallocated, uninitialized memorysize_t outputLength; // should be N_x in our context
};float* applyFirFilterSingle(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;for (auto i = 0u; i < input.outputLength; ++i) {y[i] = 0.f;for (auto j = 0u; j < input.filterLength; ++j) {y[i] += x[i + j] * c[j];}}return y;
}

正如你所看到的，这段代码的效率不高；我们一个一个地迭代样本。

它的算法时间复杂度是 O(NhNx)O(N_hN_x)O(Nh​Nx​)，让我们来看看如何对这份代码进行向量化。

循环向量化

FIR 滤波器场景下，有 3 中类型的循环向量化方法：

1. Inner loop vectorization (VIL),
2. Outer loop vectorization (VOL),
3. Outer and inner loop vectorization (VOIL).

它们的名字指明了我们将在哪一行将数据加载到 SIMD 寄存器中。其中最容易理解是 VIL

Inner loop vectorization (VIL)

在 VIL 中，我们将在内循环中进行向量化操作。
我们首先以一种粗略的方式来重写之前的代码。我们假设我们的向量长度为 4，这将对应于可以容纳 4 个浮点的寄存器（例如ARM的Neon寄存器）。

float* applyFirFilterInnerLoopVectorization(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;for (auto i = 0u; i < input.outputLength; ++i) {y[i] = 0.f;// Note the increment by 4for (auto j = 0u; j < input.filterLength; j += 4) {y[i] += x[i + j] * c[j] + x[i + j + 1] * c[j + 1] +x[i + j + 2] * c[j + 2] + x[i + j + 3] * c[j + 3];}}return y;
}

上述代码中，在内循环的每一次迭代中，我们做两个 4 元素向量的内积。就这样，我们计算出公式（3）中卷积和的一部分。

请注意，我们假设传入的向量已经是零填充的，并且长度是4的倍数。

下图显示了 VIL 的情况

当然，这份代码并不比之前的代码更优化。它只是以向量形式重写了代码。但这种向量形式现在很容易用向量指令来实现。

这个实现在真正的SIMD代码中看起来如何呢？

VIL NEON 实现

下面我将使用 NEON intrinsic 函数来实现 VIL。NEON 相关教程请参考Neon intrinsics 简明教程。

#ifdef __aarch64__
float* applyFirFilterInnerLoopVectorizationARM(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;assert(input.inputLength % 4 == 0);assert(input.filterLength % 4 == 0);for (auto i = 0u; i < input.outputLength; ++i) {y[i] = 0.f;float32x4_t outChunk = vdupq_n_f32(0.0f);for (auto j = 0u; j < input.filterLength; j += 4) {float32x4_t xChunk = vld1q_f32(x + i + j);float32x4_t cChunk = vld1q_f32(c + j);float32x4_t temp = vmulq_f32(xChunk, cChunk);outChunk = vaddq_f32(outChunk, temp);}y[i] = vaddvq_f32(outChunk);}return y;
}
#endif

同样，我们假设传入的向量已经是零填充，并且长度是4的倍数.

其算法时间复杂度是 O(NxNh/4)O(N_xN_h/4)O(Nx​Nh​/4)。当然，在复杂性理论中，这和非向量化算法是一样的。但是请注意，在内循环中，我们的迭代次数减少了 4 倍。这是因为我们可以用单条 NEON 指令对 4 个浮点数的向量进行操作。

Outer Loop Vectorization (VOL)

VOL 是一个更疯狂的方法。在这种方法中，我们试图在一次外循环中一次性计算出一个输出的向量。同样地，我们先给出粗略的 VOL 代码实现

float* applyFirFilterOuterLoopVectorization(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;// Note the increment by 4for (auto i = 0u; i < input.outputLength; i += 4) {y[i] = 0.f;y[i + 1] = 0.f;y[i + 2] = 0.f;y[i + 3] = 0.f;for (auto j = 0u; j < input.filterLength; ++j) {y[i] += x[i + j] * c[j];y[i + 1] += x[i + j + 1] * c[j];y[i + 2] += x[i + j + 2] * c[j];y[i + 3] += x[i + j + 3] * c[j];}}return y;
}

同样，我们假设传入的向量已经是零填充的，并且长度是4的倍数。
下图显示了 VOL 的情况

我们不是在内循环中从卷积和中计算4个元素（如VIL），而是从4个卷积和中计算1个元素。在VIL中，我们的内循环迭代次数减少了4次，在VOL中，我们的外循环迭代次数减少了4次。

因此，VOL并不比VIL更优化。

这个代码现在很容易用SIMD指令来实现。

VOL NEON Implementation

下面代码说明了如何使用 NEON 实现 VOL

float* applyFirFilterOuterLoopVectorizationARM(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;// Note the increment by 4for (auto i = 0u; i < input.outputLength; i += 4) {float32x4_t yChunk{0.0f, 0.0f, 0.0f, 0.0f};for (auto j = 0u; j < input.filterLength; ++j) {float32x4_t xChunk = vld1q_f32(x + i + j);float32x4_t temp = vmulq_n_f32(xChunk, c[j]);yChunk = vaddq_f32(yChunk, temp);}// store to memoryvst1q_f32(y + i, yChunk);}return y;
}

这段代码应该比最初的代码快 4 倍。在实践中，由于额外的代码，速度提升会更小。

一个问题出现了：我们可以做得更好吗？是的，我们可以!

Outer and Inner Loop vectorization (VOIL)

真正的突破来自于我们结合两种类型的矢量化。

在这种方法中，我们在外循环中计算一个向量，在内循环中使用向量内积，粗略的代码如下：

float* applyFirFilterOuterInnerLoopVectorization(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;// Note the incrementfor (auto i = 0u; i < input.outputLength; i += 4) {y[i] = 0.f;y[i + 1] = 0.f;y[i + 2] = 0.f;y[i + 3] = 0.f;// Note the incrementfor (auto j = 0u; j < input.filterLength; j += 4) {y[i] += x[i + j] * c[j] +x[i + j + 1] * c[j + 1] +x[i + j + 2] * c[j + 2] +x[i + j + 3] * c[j + 3];y[i + 1] += x[i + j + 1] * c[j] +x[i + j + 2] * c[j + 1] +x[i + j + 3] * c[j + 2] +x[i + j + 4] * c[j + 3];y[i + 2] += x[i + j + 2] * c[j] +x[i + j + 3] * c[j + 1] +x[i + j + 4] * c[j + 2] +x[i + j + 5] * c[j + 3];y[i + 3] += x[i + j + 3] * c[j] +x[i + j + 4] * c[j + 1] +x[i + j + 5] * c[j + 2] +x[i + j + 6] * c[j + 3];}}return y;
}

为什么 VOIL 更优？

你可能会对自己说 “好吧，但这只是一个手动展开的循环! 为什么它更快？”。

这是因为 SIMD 指令同时使用多个寄存器，如VOIL情况下，为处理器提供了更多的空间，可以更快地执行它们。这与 VIL 或 VOL 情况下只使用一个或两个寄存器形成对比。

当我说 "同时 "时，我并不是指多线程。我的意思是保持对各种寄存器的引用。这使处理器能够最有效地处理事情。

数据对齐

VOIL有如此大的优化潜力的另一个原因是，我们可以使用对齐的加载/存储SIMD指令来实现它。如何实现？这将是下一篇文章的主题!

让我们看看如何用 NEON 指令实现 VOIL。

VOIL NEON Implementation

float* applyFirFilterOuterInnerLoopVectorizationARM(FilterInput& input) {const auto* x = input.x;const auto* c = input.c;auto* y = input.y;std::array outChunk{};for (auto i = 0u; i < input.outputLength; i += 4) {for(auto k = 0; k < 4; ++k){outChunk[k] = vdupq_n_f32(0.0f);}for (auto j = 0u; j < input.filterLength; j += 4) {float32x4_t cChunk = vld1q_f32(c + j);for(auto k = 0; k < 4; ++k){float32x4_t xChunk = vld1q_f32(x + i + j +k);float32x4_t temp = vmulq_f32(cChunk, xChunk);outChunk[k] = vaddq_f32(temp, outChunk[k]);}}for(auto k = 0; k < 4; ++k){y[i + k] = vaddvq_f32(outChunk[k]);}}return y;
}

总结

在这篇文章中，我们讨论了什么是FIR滤波器以及如何有效地实现它；要么选择快速卷积算法，要么使用现代处理器的单指令、多数据指令。当然，你可以两者兼得

我们重新定义了卷积和，以方便讨论和实现。

我们研究了使用一种叫做循环矢量化的技术来实现FIR滤波器。我们展示了 VIL、VOL 和 VOIL 的普通C语言实现，讨论了它们的可视化，并使用 NEON 指令集展示了它们的SIMD等价物。

最后，我们指出，我们可以用对齐的数据做得更好。这将在下一篇文章中讨论。

请看看下面的有用的参考资料。整个代码可以在我的GitHub仓库中找到。

和往常一样，如果你有任何问题，请随时在下面发表。

参考

• 【音频处理】Fast Convolution 快速卷积算法简介
• Convolution: The secret behind filtering