/* can.c - IIR Filtering in canonical form */

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

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