Scilab bilinear transform

Started by Andrew Holme May 29, 2007
I'm using Scilab to convert an s-domain filter design to the z-domain for a 
PLL simulation.  It works fine with a time-step of one micro-second 
(Ts=1e-6): the denominator contains z^4 terms, which I was expecting; and 
the PLL simulation agrees with the real hardware; but I only get z^3 terms 
if I reduce the time-step to 1e-7; and the simulation is unstable (it 
oscillates at 3 KHz).  What am I doing wrong?  Am I asking too much of 
Scilab or the bilinear transform?  Setting various pre-warping frequencies 
doesn't seem to make much difference.  Here's my code:

s = poly(0,'s');

// Active feedback network
c15 = 3.3e-6;
c12 = 330e-9;
R7 = 180;

z1 = 1/s/c15 + R7;
z2 = 1/s/c12;

// Input RC (pole)
R4=1.8e3;
R5=1.5e3;
C13=47e-9;

// Output RC (pole)
R6=220;
C14=100e-9;

// Loop filter transfer function
f = 1 / (R4+R5) / (1+s*C13/(1/R4+1/R5)) / (1/z1+1/z2) / (1+s*R6*C14);

clc;
format(20,'e');
f

z=poly(0,'z');
Ts=1e-6;            // Try 1e-7

// One way of doing it
horner(f, (2/Ts)*(z-1)/(z+1))

// A different way
ctf=syslin('c', f);
css=tf2ss(ctf, 1e-20);
dss=cls2dls(css, Ts);
dtf=ss2tf(dss, 1e-20)



On Tue, 29 May 2007 12:00:06 +0100, Andrew Holme wrote:

> I'm using Scilab to convert an s-domain filter design to the z-domain for a > PLL simulation. It works fine with a time-step of one micro-second > (Ts=1e-6): the denominator contains z^4 terms, which I was expecting; and > the PLL simulation agrees with the real hardware; but I only get z^3 terms > if I reduce the time-step to 1e-7; and the simulation is unstable (it > oscillates at 3 KHz). What am I doing wrong? Am I asking too much of > Scilab or the bilinear transform? Setting various pre-warping frequencies > doesn't seem to make much difference. Here's my code: >
-- specifics snipped -- Why, oh why, are you converting this to discrete time just for a simulation? I'm going to rant at you about this, in a separate response, shortly. I played around with it about as much as I want to any morning before the tea kicks in, and didn't learn too much. The bilinear transform, by itself, should always cough up a stable system in the z domain for a stable system in the s domain. It may be the _wrong_ system, but stability properties should be preserved. Polynomial math gets very nasty very quickly as the order goes up, and Scilab is sensitive to this. I have found that many of the Scilab functions will do poorly if you give them transfer functions with poorly conditioned polynomials. On the other hand, it does quite well manipulating systems that are represented in state space, and it supports the '*' and '/' operators for state-space systems, treating them as the appropriate control-system operations. I tried out your code, _except_ that I wrapped all the transfer functions with 'tf2ss' in the calculation of f, and things appeared to work just fine. I could use 'cls2dls' to get the right answer (and I could use 'dscr' to good advantage). I suggest you give that a try. -- Tim Wescott Control systems and communications consulting http://www.wescottdesign.com Need to learn how to apply control theory in your embedded system? "Applied Control Theory for Embedded Systems" by Tim Wescott Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html
On Tue, 29 May 2007 12:00:06 +0100, Andrew Holme wrote:

