用 NEON 实现高效的 FIR 滤波器

Posted 2023-02-15 芥末的无奈

系列文章目录

• 数字信号处理中的 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]$ 定义，FIR 滤波器的输出 $y[n]$ 是输入信号 $x[n]$ 与脉冲响应的卷积。我们把它写成：
$$y[n]=x[n]\ast h[n]=\sum _{k=-\infty }^{\infty }x[k]h[n-k]$$

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

优化方案

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

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

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

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

---

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

你或许会这么问自己：

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

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

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

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

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

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

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

简而言之：通过 SIMD。

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

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

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

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

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

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

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

初步假设

线性卷积公式如下：
$$x[n]\ast h[n]=\sum _{k=-\infty }^{\infty }x[k]h[n-k]=y[n],n\in \mathbb{Z}.\text{\hspace{1em}(1)}$$

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

有限长度信号

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

我们用 $N_x$ 表示 $x[n]$ 的长度，用 $N_h$ 表示 $h[n]$ 的长度。

且我们假设 $N_x>N_h$。

对于下标 $n<0$ 或者 $n>=N_x$，信号都为 0。

反转的滤波器系数

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

因此，我们可以对 $h[n]$ 进行反转，我们定义长度为 $N_h$ 的 信号 $c$ 为：
$$c[n]=h\left[N_h-n-1\right],n=0,\dots ,N_h-1\text{\hspace{1em}(2)}$$

实际的卷积公式

引入上述两个假设后，我们可以将公式（1）改写为：
$$\begin{array}{c}y[n]=(x[n]\ast h[n])[n]\\ ={\sum }_{k=0}^{N_h-1}x[n+k]c[k],n=0,\dots ,N_x-1.\end{array}\text{\hspace{1em}(3)}$$

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

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

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

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

卷积过程可视化

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

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

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

线性卷积最简单的实现

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

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

struct FilterInput 
// assume that these fields are correctly initialized
  const float* x;  // input signal with (N_h-1) zeros appended
  size_t inputLength;   // N_x
  const float* c;  // reversed filter coefficients
  size_t filterLength;  // N_h
  float* y;  // output (filtered) signal; 
             // pointer to preallocated, uninitialized memory
  size_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(N_h{N}_x)$，让我们来看看如何对这份代码进行向量化。

循环向量化

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 4
    for (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以上是关于用 NEON 实现高效的 FIR 滤波器的主要内容，如果未能解决你的问题，请参考以下文章
 
 NEON Intrinsics 练习题
 NEON Intrinsics 练习题
 语音信号的“短时时域”分析
 使用 FIR 滤波器过滤光谱
 FIR IP核使用
 算法常用算法概览