4

我一直在分析我几乎完成的项目，我看到大约四分之三的 CPU 时间花在这个 IIR 过滤器功能上（目前在目标硬件上大约在一秒钟内调用了数十万次）因此，在其他一切正常的情况下，我想知道它是否可以针对我的特定硬件和软件目标进行优化。我的目标只有 iPhone 4 和更新版本，只有 iOS 4.3 和更新版本，只有 LLVM 4.x。如果有收获，一点点不精确可能是可以的。

static float filter(const float a, const float *b, const float c, float *d, const int e, const float x)
{
    float return_value = 0;

    d[0] = x;
    d[1] = c * d[0] + a * d[1];

    int j;

    for (j = 2; j <= e; j++) {
        return_value += (d[j] += a * (d[j + 1] - d[j - 1])) * b[j];
    }

    for (j = e + 1; j > 1; j--) {
        d[j] = d[j - 1];
    }
    return (return_value);
}

任何有关加速它的建议都值得赞赏，如果可以在默认编译器优化之外进行优化，也对您的意见感兴趣。我想知道 NEON SIMD 是否会有所帮助（这对我来说是新的领域），或者是否可以利用 VFP，或者 LLVM 自动矢量化是否会有所帮助。

我尝试了以下 LLVM 标志：

-ffast-math（没有显着差异）

-O4（在 iPhone 4S 上有很大的不同，时间减少了 25%，但在我的最低目标设备 iPhone 4 上没有显着差异，改进是我的主要目标）

-O3 -mllvm -unroll-allow-partial -mllvm -unroll-runtime -funsafe-math-optimizations -ffast-math -mllvm -vectorize -mllvm -bb-vectorize-aligned-only （来自 Hal Finkel 幻灯片的 LLVM 自动矢量化标志：http://llvm.org/devmtg/2012-04-12/Slides/Hal_Finkel.pdf，比 Xcode 发布目标的默认 LLVM 优化慢）

对其他标志、不同方法和功能更改开放。我宁愿单独留下输入和返回类型和值。实际上在这里讨论了使用 FIR 的 NEON 内在函数：https ://pixhawk.ethz.ch/_media/software/optimization/neon_support_in_the_arm_compiler.pdf但我没有足够的经验来成功应用信息我自己的情况。谢谢你的澄清。

编辑我很抱歉没有早点注意到这一点。在调查了 aka.nice 的建议后，我注意到为 e、a 和 c 传入的值始终是相同的值，并且我在运行前就知道它们，因此可以选择包含此信息的方法。

4

3 回答 3

7

以下是可以对代码进行的一些转换以使用 vDSP 例程。这些转换使用名为 T0、T1 和 T2 的各种临时缓冲区。每一个都是一个浮点数组，有足够的空间容纳 e-1 元素。

首先，使用临时缓冲区来计算a * b[j]. 这改变了原始代码：

for (j = 2; j <= e; j++) {
    return_value += (d[j] += a * (d[j + 1] - d[j - 1])) * b[j];
}

至：

vDSP_vsmul(b+2, 1, &a, T0, 1, e-1);
for (j = 2; j <= e; j++)
    return_value += (d[j] += (d[j+1] - d[j-1])) * T0[j-2];

然后使用 vDSP_vmul 计算d[j+1] * T0[j-2]：

