Emil Axelsson emax at chalmers.se
Wed Nov 25 08:29:39 EST 2009

```(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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: fir.c
Type: text/x-csrc
Size: 1520 bytes
Desc: not available