.. ECE 4703 

.. _lecture 3:

IIR Implementation Techniques 
=============================

.. contents:: 

The purpose of this lecture is as follows.

* To describe the direct-form I and direct-form II implementations of IIR designs
* To describe the cascade-form implementations of IIR designs
* To describe transpose-form implementations of IIR designs
* To describe parallel-form implementations of IIR designs
* To demonstrate useful functions in Matlab used for polynomial manipulation

Direct-form IIR implementation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

An IIR filter is a digital filter described by the following input/output relationship

.. math::

   y[n] = \sum_{k=0}^{M} b[k] x[n-k] - \sum_{k=1}^{N} a[l] y[n-l]

which means that the output of an IIR filter is defined by both the previous
inputs as well as the previous outputs.

We can write this equivalently as a z-transform

.. math::

   Y(z) = X(z) . \sum_{k=0}^{M} b[k] z^{-k} - Y(z) . \sum_{k=1}^{N} a[k] z^{-k}

Assume that the order of the numerator and denominator are identical, then we can write
the filter transfer function as follows.

.. math::

   G(z) =& \frac{Y(z)}{X(z)}  \\
   G(z) =& \frac{\sum_{k=0}^{N} b[k] z^{-k}}{ 1 + \sum_{k=1}^{N} a[k] z^{-k}}

Since these are N-th order polynomials in z, they contribute N zeroes and N poles.
For an IIR system to be causal and stable, all poles must lie within the unit circle.

We can find the locations of the zeroes and poles of an arbitrary :math:`G(z)` by 
multiplying the numerator and denominator by :math:`z^N` and finding the roots
of each polynomial. See further for an example on using Matlab to find these roots.

Direct-form I IIR Implementation
""""""""""""""""""""""""""""""""

If we consider a transfer function of the form

.. math::

   G(z) = \frac{B(z)}{A(z)}

and we look for the output signal :math:`Y(z)` in response to the input signal :math:`X(z)`, then we compute:

.. math::

   Y(z) &= G(z).X(z)  \\
        &= \frac{B(z)}{A(z)}.X(z)


A **direct-form I IIR** implementation first computes the zeroes, and then the poles:

.. math::

   Y(z) &= \frac{1}{A(z)}.[B(z).X(z)]

In the time domain we would implement

.. math::

   w(n) &= \sum_{k=0}^{M} b[k] x[n-k]  \\
   y(n) &= w(n) - \sum_{k=1}^{N} a[l] y[n-l]


This formulation leads to a design as shown below. This design requires
:math:`1 + M + N` multiplications per output sample, 
:math:`M + N` additions per output sample,
and :math:`M + N` delays.

.. figure:: images/directformi.jpg
   :figwidth: 500px
   :align: center

