使用矢量實現FIR濾波器

我在Haskell中實現了FIR濾波器。我不太瞭解FIR濾波器，我的代碼主要基於現有的C＃實現。因此，我感覺我的實現有太多的C＃風格，並不像Haskell那樣。我想知道是否有更實用的Haskell實現我的代碼的方式。理想情況下，我很幸運可以實現算法的高階函數（地圖，過濾器，摺疊等）的組合。使用矢量實現FIR濾波器

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

來源

2017-04-22 Thomas Vanhelden

順便說一句，這種操作稱爲[卷積（https://en.wikipedia.org/wiki/Convolution）。在_O_（_n_²）窗體中，如果兩個輸入數組都很大，則無論語言如何，這個窗體總是很慢。因此，行業標準使用基於[FFT]（https://en.wikipedia.org/wiki/Fast_Fourier_transform）的_fast卷積，其僅使用_O_（_n_·log_n_）。如果其中一個陣列很小，Daniel Martin的版本應該沒問題。 – leftaroundabout

是的，你是對的！我會考慮使用FFT進行研究。 –

首先，我不認爲你的初始向量的代碼是一個忠實的翻譯 - 也就是說，我認爲它與C＃代碼不一致。例如，假設「x」和「b」（「b」是C＃中的coeff）的長度均爲3，並且所有值均爲1.0。然後，對於y[0]，C＃代碼將生成x[0] * coeff[0]或1.0。（這將打擊continue爲bi所有其他值）

隨着你的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

來源

2017-04-22 16:12:35

我正在考慮使用切片進行實現，但由於某種原因完全忘記了它。非常感謝這個版本（和更正）！它削減了8秒的時間。 –

使用矢量實現FIR濾波器

回答

相關問題