(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
#include "feldspar.h"
void fir( signed int var0_0_0, signed int var0_0_1[10], signed int var0_1_0, signed int var0_1_1[100], signed int *out_0, signed int out_1[100] )
{
signed int var23[100];
{
int var1;
for( var1 = 0; var1 < var0_1_0; var1 += 1)
{
signed int var7;
int var8;
signed int var9;
int var10;
signed int var11;
signed int var12;
signed int var17;
signed int var22_0;
var7 = (var1 + 1);
var8 = (var0_1_0 <= var7);
if(var8)
{
var9 = var0_1_0;
}
else
{
var9 = var7;
}
var10 = (var0_0_0 <= var9);
if(var10)
{
var11 = var0_0_0;
}
else
{
var11 = var9;
}
var12 = (var11 - 1);
var17 = (var9 - 1);
var22_0 = 0;
var23[var1] = 0;
{
int var13;
var13 = (var22_0 <= var12);
while(var13)
{
var23[var1] = (var23[var1] + (var0_0_1[var22_0] * var0_1_1[(var17 - var22_0)]));
var22_0 = (var22_0 + 1);
var13 = (var22_0 <= var12);
}
}
}
}
*out_0 = var0_1_0;
copy_arrayOf_signed_int(&(var23[0]), 100, &(out_1[0]));
}
#include "feldspar.h"
void iir( signed int var0_0_0, signed int var0_0_1[10], signed int var0_1_0, signed int var0_1_1[10], signed int var0_2_0, signed int var0_2_1[100], signed int *out_0, signed int out_1[100] )
{
signed int var3;
signed int var51_0;
signed int var51_1[100];
signed int var53[100];
var3 = (var0_2_0 - 1);
var51_0 = 0;
copy_arrayOf_signed_int(&({}[0]), 100, &(var51_1[0]));
{
int var4;
var4 = (var51_0 <= var3);
while(var4)
{
signed int var12;
int var13;
signed int var14;
int var15;
signed int var16;
signed int var17;
signed int var22;
signed int var27_0;
signed int var27_1;
int var33;
signed int var34;
int var35;
signed int var36;
signed int var37;
signed int var42;
signed int var47_0;
signed int var47_1;
signed int var49[100];
var12 = (var51_0 + 1);
var13 = (var0_2_0 <= var12);
if(var13)
{
var14 = var0_2_0;
}
else
{
var14 = var12;
}
var15 = (var0_0_0 <= var14);
if(var15)
{
var16 = var0_0_0;
}
else
{
var16 = var14;
}
var17 = (var16 - 1);
var22 = (var14 - 1);
var27_0 = 0;
var27_1 = 0;
{
int var18;
var18 = (var27_0 <= var17);
while(var18)
{
var27_1 = (var27_1 + (var0_0_1[var27_0] * var0_2_1[(var22 - var27_0)]));
var27_0 = (var27_0 + 1);
var18 = (var27_0 <= var17);
}
}
var33 = (var51_0 <= var51_0);
if(var33)
{
var34 = var51_0;
}
else
{
var34 = var51_0;
}
var35 = (var0_1_0 <= var34);
if(var35)
{
var36 = var0_1_0;
}
else
{
var36 = var34;
}
var37 = (var36 - 1);
var42 = (var34 - 1);
var47_0 = 0;
var47_1 = 0;
{
int var38;
var38 = (var47_0 <= var37);
while(var38)
{
var47_1 = (var47_1 + (var0_1_1[var47_0] * var51_1[(var42 - var47_0)]));
var47_0 = (var47_0 + 1);
var38 = (var47_0 <= var37);
}
}
copy_arrayOf_signed_int(&(var51_1[0]), 100, &(var49[0]));
var49[var51_0] = (var27_1 - var47_1);
var51_0 = (var51_0 + 1);
copy_arrayOf_signed_int(&(var49[0]), 100, &(var51_1[0]));
var4 = (var51_0 <= var3);
}
}
{
int var1;
for( var1 = 0; var1 < var0_2_0; var1 += 1)
{
var53[var1] = var51_1[var1];
}
}
*out_0 = var0_2_0;
copy_arrayOf_signed_int(&(var53[0]), 100, &(out_1[0]));
}
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
