jneudorf wrote:
>> You could also use "bode((25 * %z - 23) / (%z - 1))".  It would be
> easier.
> 
> The reason I did it the way I did was that I need to have a mixed
> system--the rest of the system gain simplifies to Vo/sL, plus the ADC and a
> trailing-edge pulse width modulator.  If I multiply the digital part by
> that, substituting in 1/s = T/log(z), scilab seems to be unhappy, probably
> because the result is no longer a polynomial in z; 
> 
> Side note: I'm a newbie with scilab, so I'm not sure what I need to do
> before I call bode as above: with a fresh session I get an error:
> 
> -->bode((25 * %z - 23)/(%z-1))
>  !--error 10000
> Undefined time domain (sl(4))
> at line      23 of function bode called by :
> bode((25 * %z - 23)/(%z-1))  

Oops -- sorry about that.  When you do (25 * %z - 23) / (%z - 1) you 
make a transfer function in z, but Scilab gives you the ability to 
define the time domain, and the bode plot needs to know (many other 
things you can do with transfer functions in z don't).  So you'd use 
syslin(<your sampling interval>, (25 * %z - 23) / (%z - 1)), etc.

> 
> My intuition based on the searches I've done is that people generally use
> one transform or other (matched z, bilinear etc) to work entirely in one
> domain or other, and get approximate results; my issue with this is that I
> don't really understand the frequency and phase warping effects of this and
> I want to see actual frequency domain results.
> 
> Perhaps if I had more experience, I could just use rule-of-thumb
> approximations that everyone in the field knows are "good enough".

Unfortunately, the technique you're using also is an approximation, in 
that it ignores the filtering effect of your output process.

If you can bear doing all your work in the sampled-time domain you can 
do this with an exact conversion of the plant model.  Furthermore, if 
your ADC looks like a classical sampling process (i.e. it samples at the 
instant that it gets a conversion command) and if your output process 
looks like a zero order hold (i.e. a traditional DAC, where you write a 
value and that's where the voltage goes) then Scilab provides you with 
the conversion function, 'dscr':

Hplant = syslin('c', Vo / L * (1 / %s));

Hdp = dscr(<your sampling interval>, Hplant);

'Hdp' here is an _exact_ model of the plant from your controller's point 
of view.  This is of great advantage if you can do all of your stability 
and performance analysis in the sampled-time domain, but it does become 
an approximation if you try to go back to continuous-time.

Then you can just do arithmetic with your plant transfer function and 
your controller transfer function.

> Any further suggestions on how far off I am would be greatly appreciated,
> and I'm really thankful for the help so far.
> 
> //This is how far I've gotten so far, still need to add in ADC and PWM
> f = (1e+6^0.0001)^(0:9999);
> w  = f*2*%pi;
> T  = 1/90000;
> s  = %i*w;
> jwT = %i*w*T;
> z = %e^jwT;
> g  = (25*z-23)/32 ./(z-1) ./s *31000;
> subplot(211);
> plot2d(f,abs(g),logflag='ll');
> plot2d(f,0*f+1,logflag='ll'); //0db Gain
> subplot(212);
> plot2d(f,180/%pi*atan( imag(g) , real(g) ),logflag='ln');
> plot2d(f,0*f,logflag='ln'); //0 phase
> 

Play with these:

Hdp = dscr(T, syslin('c', Vo / L / %s));      // plant transfer
Hc  = syslin(T, (25 * %z - 23) / (%z - 1));   // controller transfer

Hol = Hdp * Hc;    // open-loop transfer

bode(Hol);         // open-loop bode plot
bode(Hol /. 1);    // closed-loop bode for unity gain

Note that if you want to build complex plant models, or complex 
controllers, you should convert your system models to state space, e.g.

Hdp = tf2ss(dscr(...));  // instead of just dscr(...) as above

This avoids a bunch of numerically challenging polynomial math that 
Scilab (or any other numerical tool) can't do well, instead using 
state-space representations that are much more numerically robust.

I taught a course in control systems for software engineers a while 
back, for which I developed some homework problems similar to the above 
to be done in Scilab.  I make ABSOLUTELY NO PROMISES, but you may find 
some useful stuff in the worked-homework examples here: 
http://www.wescottdesign.com/oit/cst407/syllabus.html, particularly in 
the homework that pertains to chapters 5 and 6.

Basically, I did not have time to include much Scilab-specific 
information in my book, so I included mini-tutorials in the homework. 
You may be able to use these as examples for what you want to do.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

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

>You could also use "bode((25 * %z - 23) / (%z - 1))".  It would be
easier.

The reason I did it the way I did was that I need to have a mixed
system--the rest of the system gain simplifies to Vo/sL, plus the ADC and a
trailing-edge pulse width modulator.  If I multiply the digital part by
that, substituting in 1/s = T/log(z), scilab seems to be unhappy, probably
because the result is no longer a polynomial in z; 

Side note: I'm a newbie with scilab, so I'm not sure what I need to do
before I call bode as above: with a fresh session I get an error:

-->bode((25 * %z - 23)/(%z-1))
 !--error 10000
Undefined time domain (sl(4))
at line      23 of function bode called by :
bode((25 * %z - 23)/(%z-1))       

My intuition based on the searches I've done is that people generally use
one transform or other (matched z, bilinear etc) to work entirely in one
domain or other, and get approximate results; my issue with this is that I
don't really understand the frequency and phase warping effects of this and
I want to see actual frequency domain results.

Perhaps if I had more experience, I could just use rule-of-thumb
approximations that everyone in the field knows are "good enough".

Any further suggestions on how far off I am would be greatly appreciated,
and I'm really thankful for the help so far.

//This is how far I've gotten so far, still need to add in ADC and PWM
f = (1e+6^0.0001)^(0:9999);
w  = f*2*%pi;
T  = 1/90000;
s  = %i*w;
jwT = %i*w*T;
z = %e^jwT;
g  = (25*z-23)/32 ./(z-1) ./s *31000;
subplot(211);
plot2d(f,abs(g),logflag='ll');
plot2d(f,0*f+1,logflag='ll'); //0db Gain
subplot(212);
plot2d(f,180/%pi*atan( imag(g) , real(g) ),logflag='ln');
plot2d(f,0*f,logflag='ln'); //0 phase

jneudorf wrote:
>> jneudorf wrote:
>>
>>> I have a digital controller that has been hand-tuned to sort-of-work.
>>> [ G[z] =3D ( 25z - 23 ) / (z-1) and sampling frequency is 90kHz ].  I
> wou=
>> ld
>>> like to know how to determine its gain at various frequencies.
>> It's very simple from first-principles:
>>
> Thanks Greg.
> Lots of implications to think about.  For example, the polynomial
> expansion of sin and cos are infinite, and exp(sT)-1 has an infinite number
> of roots.   So it becomes clear that even a simple digital filter has
> properties totally different from analog filters.  The response at the
> multiples of the sampling frequency seems to be the same as the result at
> DC--in this case, infinite gain.
> 
> I've written a scilab script to show me the results (hopefully the web
> interface I'm using won't screw it up too much--the 3Ds and B7s that got
> inserted last time were annoying ):
> 
> f = (1e+6^0.0001)^(0:9999); //10K values from 1 to 1 million
> w  = f*2*%pi;
> T  = 1/90000; //sampling rate
> jwT = %i*w*T;
> z = %e^jwT;
> g  = (25*z-23) ./(z-1);
> plot2d(f,abs(g),logflag='ll');
> 
> 
You could also use "bode((25 * %z - 23) / (%z - 1))".  It would be easier.

You are correct that the gain of the transfer function goes to infinity 
at the sampling rate, but in a real-world situation this doesn't mean 
that the system will oscillate at the sampling rate -- it means that any 
noise at the sampling rate will alias down to DC in your sampler, and 
create a DC offset in your measurement which will cause a corresponding 
offset in the point your system servos to.

Since the usual output of a digital controller is a DAC that acts like a 
zero-order hold, the controller doesn't inject any energy at the sample 
rate back into the system.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

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

On Aug 23, 5:02 am, Greg Berchin <gberc...@sentientscience.com> wrote:
> jneudorf wrote:
> > I have a digital controller that has been hand-tuned to sort-of-work.
> > [ G[z] = ( 25z - 23 ) / (z-1) and sampling frequency is 90kHz ].  I would
> > like to know how to determine its gain at various frequencies.
>
> It's very simple from first-principles:
>
> Recognize that z = exp(sT).
>
> Recognize that T = 1/(90kHz).
>
> Recognize that s = sigma + j&#4294967295;omega.
>
> Recognize that sigma equals zero and omega equals (2PI&#4294967295;frequency) when
> computing frequency response.
>
> Recognize that exp(j&#4294967295;omega) = cos(omega) + j&#4294967295;sin(omega).
>
> Plug z = cos(omega) + j&#4294967295;sin(omega) into G[z] and solve for magnitude
> and phase at whatever frequencies you desire.
>
> > Or should I have posted this in one of the sci.eng groups?
>
> Comp.dsp is perfectly appropriate.
>
> Greg

Substitute z=exp(j theta)

where theta=omega.T where T is sampling interval omega is angular freq
rad/s
Then z=cos(theta)+jsin(theta) etc
exclude omega=0 since the dc gain of an integrator is infinity.

remember that 2 integrators are always better than 1 so you could do
better with another P-I in cascade possibly. ie(1+sT2)/S


K.

>jneudorf wrote:
>
>> I have a digital controller that has been hand-tuned to sort-of-work.
>> [ G[z] =3D ( 25z - 23 ) / (z-1) and sampling frequency is 90kHz ].  I
wou=
>ld
>> like to know how to determine its gain at various frequencies.
>
>It's very simple from first-principles:
>
Thanks Greg.
Lots of implications to think about.  For example, the polynomial
expansion of sin and cos are infinite, and exp(sT)-1 has an infinite number
of roots.   So it becomes clear that even a simple digital filter has
properties totally different from analog filters.  The response at the
multiples of the sampling frequency seems to be the same as the result at
DC--in this case, infinite gain.

I've written a scilab script to show me the results (hopefully the web
interface I'm using won't screw it up too much--the 3Ds and B7s that got
inserted last time were annoying ):

f = (1e+6^0.0001)^(0:9999); //10K values from 1 to 1 million
w  = f*2*%pi;
T  = 1/90000; //sampling rate
jwT = %i*w*T;
z = %e^jwT;
g  = (25*z-23) ./(z-1);
plot2d(f,abs(g),logflag='ll');

jneudorf wrote:

> I have a digital controller that has been hand-tuned to sort-of-work.
> [ G[z] = ( 25z - 23 ) / (z-1) and sampling frequency is 90kHz ].  I would
> like to know how to determine its gain at various frequencies.

It's very simple from first-principles:

Recognize that z = exp(sT).

Recognize that T = 1/(90kHz).

Recognize that s = sigma + j�omega.

Recognize that sigma equals zero and omega equals (2PI�frequency) when
computing frequency response.

Recognize that exp(j�omega) = cos(omega) + j�sin(omega).

Plug z = cos(omega) + j�sin(omega) into G[z] and solve for magnitude
and phase at whatever frequencies you desire.

> Or should I have posted this in one of the sci.eng groups?

Comp.dsp is perfectly appropriate.

Greg

I have a digital controller that has been hand-tuned to sort-of-work. 
[ G[z] = ( 25z - 23 ) / (z-1) and sampling frequency is 90kHz ].  I would
like to know how to determine its gain at various frequencies.  If it were
an ordinary filter, I could take the fft of its impulse response.  As it
is, I've tried stimulating it with random noise (vector x), to get a
response, (vector y), then done an fft on x and y, and computed the gain as
fft(y)/fft(x).  But this is very noisy.  I could stimulate it with sine
waves of the specific frequencies I'm interested in, but this seems messy.

How would other people approach this problem?  It seems most people work
the other way, from s-domain to z-domain, but this opens up a can of worms
with respect to trapezoidal vs rectangular integration, warping, etc, and
if I work this way, I still won't be able to predict how far my digital
implementation deviates from the s domain solution.

Or should I have posted this in one of the sci.eng groups?