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

  • You already know Dokkio is an AI-powered assistant to organize & manage your digital files & messages. Very soon, Dokkio will support Outlook as well as One Drive. Check it out today!

View
 

Bilinear Transform

Page history last edited by Bill McEachen 1 year, 4 months ago

update 1/21/2015  (at bottom)

 

my work from the '90s led to the following

 

orig post:  http://en.wikibooks.org/wiki/Talk:Digital_Signal_Processing/Bilinear_Transform

 

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 , try  ( 0.85, 3.17 )

 

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

Similarly, term s^2+2.31s+2.72 (roots s= -1.155 +/-  1.177i  )       becomes z^2 - 1.64 z + 0.70   (roots z= 0.82 +/- 0.1661*I )

using    exp((-1.155+1.177*I)*T)*(1- 0.85*(4*(1.155+1.177*I)/w)^3.17)   = 0.8220 + 0.1523*I

 

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...

 

****************************

Nov 2022 update  (online links, Octave stuff )

 

https://octave-online.net/

 

https://octave.sourceforge.io/    (old but useful )

 

sys(b,a,tsam)

 

https://octave.sourceforge.io/control/function/tf.html    

 

https://octave.sourceforge.io/octave/function/freqz_plot.html    :   freqz_plot (w, h, freq_norm)

https://octave.sourceforge.io/octave/function/freqz.html

 

https://octave.sourceforge.io/control/function/step.html    step response examples

 

https://octave.sourceforge.io/control/function/@lti/c2d.html    analog-to-digital ie c2d

 

https://octave.discourse.group

 

https://stackoverflow.com/questions/67671222/how-to-compute-an-iir-filter-frequency-response-in-octave-matlab

 

ex.

a = 0.7; % example value
num = 1-a;
den = [1, -a];
freqz(num, den, 1001)

 

Sampling time: 0.2 s
Discrete-time model.
octave:3> z = tf ('z', 0.2)
octave:4> H = 0.095/(z-0.9)
octave:5> freqz_plot (Pi,H,freq_norm)
b= [x,x,x] 
a= [y,y,y ] 
freqz(b,a,10)  where b is numer vector

Transfer function 'H' from input 'u1' to output ... 0.095 y1: ------- z - 0.9 **

 

s-plane:

b = [1 2]; a = [1 1];
w = linspace (0, 4, 128);
freqs (b, a, w);

Comments (0)

You don't have permission to comment on this page.