使用向量实现 FIR 滤波器

Question

我在 Haskell 中实现了 FIR 滤波器。我对 FIR 滤波器了解不多，而且我的代码很大程度上基于现有的 C# 实现。因此，我有一种感觉，我的实现有太多的 C# 风格，而不是真正的 Haskell 风格。我想知道是否有更惯用的 Haskell 方式来实现我的代码。理想情况下，我很幸运能够实现算法的一些高阶函数（映射、过滤器、折叠等）的组合。

我的 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