algorithm - 使用向量实现 FIR 滤波器

Question

我已经在 Haskell 中实现了一个 FIR 滤波器。我不太了解 FIR 滤波器，我的代码很大程度上基于现有的 C# 实现。因此，我觉得我的实现有太多的 C# 风格，而不是真正的 Haskell 风格。我想知道是否有更惯用的 Haskell 方式来实现我的代码。理想情况下，我很幸运能够实现该算法的高阶函数(map、filter、fold 等)的某种组合。

我的 Haskell 代码如下所示:

  applyFIR :: Vector Double -> Vector Double -> Vector Double
  applyFIR b x = generate (U.length x) help
      where
        help i = if i >= (U.length b - 1) then loop i (U.length b - 1) else 0
        loop yi bi = if bi < 0 then 0 else b !! bi * x !! (yi-bi) + loop yi (bi-1)
        vec !! i = unsafeIndex vec i -- Shorthand for unsafeIndex

此代码基于以下 C# 代码:

public float[] RunFilter(double[] x)
      {
         int M = coeff.Length;
         int n = x.Length;
         //y[n]=b0x[n]+b1x[n-1]+....bmx[n-M]
         var y = new float[n];
         for (int yi = 0; yi < n; yi++)
         {
            double t = 0.0f;
            for (int bi = M - 1; bi >= 0; bi--)
            {
               if (yi - bi < 0) continue;

               t += coeff[bi] * x[yi - bi];
            }
            y[yi] = (float) t;
         }

         return y;
      }

如您所见，它几乎是一个直接的副本。如何将我的实现变成更像 Haskell 的实现？你有什么想法？我唯一能想到的就是使用 Vector.generate。

我知道 DSP 库有可用的实现。但它使用列表并且对于我的用例来说太慢了。此 Vector 实现比 DSP 中的实现快很多。

我还尝试使用 Repa 实现该算法。它比 Vector 实现更快。这是结果:

applyFIR :: V.Vector Float -> Array U DIM1 Float -> Array D DIM1 Float
applyFIR b x = R.traverse x id (\_ (Z :. i) -> if i >= len then loop i (len - 1) else 0)
  where
    len = V.length b
    loop :: Int -> Int -> Float
    loop yi bi = if bi < 0 then 0 else (V.unsafeIndex b bi) * x !! (Z :. (yi-bi)) + loop yi (bi-1)
    arr !! i = unsafeIndex arr i

Answer 1

首先，我不认为您的初始向量代码是忠实的翻译——也就是说，我认为它不符合 C# 代码。例如，假设“x”和“b”(“b”是 C# 中的 coeff)的长度均为 3，并且所有值为 1.0。然后对于 y[0]，C# 代码将生成 x[0] * coeff[0]，或 1.0。 (对于 bi 的所有其他值，它将点击 continue)

然而，对于您的 Haskell 代码，help 0 会生成 0。您的 Repa 版本似乎遇到了同样的问题。

那么让我们从一个更忠实的翻译开始:

applyFIR :: Vector Double -> Vector Double -> Vector Double
applyFIR b x = generate (U.length x) help
    where
      help i = loop i (min i $ U.length b - 1)
      loop yi bi = if bi < 0 then 0 else b !! bi * x !! (yi-bi) + loop yi (bi-1)
      vec !! i = unsafeIndex vec i -- Shorthand for unsafeIndex

现在，您基本上是在进行这样的计算，例如 y[3]:

  ... b[3]   |   b[2]   |   b[1]   |   b[0]
      x[0]   |   x[1]   |   x[2]   |   x[3]   |   x[4]   |   x[5]   | ....
           multiply
    b[3]*x[0]|b[2]*x[1] |b[1]*x[2] |b[0]*x[3] 
           sum
      y[3] = b[3]*x[0] + b[2]*x[1] + b[1]*x[2] + b[0]*x[3]

因此，考虑您正在做的事情的一种方法是“获取 b 向量，将其反转，并计算结果的 i 点，第 行>b[0] 与 x[i] 相乘，将所有对应的 x 和 b 项相乘，并计算总和".

那么让我们这样做:

applyFIR :: Vector Double -> Vector Double -> Vector Double
applyFIR b x = generate (U.length x) help
  where
    revB = U.reverse b
    bLen = U.length b
    help i = let sliceLen = min (i+1) bLen
                 bSlice = U.slice (bLen - sliceLen) sliceLen revB
                 xSlice = U.slice (i + 1 - sliceLen) sliceLen x
             in U.sum $ U.zipWith (*) bSlice xSlice