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

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

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

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

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)