Imp.inv/bilinear trans. vs. expm()

Hi,

What is the reason to that dicrete IIR filters are usually
designed using impulse invariant transformation (= Euler
integration) or bilinear transformation (= trapetzoidal
integration) instead of the closed form solution that
can be computed by the matrix exponential function?

At least in DSP books IIR filter design seems to be taught
like this:

1. Design a continuous-time (analog) filter, for example
   Butterworth, by finding the order and parameters for
   system function H(s) based on the specifications. The
   filter is then

     Y(s) = H(s) U(s)

   where U(s) and Y(s) are the Laplace transformations of
   input and output.

2. Convert the continuous-time filter from s-plane to
   z-plane by using impulse invariant or bilinear transform.
   This gives the discrete IIR system function H(z^-1).

3. The resulting filter is then implemented in time domain
   as a difference equation derived from H(z^-1).

But couldn't this be done also as follows:

1. Design a continuous-time filter, say Butterworth, and convert
   it into state space form

     dx/dt = F x + L u
         y = C x

   where u is the input and y is the output.

2. Construct the equivalent discrete model with given
   sampling period dt (here x[k] is a vector):

     x[k] = A x[k-1] + B u[k-1]
     y[k] = C x[k]

   where transition matrix A = expm(F dt) and B can be
   also computed in closed form. At least A is very easy
   to compute, because its eigenvalues are simple functions
   of the poles of H(s).

This latter solution is theoretically exact, but the impulse
invariant and bilinear transformations are only numerical
integrations. It would seem better idea to use the exact
solution, because it is available. But why do DSP books usually
talk about the imp.inv. and bilinear transformations only and
don't mention this exact solution? Is the computational
complexity the reason to this, or are there other issues? The
exact solution is computationally heavier but a significant
reduction in computations can be achieved by diagonalization
of the matrix A, I think.

-- 
- Simo

Simo S&#4294967295;rkk&#4294967295; <simo.sarkka@hut.fi> writes:

> Hi,
> 
> What is the reason to that dicrete IIR filters are usually
> designed using impulse invariant transformation (= Euler
> integration) or bilinear transformation (= trapetzoidal
> integration) instead of the closed form solution that
> can be computed by the matrix exponential function?
> 
> At least in DSP books IIR filter design seems to be taught
> like this:
> 
> 1. Design a continuous-time (analog) filter, for example
>    Butterworth, by finding the order and parameters for
>    system function H(s) based on the specifications. The
>    filter is then
> 
>      Y(s) = H(s) U(s)
> 
>    where U(s) and Y(s) are the Laplace transformations of
>    input and output.
> 
> 2. Convert the continuous-time filter from s-plane to
>    z-plane by using impulse invariant or bilinear transform.
>    This gives the discrete IIR system function H(z^-1).
> 
> 3. The resulting filter is then implemented in time domain
>    as a difference equation derived from H(z^-1).
> 
> But couldn't this be done also as follows:
> 
> 1. Design a continuous-time filter, say Butterworth, and convert
>    it into state space form
> 
>      dx/dt = F x + L u
>          y = C x
> 
>    where u is the input and y is the output.

Hi Simo,

Your ideas sound very intriguing. Can you please explain how
you would do step 1 above? Let's say we have a simple
1st-order RC lowpass with H(s) = 1/(1 + sRC). How would
you get the F and L and C from this?
-- 
Randy Yates
Sony Ericsson Mobile Communications
Research Triangle Park, NC, USA
randy.yates@sonyericsson.com, 919-472-1124

On Fri, 27 Aug 2004, Randy Yates wrote:

> Simo S&#4294967295;rkk&#4294967295; <simo.sarkka@hut.fi> writes:
>
> > 1. Design a continuous-time filter, say Butterworth, and convert
> >    it into state space form
> >
> >      dx/dt = F x + L u
> >          y = C x
> >
> >    where u is the input and y is the output.
>
> Hi Simo,
>
> Your ideas sound very intriguing. Can you please explain how
> you would do step 1 above? Let's say we have a simple
> 1st-order RC lowpass with H(s) = 1/(1 + sRC). How would
> you get the F and L and C from this?

I guess that the correponding state space model in
this case is simply (by inverse Laplace transform)

  dx/dt = (-1/(RC)) x + (1/(RC)) u
      y = x,

where x is scalar and thus the parameters are

  F = -1/(RC)
  L = 1/(RC)
  C = 1.

-- 
- Simo

Simo S&#4294967295;rkk&#4294967295; <simo.sarkka@hut.fi> wrote in message news:<Pine.OSF.4.58.0408271231210.280764@zeus.lce.hut.fi>...
> Hi,
> 
> What is the reason to that dicrete IIR filters are usually
> designed using impulse invariant transformation (= Euler
> integration) or bilinear transformation (= trapetzoidal
> integration) instead of the closed form solution that
> can be computed by the matrix exponential function?

The impulse invariant transform does not correspond to Euler
integration.

