用 NEON 实现高效的 FIR 滤波器
Posted 芥末的无奈
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了用 NEON 实现高效的 FIR 滤波器相关的知识,希望对你有一定的参考价值。
系列文章目录
文章目录
- 系列文章目录
- 写在前面
- 前言
- 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
优化方案
如果你想让任何软件尽可能快地执行,有两种方法可以实现。
- 最佳的算法
- 高效的执行
同样的原则也适用于 DSP 代码。要想在代码中拥有一个快速的有限脉冲响应(FIR)滤波器,你可以这样做:
- 使用一种时间复杂度更低的算法,例如通过傅里叶变换到频域后进行快速卷积,或者
- 或者,利用硬件和软件资源的优势来有效地实现时域卷积。这通常意味着使用 SIMD 指令对你的代码进行矢量化。
何时应该使用时域的滤波器?
你或许会这么问自己:
既然我们有一个快速卷积算法,为什么我们还需要一个时域的实现呢?
快速卷积的算法复杂度是 O ( N log N ) O(N\\log N) O(NlogN),其中 N 可能输入信号或滤波器的长度。而线性卷积的算法复杂度是 O ( N 2 ) O(N^2) O(N2)
首先,这意味着对于足够小的 N N N,线性卷积算法将比快速卷积算法更快。
其次,快速卷积在复数域上操作,而线性卷积总是使用实数。这实际上意味着,快速卷积需要处理的数据是线性卷积的两倍。
如果时域卷积可能比快速卷积更快这一事实让你感到困惑,请想想排序算法。一般来说,快排被认为是最快的排序算法。但是如果你看一下排序的实现,比如std::sort,你会发现它最初只使用quicksort。在将排序的容器划分为足够小的部分后,会使用插入排序对其中的元素进行排序。
在解释了这一点之后,我们可以研究如何在时域中有效地实现FIR滤波。
如何加快FIR滤波器的实现速度?
简而言之:通过 SIMD。
加快滤波器执行速度的最好方法是使用 SIMD 指令一次处理多个样本。为了实现这一点,我们需要重写线性卷积算法,使我们的代码在向量上运行。
这个过程被称为循环向量化。
循环向量化通常由编译器完成,但这种自动向量化的程度通常对实时 DSP 来说是不够的。相反,我们需要准确地指示编译器要做什么。
SIMD 指令在对对齐的数据进行操作时能达到最佳性能。因此,数据对齐是我们应该考虑的另一个因素。
总之,一个高效的FIR滤波器的实现要同时使用两种策略。
- 循环矢量化
- 数据对齐
现在我们将详细讨论这两种策略。
初步假设
线性卷积公式如下:
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 \\mathbbZ . \\tag1
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]不一定是如此。
我们用 N x N_x Nx 表示 x [ n ] x[n] x[n] 的长度,用 N h N_h Nh 表示 h [ n ] h[n] h[n] 的长度。
且我们假设 N x > N h N_x \\gt N_h Nx>Nh。
对于下标 n < 0 n < 0 n<0 或者 n > = N x n >=N_x n>=Nx,信号都为 0。
反转的滤波器系数
在实际的实时音频场景中,如虚拟现实、计算机游戏或数字音频工作站,我们知道 h [ n ] h[n] h[n] 但不知道 x [ n ] x[n] x[n]
因此,我们可以对
h
[
n
]
h[n]
h[n] 进行反转,我们定义长度为
N
h
N_h
Nh 的 信号
c
c
c 为:
c
[
n
]
=
h
[
N
h
−
n
−
1
]
,
n
=
0
,
…
,
N
h
−
1
(2)
c[n]=h\\left[N_h-n-1\\right], \\quad n=0, \\ldots, N_h-1 \\tag2
c[n]=h[Nh−n−1],n=0,…,Nh−1(2)
实际的卷积公式
引入上述两个假设后,我们可以将公式(1)改写为:
y
[
n
]
=
(
x
[
n
]
∗
h
[
n
]
)
[
n
]
=
∑
k
=
0
N
h
−
1
x
[
n
+
k
]
c
[
k
]
,
n
=
0
,
…
,
N
x
−
1.
(3)
\\beginarrayc 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 . \\endarray \\tag3
y[n]=(x[n]∗h[n])[n]=∑k=0Nh−1x[n+k]c[k],n=0,…,Nx−1.(3)
这个公式与相关公式非常相似,但请记住,它仍然是卷积,尽管写法不同。
在这个新的公式中,一个卷积输出只是向量 c c c 和 x x x 的内积。
注意公式(3)的卷积使用的是 “same” 模式,如果我们往 x x x 中追加 N h − 1 N_h-1 Nh−1个0的话,那就变成了 “full” 模式。
正如你将看到的,full 模式将大大简化我们的讨论。
卷积过程可视化
下图说明了卷积是如何计算的。
橙色的框框表明哪些元素进行相乘,然后将乘法的结果累加得到
y
[
0
]
y[0]
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 ) O(N_hN_x) O(NhNx),让我们来看看如何对这份代码进行向量化。
循环向量化
FIR 滤波器场景下,有 3 中类型的循环向量化方法:
- Inner loop vectorization (VIL),
- Outer loop vectorization (VOL),
- 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 滤波器的主要内容,如果未能解决你的问题,请参考以下文章