ベクトルを使用したFIRフィルタの実装

Question

私はHaskellでFIRフィルタを実装しました。私はそれほどFIRフィルターについて知らず、私のコードは既存のC＃実装に大きく基づいています。したがって、私の実装はC＃スタイルがあまりにも多く、実際にはハスケルのようなものではないと感じています。私は、コードを実装するためにもっと慣用的な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となります。あなたのHaskellコードで

（それはbiの他のすべての値についてcontinueを打つだろう）、しかし、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