Direct-form II IIR Implementation
"""""""""""""""""""""""""""""""""

A **direct-form II IIR** implementation first computes the poles, and then the zeroes:

.. math::

   Y(z) &= B(z) . [\frac{1}{A(z)}.X(z)]

This structure is created by swapping the order of the operations in a direct-form I
structure. However, an important optimization is possible because the taps can be
shared by the feedback part as well as by the feedforward part.

.. figure:: images/directformii.jpg
   :figwidth: 400px
   :align: center

This design requires :math:`1 + M + N` multiplications per output sample, 
:math:`M + N` additions per output sample,
and :math:`max(M,N)` delays.

This implementation has the mimimum number of delays possible for a given transfer function.
Hence, an N-th order FIR, and an N-th order IIR, both can be implemented with only N delays.

Cascade IIR Implementation
""""""""""""""""""""""""""

Just as with the FIR implementation, a cascade IIR is created by factoring the
nominator and denominator polynomials in :math:`G(z)`. Assuming both
nominator and denominator have order N, then

.. math::

   G(z) =& \frac{\sum_{k=0}^{N} b[k] z^{-k}}{ 1 + \sum_{k=1}^{N} a[k] z^{-k}} \\
        =& A . \prod_{k=1}^{N} \frac{1 - \beta_k . z^{-1}}{1 - \alpha_k . z^{-1}}


It's common to formulate this as *second-order* sections, in order to ensure
that complex (conjugate) poles and complex (conjugate) zeroes can be represented
with real coefficients. Each such a section contains a conjugate pole pair and/or conjugate zero pair.

.. math::

   G_k(z) =& \frac{(1 - \beta_k . z^{-1})(1 - \beta_k^* . z^{-1})}{(1 - \alpha_k . z^{-1})(1 - \alpha_k^* . z^{-1})} \\
   G_k(z) =& \frac{1 - 2.Re(\beta_k).z^{-1} + |\beta_k|^2.z^{-2}}{1 - 2.Re(\alpha_k).z^{-1} + |\alpha_k|^2.z^{-2}}

.. figure:: images/cascadeiir.jpg
   :figwidth: 600px
   :align: center

Cascade IIR Example
"""""""""""""""""""

We implement the following cascade IIR design in C.
It has four poles, located at :math:`\pm 0.5 . e^{\pm j \pi/4}`, as well
as three zeroes, located at :math:`j, -1, -j`.

.. figure:: images/cascadeiirexample.jpg
   :figwidth: 300px
   :align: center

The first step is to derive the proper filter coefficients. We group the poles and zeroes 
as follows.

.. math::

   C1(z) &= \frac{(1 - j.z^{-1}).(1 + j.z^{-1})}{(1 - e^{j.\pi/4}.z^{-1}/2).(1 - e^{-j.\pi/4}.z^{-1}/2)} \\
   C2(z) &= \frac{(1  + z^{-1})}{(1 + e^{j.\pi/4}.z^{-1}/2).(1 + e^{-j.\pi/4}.z^{-1}/2)} 

Which can be multiplied out to derive the filter coefficients:

.. math::

   C1(z) &= \frac{(1 + z^{-2})}{   (1 - \frac{1}{\sqrt 2}.z^{-1} + \frac{1}{4}.z^{-2})} \\
   C2(z) &= \frac{(1  + z^{-1})}{   (1 + \frac{1}{\sqrt 2}.z^{-1} + \frac{1}{4}.z^{-2})}


The construction of the filter proceeds similar to the cascade FIR filter. 
First, we create a data structure for the cascade IIR stage, containing
filter coefficients and filter state. We will use a direct-form Type II design, which
gives us a state of only two variables.

.. code:: c
   :number-lines: 1

   typedef struct cascadestate {
       float32_t s[2];   // state
       float32_t b[3];  // nominator coeff  b0 b1 b2
       float32_t a[2];  // denominator coeff   a1 a2
   } cascadestate_t;
   
   float32_t cascadeiir(float32_t x, cascadestate_t *p) {
       float32_t v = x - (p->s[0] * p->a[0]) -  (p->s[1] * p->a[1]);
       float32_t y = (v * p->b[0]) + (p->s[0] * p->b[1]) + (p->s[1] * p->b[2]);
       p->s[1] = p->s[0];
       p->s[0] = v;
       return y;
   }
   
   void createcascade(float32_t b0,
                      float32_t b1,
                      float32_t b2,
                      float32_t a1,
                      float32_t a2,
                      cascadestate_t *p) {
       p->b[0] = b0;
       p->b[1] = b1;
       p->b[2] = b2;
       p->a[0] = a1;
       p->a[1] = a2;
       p->s[0] = p->s[1] = 0.0f;
   }


The ``cascadestate_t`` now contains both feed-forward and feed-back coefficients.
The ``cascadeiir`` function is likely the most complicated element
for this design. The order of evaluation of each expression is chosen
so that we don't update filter state until it has been used for all expressions
required. In ``cascadeiir`` we first compute the intermediate variable ``v`` which
corresponds to the center-tap of the filter. Next, we compute the
feed-forward part based in this intermediate ``v``. Finally, we update the filter
state and return the filter output ``y``.

The filter program then consists of filter initialization, followed by the
chaining of filter cascades:

.. code:: c
   :number-lines: 1

   cascadestate_t stage1;
   cascadestate_t stage2;
   
   void initcascade() {
       createcascade(  /* b0 */  1.0f,
                       /* b1 */  0.0f,
                       /* b2 */  1.0f,
                       /* a1 */ -0.7071f,
                       /* a2 */ 0.25f,
                       &stage1);
       createcascade(  /* b0 */  1.0f,
                       /* b1 */  1.0f,
                       /* b2 */  0.0f,
                       /* a1 */ +0.7071f,
                       /* a2 */ 0.25f,
                       &stage2);
   }
   
   uint16_t processCascade(uint16_t x) {
   
       float32_t input = xlaudio_adc14_to_f32(0x1800 + rand() % 0x1000);
       float32_t v;
   
       v = cascadeiir(input, &stage1);
       v = cascadeiir(v, &stage2);
   
       return xlaudio_f32_to_dac14(v*0.125);
   }

Transposed Structures
^^^^^^^^^^^^^^^^^^^^^

The transposition theorem says that the input-output properties of a network remain unchanged after the following sequence of transformations:

1. Reverse the direction of all branches

2. Change branch points into summing nodes and summing nodes into branch points

3. Interchange the input and output

.. figure:: images/transposetheorem.jpg
   :figwidth: 500px
   :align: center

Once you have the defined the transposed-form structure, you can develop a C program for that implementation.
The transposed-form IIR can be implemented using the same data structure as the direct-form IIR; only the filter operation would be rewritten. The following is the implementation of a second-order
section of a transposed direct-form II design.

.. code:: c
   :number-lines: 1

   typedef struct cascadestate {
       float32_t s[2];   // state
       float32_t b[3];  // nominator coeff  b0 b1 b2
       float32_t a[2];  // denominator coeff   a1 a2
   } cascadestate_t;
   
   float32_t cascadeiir_transpose(float32_t x, cascadestate_t *p) {
       float32_t y = (x * p->b[0]) + p->s[0];
       p->s[0]     = (x * p->b[1]) - (y * p->a[0]) + p->s[1];
       p->s[1]     = (x * p->b[2]) - (y * p->a[1]);
       return y;
   }

Parallel Structures
^^^^^^^^^^^^^^^^^^^

IIR filters are rarely implemented as single, monolithic structures; the risk for
instability is too high. One way to achieve this, is to use a cascade expansion
of an IIR filter.

An alternate implementation
is the parallel-form implementation. The idea of a parallel-from implementation
is to build two or more parallel filters, that can be summed up together to form the
overall transfer function.

.. figure:: images/paralleldesign.jpg
   :figwidth: 600px
   :align: center

Let's consider this for the pole plot shown in the figure. This design has two
poles, so it's transfer function would be:

.. math::

   G(z) =& \frac{1}{(z + 0.25).(z - 0.5)} \\
        =& \frac{z^{-2}}{(1 + 0.25.z^{-1})(1 - 0.5.z^{-1})}

To build a parallel design, we have to decompose :math:`G(z)`

.. math::

   G(z) =& G1(z) + G2(z)

The design of these partial function proceeds by *partial fraction expansion*.
We split the poles of the overall function :math:`G(z)` in two.
In this case, one pole goes with :math:`G1(z)` and the other pole goes with
:math:`G2(z)`. Partial fraction expansion will now proceed by looking for terms
*A* and *B* such that:

.. math::

   \frac{A}{(1 + 0.25.z^{-1})} + \frac{B}{(1 - 0.5.z^{-1})} =  \frac{z^{-2}}{(1 + 0.25.z^{-1})(1 - 0.5.z^{-1})}  

Solving for A and B we find:

.. math::

   G(z) &=  \frac{-1.33333}{(1 + 0.25.z^{-1})} + \frac{1.3333}{(1 - 0.5.z^{-1})}

The meaning of :math:`G(z) = G1(z) + G2(z)` is that a second order system can be implemented by
two independent first-order systems, whose output is added together.

Using Matlab
^^^^^^^^^^^^

In the following, we illustrate a few useful Matlab functions that help with the analysis and
decomposition of IIR filters.

Finding Poles and Zeroes
""""""""""""""""""""""""

.. attention::

   To identify pole and zero locations, use the **roots** function in Matlab,
   which finds the roots of a polynomial. For a polynomial 
   :math:`P(z) = a_0 + a_1 . z^{-1} + ... + a_n . z^{-n}`,
   the roots can be found through ``roots([a0 a1 ... an])``

Example: Find the zeroes and poles of

.. math::

   G(z) = \frac{1 + 0.2 . z^{-1} - z^{-2}}{ 1 - 0.5 . z^{-2}}


The *roots* function in Matlab will return the locations of zeroes and poles:

.. code::

   >> roots([1 0.2 -1])  % numerator coefficients

   ans =

      -1.1050
       0.9050

   >> roots([1 0 -0.5])  % denominator coefficients

   ans =
   
       0.7071
      -0.7071

   >> zplane([1 0.2 -1],[1 0 -0.5])  % create a pole-zero plot


.. figure:: images/zplaneiirexample.png
   :figwidth: 400px
   :align: center

We can write G(z) in factored form as follows

.. math::

   G(z) =& \frac{(1 - n_0 . z^{-1})(1 - n_1 . z^{-1})}{(1 - p_0 . z^{-1})(1 - p_1 . z^{-1}])} \\
   G(z) =& \frac{(1 + 1.1050.z^{-1})(1 - 0.9050.z^{-1})}{(1 + 0.7071.z^{-1})(1 - 0.7071.z^{-1})} 

Finding Partial Fraction Expensions
"""""""""""""""""""""""""""""""""""

