/* dir.c - IIR filtering in direct form */

double dir(M, a, L, b, w, v, x)           /* usage: y = dir(M, a, L, b, w, v, x); */
double *a, *b, *w, *v, x;                 /* \(v,w\) are internal states */
int M, L;                                 /* denominator and numerator orders */
{
       int i;

       v[0] = x;                          /* current input sample */
       w[0] = 0;                          /* current output to be computed */

       for (i=0; i<=L; i++)               /* numerator part */
              w[0] += b[i] * v[i];

       for (i=1; i<=M; i++)               /* denominator part */
              w[0] -= a[i] * w[i];

       for (i=L; i>=1; i--)               /* reverse-order updating of \(v\) */
              v[i] = v[i-1];

       for (i=M; i>=1; i--)               /* reverse-order updating of \(w\) */
              w[i] = w[i-1];

       return w[0];                       /* current output sample */
}