/*-------------- Telecommunications & Signal Processing Lab ---------------
                             McGill University

Routine:
  void FIfiltIIR (const float x[], float y[], int Nout, const float h[][],
		  int Nsec)

Purpose:
  Filter a signal using a cascade of biquadratic IIR filters

Description:
  This procedure forms the output of a cascade of IIR filters.  Each IIR filter
  is a biquadratic section with the z-transform,

             h(i,0)*z^2 + h(i,1)*z + h(i,2)
    H(i,z) = ------------------------------ .
                z^2 + h(i,3)*z + h(i,4)

  Each filter section needs an input sample, two past input samples and two
  past output samples to calculate each output sample.  In the cascade
  connection, the output of one section is the input to next section.  This
  means that past output samples for one section are also the past input
  samples for the next section.  For Nsec sections, the total filter "memory"
  required is mem=2*(Nsec+1) values.  In this routine, the input array has 2
  extra elements at the beginning to store the two past inputs.  The output
  array must be primed with 2 past outputs for each of the Nsec sections.

  On input, the first two elements of the input array (x[0] and x[1]) are the
  two previous inputs to the filter; the next Nout elements are new input
  values.  The first 2*Nsec elements of the output array y[] prime the
  calculation.  Elements y[0] and y[1] are the previous outputs of the filter.
  The next two elements are the outputs of the previous to last filter section,
  and so on up to y[2*Nsec-2] and y[2*Nsec-1] which are previous outputs of the
  first section.

  The filter is processed one section at a time.  For the first section, the
  input array is the input.  For the second section, the output from the first
  section is the input and so on.  The output of these sections overlay each
  other.  The computation is arranged such that the overall output of the
  filter appears as Nout samples starting with y[2].  On return, the last
  2*Nsec elements of this array are suitable for priming the array for the
  next invocation.  This priming can be accomplished by moving the top 2*Nsec
  elements of the output array to the bottom of the array.  The input array can
  be handled similarly; the top two elements of x[], serve as the needed memory
  values for the next invocation of this routine.

  The input and output arrays can also be shared.  As output values are
  computed, they can overlay the input values that are no longer needed.
  Consider an array x[] with mem+Nout elements, where mem=2*(Nsec+1).  The
  first 2 elements are previous outputs of the overall filter, the next two
  elements are previous outputs of the penultimate section, and so on up to
  elements x[2*Nsec] and x[2*Nsec+1] which are previous inputs.  The remaining
  Nout values are new input values.  Invoke this routine as
    FIfiltIIR (&x[mem-2], x, Nout, h, Nsec)
  On return, the Nout new output values start at x[2].  The top mem elements
  of x[] are suitable for priming the array for the next invocation.

Parameters:
   -> const float x[]
      Input array with Nout+ elements.  The first two elements are previous
      input values.  The next Nout values are new input values.
  <-> float y[]
      Input/output array with Nout+2*Nsec elements.  On input, the first 2*Nsec
      elements represent filter memory values.  On output, Nout outputs are
      placed starting at y[2].  The top 2*Nsec elements of the array contain
      new values which can be used as the filter memory for the next invocation
      of this routine.
   -> int Nout
      Number of output samples
   -> const float h[][5]
      Array of filter coefficients for the filter sections.  Each filter
      section is defined by 5 filter coefficients.
   -> int Nsec
      Number of filter sections

Author / revision:
  P. Kabal  Copyright (C) 1995
  $Revision: 1.8 $  $Date: 1995/06/08 05:45:41 $

-------------------------------------------------------------------------*/

static char rcsid[] = "$Id: FIfiltIIR.c 1.8 1995/06/08 AFsp-V2R1 $";

#include <libtsp.h>

void
FIfiltIIR (x, y, Nout, h, Nsec)

     const float x[];
     float y[];
     int Nout;
     const float h[][5];
     int Nsec;

{
  int k;
  int i;
  float *xp;

/* Process the input, one section at a time */
  k = 2 * Nsec;
  xp = (float *) &x[0];
  for (i = 0; i < Nsec; ++i) {
    FIbiquad (xp, &y[k-2], Nout, h[i]);
    k = k - 2;
    xp = &y[k];
  }
}