Forums

characterizing frequency domain performance of a digital controller

Started by jneudorf August 22, 2008
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?  


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
>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');
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&#2013266103;omega. > > Recognize that sigma equals zero and omega equals (2PI&#2013266103;frequency) when > computing frequency response. > > Recognize that exp(j&#2013266103;omega) = cos(omega) + j&#2013266103;sin(omega). > > Plug z = cos(omega) + j&#2013266103;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:
>> 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
>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:
>> 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