> But couldn't this be done also as follows:
> 
> 1. Design a continuous-time filter, say Butterworth, and convert
>    it into state space form
> 
>      dx/dt = F x + L u
>          y = C x
> 
>    where u is the input and y is the output.
> 
> 2. Construct the equivalent discrete model with given
>    sampling period dt (here x[k] is a vector):
> 
>      x[k] = A x[k-1] + B u[k-1]
>      y[k] = C x[k]
> 
>    where transition matrix A = expm(F dt) and B can be
>    also computed in closed form. At least A is very easy
>    to compute, because its eigenvalues are simple functions
>    of the poles of H(s).
> 
> This latter solution is theoretically exact, but the impulse
> invariant and bilinear transformations are only numerical
> integrations. It would seem better idea to use the exact

This _is_ the impulse invariant transform. Furthermore, it is not
really closed form because you have to compute the poles of H(s). The
ordinary way of computing the impulse invariant transform is just as
much a closed form once you have the poles of H(s).

-Frederick Umminger

On Fri, 27 Aug 2004, Frederick Umminger wrote:

> > 2. Construct the equivalent discrete model with given
> >    sampling period dt (here x[k] is a vector):
> >
> >      x[k] = A x[k-1] + B u[k-1]
> >      y[k] = C x[k]
> >
> >    where transition matrix A = expm(F dt) and B can be
> >    also computed in closed form. At least A is very easy
> >    to compute, because its eigenvalues are simple functions
> >    of the poles of H(s).
> >
> > This latter solution is theoretically exact, but the impulse
> > invariant and bilinear transformations are only numerical
> > integrations. It would seem better idea to use the exact
>
> This _is_ the impulse invariant transform. Furthermore, it is not
> really closed form because you have to compute the poles of H(s). The
> ordinary way of computing the impulse invariant transform is just as
> much a closed form once you have the poles of H(s).

Doesn't this mean that we can compute matrix exponentials by using
impulse invariant transformations? What is the relationship?

It is true that computing the matrix exponential can be done easily if the
poles are known. But that is not the only way of computing the matrix
exponential, and generally the poles are not required to be known. Or what
did you mean by that the poles have to be known?

-- 
- Simo

Simo S&#4294967295;rkk&#4294967295; <simo.sarkka@hut.fi> wrote in message news:<Pine.OSF.4.58.0408271231210.280764@zeus.lce.hut.fi>...
> Hi,
> 
> What is the reason to that dicrete IIR filters are usually
> designed using impulse invariant transformation (= Euler
> integration) or bilinear transformation (= trapetzoidal
> integration) instead of the closed form solution that
> can be computed by the matrix exponential function?

You need to keep an eye on the history of DSP to see the full reasons. 
When DSP was in its infancy (form an operational, not necessarily 
theoretical, point of view), the engineers who designed and implemented
the filters were trained as "analog" engineers, and had little if any 
formal or theoretical training in designing discrete-time systems. 
These engineers were trained to think in terms of Butterorth, Chebychev, 
etc., filters, and knew how to design the filters from an *analog* spec. 

The efficient way of having people with this background and training 
design digital filters, was to let them design an analog counterpart 
and then convert it, e.g. by the bilinear transform, to discrete time.
In this setting, one needs something that works, that is understandable, 
and that does not take forever to compute with pen, paper and a slide 
rule as tools. Theoretical or numerical exactness is not an issue. 
Getting the job done, is. 

Once DSP "matured" in the sense that DSP was taught in courses in 
universities and engineering schools, IIR design by analog template 
was kept as a link to the world of continuous time, while purely 
discrete-time design methods were developed. Unfortunately, most 
of the latter methods I know of are used for FIR filters, not IIR 
filters. Perhaps some of the "old foxes" that hang out here, would 
be kind enough to fill us in on the matter...

Anyway, do note that this method of "analog template" does restrict 
the digital filter in certain ways. The analog filters need to obey 
the rules of physics. Digital filters only have to obey the rules of 
maths, which is not at all the same as physics. 

As for your actual question, matrix exponential functions might look 
good on paper but are usually quite cumbersome in real life. The naive 
matrix exponential would be implemented as something like

                    inf
  E = expm(A) = I + sum  A^n/n!
                    n=1

which is both cumbersome, from a computational point of view, and also 
questionable from a numerical point of view. Another, more sophisticated
implementation would require A to be symmetrical (are the As you think 
of symmetrical?) but requires an eigenvalue decomposition. Eigenvalues
are not always easy to compute.  

> But why do DSP books usually
> talk about the imp.inv. and bilinear transformations only and
> don't mention this exact solution?

There are quite a few books on DSP that do not emphasize the state-space 
model. I am not all that sure you are right in that the exponential 
matrix is exact, and not the bilinear transform. Even if you were, I 
think most DSP engineers would find that using the state-space method 
to design IIR filters from analog templates would amount to use a 
Tomahwak missile to get rid of a rodent in the back garden.