> I'm using Scilab to convert an s-domain filter design to the z-domain for a > PLL simulation. It works fine with a time-step of one micro-second > (Ts=1e-6): the denominator contains z^4 terms, which I was expecting; and
-- snip -- Whoa, Nelly! <rant> For a PLL simulation? Do you _know_ that your simulation is good enough? Why, when you're handed a sophisticated math package that can do all sorts of cool nonlinear variable-step differential equation solving, are you taking a part of your system and converting it to a fixed-step discrete time using an approximation that's known to be inaccurate?!? I will grant you, the Scicos editor can be a pain to work with. In fact, I've heard rumors that the Scilab team is addressing this (finally -- and thanks). But a PLL is such a supremely nonlinear thing, wouldn't you _want_ to simulate it right? </rant> -- Tim Wescott Control systems and communications consulting http://www.wescottdesign.com Need to learn how to apply control theory in your embedded system? "Applied Control Theory for Embedded Systems" by Tim Wescott Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html
"Tim Wescott" <tim@seemywebsite.com> wrote in message 
news:OdKdnaV7odbgo8HbnZ2dnUVZ_gOdnZ2d@web-ster.com...
> On Tue, 29 May 2007 12:00:06 +0100, Andrew Holme wrote: > >> I'm using Scilab to convert an s-domain filter design to the z-domain for >> a >> PLL simulation. It works fine with a time-step of one micro-second >> (Ts=1e-6): the denominator contains z^4 terms, which I was expecting; and >> the PLL simulation agrees with the real hardware; but I only get z^3 >> terms >> if I reduce the time-step to 1e-7; and the simulation is unstable (it >> oscillates at 3 KHz). What am I doing wrong? Am I asking too much of >> Scilab or the bilinear transform? Setting various pre-warping >> frequencies >> doesn't seem to make much difference. Here's my code: >> > -- specifics snipped -- > > Why, oh why, are you converting this to discrete time just for a > simulation? I'm going to rant at you about this, in a separate response, > shortly. > > I played around with it about as much as I want to any morning before the > tea kicks in, and didn't learn too much. > > The bilinear transform, by itself, should always cough up a stable system > in the z domain for a stable system in the s domain. It may be the > _wrong_ system, but stability properties should be preserved. > > Polynomial math gets very nasty very quickly as the order goes up, and > Scilab is sensitive to this. I have found that many of the Scilab > functions will do poorly if you give them transfer functions with poorly > conditioned polynomials. On the other hand, it does quite well > manipulating systems that are represented in state space, and it supports > the '*' and '/' operators for state-space systems, treating them as the > appropriate control-system operations. > > I tried out your code, _except_ that I wrapped all the transfer functions > with 'tf2ss' in the calculation of f, and things appeared to work just > fine. I could use 'cls2dls' to get the right answer (and I could use > 'dscr' to good advantage). I suggest you give that a try. > > -- > Tim Wescott > Control systems and communications consulting > http://www.wescottdesign.com > > Need to learn how to apply control theory in your embedded system? > "Applied Control Theory for Embedded Systems" by Tim Wescott > Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html
Hi Tim, Thanks for taking the time to look at this. Why am I doing it this way? Probably because I don't know any better. My PLL simulator is actually written in C++ not Scilab. See http://www.holmea.demon.co.uk/Frac2/Simulate.htm for the full story. I'm not sure what you meant by wrapping the transfer functions with tf2ss. I tried this and got the same result as before: // Loop filter transfer function f1 = 1 / (R4+R5); f2 = 1 / (1+s*C13/(1/R4+1/R5)); f3 = 1 / (1/z1+1/z2); f4 = 1 / (1+s*R6*C14); f1 = tf2ss(f1, 1e-40); f2 = tf2ss(f2, 1e-40); f3 = tf2ss(f3, 1e-40); f4 = tf2ss(f4, 1e-40); f = f1 * f2 * f3 * f4; format(20,'e'); Ts=1e-7; dss=cls2dls(f, Ts); dtf=ss2tf(dss, 1e-40) Just to re-iterate: I only have a problem when Ts=1e-7. Everything works fine with Ts=1e-6.
Andrew Holme wrote:
> "Tim Wescott" <tim@seemywebsite.com> wrote in message > news:OdKdnaV7odbgo8HbnZ2dnUVZ_gOdnZ2d@web-ster.com... >> On Tue, 29 May 2007 12:00:06 +0100, Andrew Holme wrote: >> >>> I'm using Scilab to convert an s-domain filter design to the z-domain for >>> a >>> PLL simulation. It works fine with a time-step of one micro-second >>> (Ts=1e-6): the denominator contains z^4 terms, which I was expecting; and >>> the PLL simulation agrees with the real hardware; but I only get z^3 >>> terms >>> if I reduce the time-step to 1e-7; and the simulation is unstable (it >>> oscillates at 3 KHz). What am I doing wrong? Am I asking too much of >>> Scilab or the bilinear transform? Setting various pre-warping >>> frequencies >>> doesn't seem to make much difference. Here's my code: >>> >> -- specifics snipped -- >> >> Why, oh why, are you converting this to discrete time just for a >> simulation? I'm going to rant at you about this, in a separate response, >> shortly. >> >> I played around with it about as much as I want to any morning before the >> tea kicks in, and didn't learn too much. >> >> The bilinear transform, by itself, should always cough up a stable system >> in the z domain for a stable system in the s domain. It may be the >> _wrong_ system, but stability properties should be preserved. >> >> Polynomial math gets very nasty very quickly as the order goes up, and >> Scilab is sensitive to this. I have found that many of the Scilab >> functions will do poorly if you give them transfer functions with poorly >> conditioned polynomials. On the other hand, it does quite well >> manipulating systems that are represented in state space, and it supports >> the '*' and '/' operators for state-space systems, treating them as the >> appropriate control-system operations. >> >> I tried out your code, _except_ that I wrapped all the transfer functions >> with 'tf2ss' in the calculation of f, and things appeared to work just >> fine. I could use 'cls2dls' to get the right answer (and I could use >> 'dscr' to good advantage). I suggest you give that a try. >> >> -- >> Tim Wescott >> Control systems and communications consulting >> http://www.wescottdesign.com >> >> Need to learn how to apply control theory in your embedded system? >> "Applied Control Theory for Embedded Systems" by Tim Wescott >> Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html > > Hi Tim, > > Thanks for taking the time to look at this. Why am I doing it this way? > Probably because I don't know any better. My PLL simulator is actually > written in C++ not Scilab. See > http://www.holmea.demon.co.uk/Frac2/Simulate.htm for the full story. > > I'm not sure what you meant by wrapping the transfer functions with tf2ss. > I tried this and got the same result as before: > > // Loop filter transfer function > f1 = 1 / (R4+R5); > f2 = 1 / (1+s*C13/(1/R4+1/R5)); > f3 = 1 / (1/z1+1/z2); > f4 = 1 / (1+s*R6*C14); > > f1 = tf2ss(f1, 1e-40); > f2 = tf2ss(f2, 1e-40); > f3 = tf2ss(f3, 1e-40); > f4 = tf2ss(f4, 1e-40); > > f = f1 * f2 * f3 * f4; > > format(20,'e'); > Ts=1e-7; > dss=cls2dls(f, Ts); > dtf=ss2tf(dss, 1e-40) > > Just to re-iterate: I only have a problem when Ts=1e-7. Everything works > fine with Ts=1e-6. > > >
When I did this the state space version had all four states, just like you'd want to see. I suspect that you're having numerical problems in the polynomials -- each time you bump up the order of your polynomial, you drastically increase the amount of precision required to accurately locate the roots. If you have four 1st-order transfer functions, should maintain them as independent entities for the simulation and cascade them numerically, or you should do your computations on the state-space version. Either of these will not only retain the full polynomial order, but avoid the exact same numerical difficulties in computing the filter updates as you're running into computing the polynomials themselves. For just a hair more on this issue, see my paper http://www.wescottdesign.com/articles/zTransform/z-transforms.html. Skip down to section 3: "Z Transforms and Software". Alternately, look in any book on DSP methods (like Lyons: "Understanding Digital Signal Processing"). You'll be told there what a bad idea it is to use filters that are one degree higher than you need. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" gives you just what it says. See details at http://www.wescottdesign.com/actfes/actfes.html
"Tim Wescott" <tim@seemywebsite.com> wrote in message 
news:q86dncxInOPH9MHbnZ2dnUVZ_jednZ2d@web-ster.com...
> Andrew Holme wrote: >> "Tim Wescott" <tim@seemywebsite.com> wrote in message >> news:OdKdnaV7odbgo8HbnZ2dnUVZ_gOdnZ2d@web-ster.com... >>> On Tue, 29 May 2007 12:00:06 +0100, Andrew Holme wrote: >>> >>>> I'm using Scilab to convert an s-domain filter design to the z-domain >>>> for a >>>> PLL simulation. It works fine with a time-step of one micro-second >>>> (Ts=1e-6): the denominator contains z^4 terms, which I was expecting; >>>> and >>>> the PLL simulation agrees with the real hardware; but I only get z^3 >>>> terms >>>> if I reduce the time-step to 1e-7; and the simulation is unstable (it >>>> oscillates at 3 KHz). What am I doing wrong? Am I asking too much of >>>> Scilab or the bilinear transform? Setting various pre-warping >>>> frequencies >>>> doesn't seem to make much difference. Here's my code: >>>> >>> -- specifics snipped -- >>> >>> Why, oh why, are you converting this to discrete time just for a >>> simulation? I'm going to rant at you about this, in a separate >>> response, >>> shortly. >>> >>> I played around with it about as much as I want to any morning before >>> the >>> tea kicks in, and didn't learn too much. >>> >>> The bilinear transform, by itself, should always cough up a stable >>> system >>> in the z domain for a stable system in the s domain. It may be the >>> _wrong_ system, but stability properties should be preserved. >>> >>> Polynomial math gets very nasty very quickly as the order goes up, and >>> Scilab is sensitive to this. I have found that many of the Scilab >>> functions will do poorly if you give them transfer functions with poorly >>> conditioned polynomials. On the other hand, it does quite well >>> manipulating systems that are represented in state space, and it >>> supports >>> the '*' and '/' operators for state-space systems, treating them as the >>> appropriate control-system operations. >>> >>> I tried out your code, _except_ that I wrapped all the transfer >>> functions >>> with 'tf2ss' in the calculation of f, and things appeared to work just >>> fine. I could use 'cls2dls' to get the right answer (and I could use >>> 'dscr' to good advantage). I suggest you give that a try. >>> >>> -- >>> Tim Wescott >>> Control systems and communications consulting >>> http://www.wescottdesign.com >>> >>> Need to learn how to apply control theory in your embedded system? >>> "Applied Control Theory for Embedded Systems" by Tim Wescott >>> Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html >> >> Hi Tim, >> >> Thanks for taking the time to look at this. Why am I doing it this way? >> Probably because I don't know any better. My PLL simulator is actually >> written in C++ not Scilab. See >> http://www.holmea.demon.co.uk/Frac2/Simulate.htm for the full story. >> >> I'm not sure what you meant by wrapping the transfer functions with >> tf2ss. I tried this and got the same result as before: >> >> // Loop filter transfer function >> f1 = 1 / (R4+R5); >> f2 = 1 / (1+s*C13/(1/R4+1/R5)); >> f3 = 1 / (1/z1+1/z2); >> f4 = 1 / (1+s*R6*C14); >> >> f1 = tf2ss(f1, 1e-40); >> f2 = tf2ss(f2, 1e-40); >> f3 = tf2ss(f3, 1e-40); >> f4 = tf2ss(f4, 1e-40); >> >> f = f1 * f2 * f3 * f4; >> >> format(20,'e'); >> Ts=1e-7; >> dss=cls2dls(f, Ts); >> dtf=ss2tf(dss, 1e-40) >> >> Just to re-iterate: I only have a problem when Ts=1e-7. Everything works >> fine with Ts=1e-6. >> >> >> > When I did this the state space version had all four states, just like > you'd want to see. > > I suspect that you're having numerical problems in the polynomials -- > each time you bump up the order of your polynomial, you drastically > increase the amount of precision required to accurately locate the roots. > > If you have four 1st-order transfer functions, should maintain them as > independent entities for the simulation and cascade them numerically, or > you should do your computations on the state-space version. Either of > these will not only retain the full polynomial order, but avoid the exact > same numerical difficulties in computing the filter updates as you're > running into computing the polynomials themselves. > > For just a hair more on this issue, see my paper > http://www.wescottdesign.com/articles/zTransform/z-transforms.html. Skip > down to section 3: "Z Transforms and Software". Alternately, look in any > book on DSP methods (like Lyons: "Understanding Digital Signal > Processing"). You'll be told there what a bad idea it is to use filters > that are one degree higher than you need.
Thanks again Tim. I've split the filter into 3 parts and the PLL sim is now working with Ts=1e-7. What's more, the background noise level due to the simulator itself has dropped 40dB. I was hoping it might, due to the finer timestep, but I suppose filter precision could have as much to do with it.