我已经在 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
最佳答案
首先,我不认为您的初始向量代码是忠实的翻译——也就是说,我认为它不符合 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
关于algorithm - 使用向量实现 FIR 滤波器,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/43558345/