• If you are citizen of an European Union member nation, you may not use this service unless you are at least 16 years old.

View

# Bilinear Transform

last edited by 5 years, 8 months ago

update 1/21/2015  (at bottom)

my work from the '90s led to the following

duplicated:

== Bilinear Conversion ==

Provided sampling rate is not far in excess of the largest frequency component,

a quick term-to-term conversion can use as approximation:

z − term = e aT * [1 − A * (4 * a / ws)B]

where Ws=sampling rate <Nyquist rate (rad/s)

T=sampling period, sec =2PI/Ws

a=s-term, real or complex

4w/Ws <=1        (A,B)=specific to system

for example, the s-term  (s+5.0) with T=2PI/40 yields z-term (z-0.4129)

Similarly, term s^2+2.31s+2.72 becomes z^2-1.64z+0.70

there are existing rules for introducing additional z-terms as needed to stabilize and improve magnitude match.

For moderate sampling rate the bilinear transform sees a divergence region for frequency response towards the Nyquist limit.

This can be improved by applying a slightly increased sampling rate up to the original Nyquist limit, at least with well conditioned transfer functions ie positive real.

ref my '88 graduate paper The Role of Positive Real Functions in A/D Transfer Function Matching.I have converted some of that code to C++ now.

Jan 21 2015 update  (in progress, patience)

ok my research paper is now online (see thesis.gmu.edu).  search for mceachen

I will hopefully add a "new" App D expanding on a alternate freq response matching.  Specifically, it demonstrates how to add a "lnz" term to effect.  Essentially, we have to manually compute the change to magnitude and phase from the new term, but in something like PariGP it is very quick.

Rather than the lousy approx used by the bilinear transform, we make use of a cubic fit.

I made an arbitrary transfer function (tf) of order 8.  The additional term nicely corrects the original phase error as can be seen in this file:  \$\$pending.  I chose a sampling rate of 50 rad/s

The basic trick is to realize the substitution exp(jwT) for each z variable  done to lnz merely leaves jWT.  The manual calc is shown here:  \$\$pending.  The PariGP script one can use is here:  \$\$pending

Now, depending on (Q,a) for the added term  (a + (lnz)^Q )  where we expect both  0< a,Q <1  generally.

If one iterates with the Pari script, or uses combinations of such terms, virtually any phase error can be eliminated, provided the error varies as a f{frequency}.  Recall discrete response is taken over 0 thru Pi radians.  Recall numerator term has its magnitude multiplying the overall, while adding phase (lead).  Denominator magnitude divides the overall, and subtracts phase (lags).

For example, using Q=0.95 and a=1, taking freq steps of Pi/10, we have (as a pole):

w   , magn  , phase (deg)

0   , 1  , 0

Pi/10  , 0.9968  , -1.897

Pi/5  , 0.9923  , -3.918

3Pi/10  , 0.9864  , -5.975

2Pi/5  , 0.9792  , -8.041

Pi/2  , 0.9707  , -10.1

3Pi/5  , 0.9611  , -12.14

7Pi/10  , 0.9503  , -14.16

4Pi/5  , 0.9386  , -16.15

9Pi/10  , 0.9261  , -18.1

Pi  , 0.9127  , -20

Let's verify the step response using the substitution lnz ~ 6.3853z^3  -13.476z^2 + 10.481 * z − 3.3064.

We will do the step response manually via the standard manipulation eg    Y/X = z^3+2z^2+4z+5 /  z^3+5z^2+7z+1  becoming

Y(1+5z^-1 + 7z^-2 + z^-3) = X (1 + 2z^-1 + 4z^-2+5z^-3)   where z^-1 is a delay.  This yields :

Y(k)+5Y(k-1)+7Y(k-2)+Y(k-3) = X +2X(k-1)+4X(k-2)+5X(k-3)

So, we pgm our step difference eqn (again into Pari) and can generate the step response.  It can be seen here:  \$\$Pending

By the way, in lieu of 2(z-1)/(z+1), via Tustin, there is 3(z^2-1)/(z^2+4z+1)  which approximates the ln much better.  Details to follow...