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

double dir(int M, double* a,
           int L, double* b,
           double *w, double *v,
           double x) {

   v[0] = x;                      // current input sample
   w[0] = 0;                      // current output
   
   for (int i=0; i<=L; i++) {
      w[0] += b[i] * v[i];        // numerator
   }
   
   for (int i=0; i<=M; i++) {
      w[0] -= a[i] * w[i];        // denominator
   }
   
   for (int i=L; i>=1; i--) {
      v[i] = v[i-1];              // shift numerator memory
   }
   
   for (int i=M; i>=1; i--) {
      w[i] = w[i-1];              // shift denominator memory
   }
   
   return w[0];                   // current output sample
}