vDSP_vsmul(b+2, 1, &a, T0, 1, e-1);
vDSP_vmul(d+3, 1, T0, 1, T1, 1, e-1);
for (j = 2; j <= e; j++)
    return_value += (d[j] += T1[j-2] - d[j-1] * T0[j-2];

接下来，将 vDSP_vmul 提升为 vDSP_vma（向量乘加）以计算d[j] + d[j+1] * T0[j-2]：

vDSP_vsmul(b+2, 1, &a, T0, 1, e-1);
vDSP_vma(d+3, 1, T0, 1, d+2, 1, T1, 1, e-1);
for (j = 2; j <= e; j++)
    return_value += (d[j] = T1[j-2] - d[j-1] * T0[j-2];

我想我会计时，看看是否有任何改进。有一些问题：

  • SIMD 代码在数据为 16 字节对齐时效果最佳。使用数组索引，例如j-1和j+1防止这种情况。手机中的 ARM 处理器在数据未对齐方面并不像其他一些处理器那样糟糕，但性能会因型号而异。
  • 如果 e 很大（超过几千个），那么在 vDSP_vma 操作期间 T0 和 d 可能会从缓存中清除，并且下面的循环将不得不重新加载它们。有一种称为露天采矿的技术可以减少这种影响。我现在不会详细说明它，但本质上，操作被划分为更小的数组条带。
  • 最后循环中的 IIR 可能仍会成为处理器的瓶颈。vDSP 中有用于执行某些 IIR 的例程（例如 vDSP_deq22），但尚不清楚此过滤器是否可以以与 vDSP 例程充分匹配的方式表示，以获得比转换可能损失的性能更多的性能.
  • 最后循环中用于计算 return_value 的求和也可以从循环中删除并替换为 vDSP 例程（可能是 vDSP_sve），但我怀疑由 IIR 引起的松弛将允许在不增加显着执行时间的情况下完成添加环形。

以上是我的头顶；我没有测试过代码。我建议逐一进行转换，以便您可以在每次更改后测试代码并在继续之前识别任何错误。

如果您能找到不是 IIR 的令人满意的滤波器，则可能会有更多的性能优化。

于 2012-10-15T19:06:22.830 回答
0

我宁愿单独留下输入和返回类型和值……</p>

尽管如此，将渲染从浮点数转换为整数会有很大帮助。

将更改本地化到您提供的实现将没有用。但是，如果您将其扩展为仅将 FIR 重新实现为整数，它可以很快得到回报（除非保证大小总是非常小——然后转换/移动时间成本更高）。当然，将渲染图的较大部分移动到整数将带来更大的收益，并且需要更少的转换。

另一个考虑因素是查看 Accelerate.framework 中的实用程序（可能使您免于编写自己的 asm）。

于 2012-10-15T18:43:56.840 回答
0

我尝试了这个小练习来用延迟运算符 z 重写你的过滤器

例如，对于 e=4，我将输入 u 和输出 y 重命名

d1*z= u
d2*z= c*u + a*d1
d3*z= d2 + a*(d3-d1*z)
d4*z= d3 + a*(d4-d2*z)
d5*z= d4 + a*(d5-d3*z)
y = (b2*d3*z + b3*d4*z + b4*d5*z)

请注意，di 是过滤器状态。
d3*z 是 d3 的下一个值（它在您的代码中似乎是变量 d2）
然后您可以消除 di 以在 z 中编写传递函数 y/u。
然后，您会发现通过分解/简化上述传递函数，最小表示只需要 e 个状态。
分母是z*(z-a)^3，即 0 处的一个极点，另一个具有多重性 (e-1) 的 a 处。

然后，您可以将过滤器置于标准状态空间矩阵表示中：

z*X = A*X + B*u
y = C*X + d*u

使用特定形式的极点，您可以分解部分分数展开中的传递并以这种特殊形式获得矩阵 A 和 B（类似于 matlab 的符号）

A = [0 1 0 0;    B=[0;
     0 a 1 0;       0;
     0 0 a 1;       0;
     0 0 0 a]       1]

不过，C 和 d 不太容易......
它们是从分子和部分分数展开的直接项中提取的
它们是 bi、c（1 度）和 a（e 度）中的多项式
对于 e=4，我有

C=[   a^3*b2            - a^2*b3                      + a*b4 ,
     -a^2*b2              + a*b3                + (c-a^2)*b4 ,
        a*b2        + (c-a^2)*b3 + (-2*a^2*c-2*a^2-a+a^4)*b4 ,
  (c-a^2)*b2 + (-a^2-a^2*c-a)*b3 +         (-2*a*c+2*a^3)*b4 ]
d=     -a*b2            - a*c*b3                    + a^2*b4

如果你能在 e 中找到控制 C & d 的递归，并预先计算它们，
那么过滤器可以简化为那些简单的向量操作：

z*X = a*[ 0; x2 ; x3 ; x4 ... xe ] + [x2 ; x3 ; x4 ... xe ; u ];
Y = C*[ x1 ; x2 ; x3 ; x4 ... xe ] + d*u

或者表示为函数(Xnext,y)=filter(X,u,a,C,d,e)伪代码：

y = dot_product( C , X) + d*u; // (like BLAS _DOT)
Xnext(1:e-1) = X(2:e); // this is a memcopy (like BLAS _COPY)
Xnext(e)=u;
X(1)=0;
Xnext=a*X+Xnext; // this is an inplace vector muladd (like BLAS _AXPY)

X=Xnext; // another memcopy outside the function (can be moved inside).

请注意，如果您使用 BLAS 函数，您的代码将可移植到许多硬件，而不仅仅是以 Apple 为中心，我猜性能不会有太大不同。

编辑：关于部分分数扩展

纯部分分数展开将给出对角状态空间表示，以及一个充满 1 的矩阵 B。这也可能是一个有趣的变体。（并行过滤器）
我上面使用的变体更像是级联或梯子（串联过滤器）。

于 2012-10-16T22:48:50.593 回答