> Is the computational
> complexity the reason to this, or are there other issues? The
> exact solution is computationally heavier but a significant
> reduction in computations can be achieved by diagonalization
> of the matrix A, I think.

Theoretical complexity and comutational drudgery are the two main 
reason I would guess at. Insignificant return for the tons of extra work, 
is the third. It might very well be that some computational savings can 
be obtained for very special A matrixes. I suspect, however, that these 
savings can only be found after intense scrutiny of the particular state 
space model in question. 

Most people would not want to first spend days or weeks on scrutinising 
a linear algebra model when they can use a well-known cook-book receipt
in a matter of hours, or even minutes with the appropriate software.

If one decides to spend a lot of time and computational resources on 
designing an IIR filter, it makes no sense to me that one would work 
inside the *physical* limits of analog template filters. I would expect 
that one threw all of that away and worked in a setting based completely 
on the Z transform and the *mathematical* discrete-time "universe".

Rune

On Mon, 29 Aug 2004, Rune Allnor wrote:

> Anyway, do note that this method of "analog template" does restrict
> the digital filter in certain ways. The analog filters need to obey
> the rules of physics. Digital filters only have to obey the rules of
> maths, which is not at all the same as physics.

In tracking/navigation applications which I have been working on (or
seen), the filters are used for modeling the physical world and the models
are most naturally formulated as (stochastic) continuous-time models.
This formulation also solves the problem of non-uniform sampling (=
asynchronous sensors). The discretization is often done by computing the
matrix exponential in one way or another, which results in model that is
suitable for discrete Kalman filters. This is the background why I asked
about the connection to z-transform approaches.

> As for your actual question, matrix exponential functions might look
> good on paper but are usually quite cumbersome in real life. The naive
> matrix exponential would be implemented as something like
>
>                     inf
>   E = expm(A) = I + sum  A^n/n!
>                     n=1
>
> which is both cumbersome, from a computational point of view, and also
> questionable from a numerical point of view. Another, more sophisticated
> implementation would require A to be symmetrical (are the As you think
> of symmetrical?) but requires an eigenvalue decomposition. Eigenvalues
> are not always easy to compute.

Matlab has three demo methods for computing the matrix exponential.
They are Pade approximation, Taylor series and eigenvalue decomposition.
The first two don't require computing the eigenvalues. I guess that
one feasible method would be to use, for example, Runge-Kutta ingeration
of the differential equation model.

In the most cases A is not symmetrical - in fact I have never seen
a model which would have symmetrical A. But in some cases ariring
in practise, A can be computed in closed form by the Taylor series
expansion, for example if

  A = [0 1 0 0;
       0 0 0 0
       0 0 0 1
       0 0 0 0]

the matrix exponential is easy to compute, because A^2 = 0.
This A arises in 2d-almost-constant-velocity model which is
often used in tracking applications.

-- 
- Simo

Simo S&#4294967295;rkk&#4294967295; <simo.sarkka@hut.fi> wrote in message news:<Pine.OSF.4.58.0408301127360.152496@zeus.lce.hut.fi>...
> On Mon, 29 Aug 2004, Rune Allnor wrote:
> 
> > Anyway, do note that this method of "analog template" does restrict
> > the digital filter in certain ways. The analog filters need to obey
> > the rules of physics. Digital filters only have to obey the rules of
> > maths, which is not at all the same as physics.
> 
> In tracking/navigation applications which I have been working on (or
> seen), the filters are used for modeling the physical world and the models
> are most naturally formulated as (stochastic) continuous-time models.
> This formulation also solves the problem of non-uniform sampling (=
> asynchronous sensors). The discretization is often done by computing the
> matrix exponential in one way or another, which results in model that is
> suitable for discrete Kalman filters. This is the background why I asked
> about the connection to z-transform approaches.

Well, in tracking/navigation aooplications, the link to the physical 
world is obvious, since the output of a filtering algorithm attempts 
to predict where a vessel is, or will be in the near future. Since the 
vessel has to obey the laws of physics, a prediction algorithm must 
necessarily take the laws of physics into account. 

Not so with the digital filter. There is no reason to take for granted 
that the data we try to filter has some physical orign. We could, for 
instance, be talking about financial data, population statistics, or 
the number of books in a bookshelf as function of tage of the book owner.

A digital filter basically operates on sequences of numbers. The filter 
only have to make sense from a mathematical point of view, not a physical 
point of view. Some times one chooses to design a *mathematical* filter 
from a *physical* template. The bilinear transform and the method of 
impulse invariance are the tools to do that. But there are other methods
to design filters, that are not based in physics, only in maths, that 
work equally well. 

I am not to say that you are wrong in that matrix exponentials can be 
used in the way you suggest. I am only saying that they require lots 
of work to understand and implement, while being unecessarily limited 
by physical assumptions and without offering sufficient advantages 
to make up for the extra work.

Rune