.. attention::

   To perform partial fraction expansion, use the **residue** function in Matlab.
   For a transfer function :math:`G(z) = \frac{b_n.z^{-n}}{a_0 + a_1.z^{-1} + ... +  a_n.z^{-n}}`
   the partial fraction is found as ``[r,p,k] = residue([0 0 ... bn], [a0 a1 .. an])``.

The parallel form of an IIR requires partial fraction expansion. Consider the following example,
where the objective is to find *A* and *B*:

.. math::

   G(z) =& \frac{z^{-2}}{(1 + 0.25.z^{-1})(1 - 0.5.z^{-1})} \\
        =& \frac{A}{(1 + 0.25.z^{-1})} + \frac{B}{(1 - 0.5.z^{-1})}

.. code::

   >> % compute G(z)
   >> a = conv([1 0.25],[1 -0.5])

   a =

      1.0000   -0.2500   -0.1250

   % In other words, G(z) = z^-2 / (1 - 0.25.z^-1 - 0.125.z^-2)

   >> % compute terms A and B using 'residue'
   >> [r,p,k] = residue([0 0 1],a)

   r =
   
       1.3333
      -1.3333
   
   
   p =
   
       0.5000
      -0.2500
   
   
   k =
   
        []


From the r term in the Matlab output, we conclude that

.. math::

   G(z) &=  \frac{-1.33333}{(1 + 0.25.z^{-1})} + \frac{1.3333}{(1 - 0.5.z^{-1})}


Conclusions
^^^^^^^^^^^

We discussed the following filter implementation techniques:

  * Direct Form I filters
  * Direct Form II filters
  * Transposed Direct Form II filters

We also discussed two decomposition techniques:

  * Cascade-form decomposition
  * Parallel-form decomposition