(In response to Tom Hawkins' posting of an IIR filter in Atom) We're still experimenting with how to best describe streaming computations with feedback in Feldspar. But for completeness, here one possible implementation of an IIR filter:
iir :: forall m n o a . (NaturalT m, NaturalT n, NaturalT o, Num a , Primitive a) => VectorP m a -> VectorP n a -> VectorP o a -> VectorP o a
iir as bs = feedback f where f :: VectorP o a -> VectorP o a -> Data a f inPrev outPrev = dotProd as (resize inPrev) - dotProd bs (resize outPrev)
(Please don't mind the type clutter -- we hope to get rid of most of it in the future.) The local function `f` computes a single output, and the `feedback` combinator applies `f` across the input stream. You can find the resulting C code attached. As you can see, the generated C has lots of room for optimization, but the time complexity is right (one top-level loop with two inner loops in sequence). We plan to tackle the more small-scale optimizations in the future. The dot product is defined in standard Haskell style:
dotProd :: (Num a, Primitive a) => VectorP n a -> VectorP n a -> Data a dotProd as bs = fold (+) 0 (zipWith (*) as bs)
Interestingly, `feedback` is also defined within Feldspar:
feedback :: forall n a . (NaturalT n, Storable a) => (VectorP n a -> VectorP n a -> Data a) -> VectorP n a -> VectorP n a
feedback f inp = unfreezeVector (length inp) outArr' where outArr :: Data (n :> a) outArr = array []
outArr' = for 0 (length inp - 1) outArr $ \i arr -> let prevInps = reverse $ take (i+1) inp prevOutps = reverse $ take i $ unfreezeVector i arr a = f prevInps prevOutps in setIx arr i a
This definition uses low-level data structures and loops, and this is not something that ordinary Feldspar users should write. It is our hope that a few combinators like this one can be defined once and for all, and then reused for a wide range of DSP applications. It turns out that FIR filters are much nicer :)
fir :: (NaturalT m, Num a , Primitive a) => VectorP m a -> VectorP n a -> VectorP n a
fir coeffs = map (dotProd coeffs . resize . reverse) . inits
C code attached. / Emil