home *** CD-ROM | disk | FTP | other *** search
- ## Copyright (C) 1996,1998 Auburn University. All Rights Reserved
- ##
- ## This file is part of Octave.
- ##
- ## Octave is free software; you can redistribute it and/or modify it
- ## under the terms of the GNU General Public License as published by the
- ## Free Software Foundation; either version 2, or (at your option) any
- ## later version.
- ##
- ## Octave is distributed in the hope that it will be useful, but WITHOUT
- ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- ## for more details.
- ##
- ## You should have received a copy of the GNU General Public License
- ## along with Octave; see the file COPYING. If not, write to the Free
- ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
-
- ## -*- texinfo -*-
- ## @deftypefn {Function File } { } frdemo ( )
- ## Octave Controls toolbox demo: Frequency Response demo
- ## @end deftypefn
-
- function frdemo ()
-
- ## Written by David Clem August 15, 1994
-
- ## a s hodel: updated to match new order of ss2zp outputs
- ## J Ingram: updated for system data structure format August 1996
-
- disp("")
- clc
- j = 0;
- while (j != 4)
- disp("");
- j = menu("Octave Controls Systems Toolbox Frequency Response Demo",...
- 'Bode analysis (bode)',...
- 'Nyquist analysis (nyquist)',...
- "Nichols analysis (nichols)", ...
- 'Return to main demo menu');
-
- if (j == 1)
- k1 = 0;
- while (k1 != 4)
- disp("\n");
- clc
-
- k1 = menu("Bode analysis (bode)",...
- 'Continuous system bode analysis',...
- 'Discrete system bode analysis',...
- 'Bode command description', ...
- 'Return to frdemo menu');
-
- if( k1 == 1 )
- disp(" ")
- clc
- disp("\nContinuous system bode analysis\n");
- disp("Example #1:")
- disp("\nConsider the system sys1=");
- sys1=tf2sys([1, 1], [1, 0, -1]);
- sysout(sys1);
- disp("\nPole-zero form can be obtained as follows:")
- cmd = "sysout(sys1,""zp"");";
- run_cmd;
- disp("The systems bode plot is obtained as follows:");
- cmd = "bode(sys1);";
- run_cmd;
- disp("\nNotice that bode automatically labels the plots according to")
- disp("the selected input/output combinations.")
- disp(" ")
- disp("If the frequency range is not specified, bode automatically")
- disp("selects a frequency range based on the natural frequencies of")
- disp("of all poles away from s=0 (or z=1 in discrete time). Bode")
- disp("then checks to make sure that the phase plot is sufficiently")
- disp("smooth that relevant plot behavior is captured.")
- disp("")
- disp("Bode exits with an error if the system is mixed (both continuous")
- disp("and discrete; see is_digital for conditions)")
- prompt
- disp("\nIf the plot magnitude, phase and frequency data is desired, the");
- disp("user can enter the following command:");
- disp("\n[Mag,Phase,w] = bode(sys);");
- disp("\nThis will return three vectors containing the magnitude,");
- disp("phase and frequency.\n");
- prompt;
-
- disp("")
- clc
- disp("Example #2, sys2=")
- cmd = "sys2=zp2sys(1, [-1, -5], 10);";
- eval(cmd);
- cmd = "sysout(sys2);";
- eval(cmd);
- disp("\nThe bode plot command is identical to the tf form:")
- cmd = "bode(sys2);";
- run_cmd;
- disp("\nThe internal representation of the system is not important;")
- disp("bode automatically sorts it out internally.")
- prompt;
-
- disp("")
- clc
- disp("Example #3, Consider the following state space system sys3=:\n");
- cmd = "sys3=ss2sys([0, 1; -1000, -1001], [0; 1], [0, -891], 1);";
- eval(cmd);
- cmd = "sysout(sys3);";
- eval(cmd);
- disp("\nOnce again, the bode plot command is the same:");
- cmd = "bode(sys3);";
- run_cmd;
- disp("\nSuppose the user is interested in the response of the system");
- disp("defined over the input frequency range of 1 - 1000 rad/s.\n");
- disp("First, a frequency vector is required. It can be created");
- disp("with the command:\n");
- cmd = "wrange = logspace(log10(1),log10(1000),100);";
- disp(cmd);
- eval(cmd);
- disp("\nThis creates a logarithmically scaled frequency vector with");
- disp("100 values between 1 and 1000 rad/s\n");
- disp("Then, the bode command includes wrange in the input arguments");
- disp("like this:");
- cmd = "bode(sys3,wrange);";
- run_cmd;
- prompt;
-
- disp("")
- clc
- disp("\nExample #4, The state-space system from example 3 will be");
- disp("grouped with the system from example 2 to form a MIMO system");
- disp("The commands to do this grouping are as follows (changing signal");
- disp("names for clarity):");
- cmd = "sys2 = syssetsignals(sys2,\"out\",\"y_sys2\");";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"in\",\"u_sys2\");";
- disp(cmd); eval(cmd);
- cmd = "nn = sysdimensions(sys2);";
- disp(cmd); eval(cmd);
- cmd = "[nn,nz] = sysdimensions(sys2);";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"st\",sysdefioname(nn+nz,\"x_sys2\"));";
- disp(cmd); eval(cmd);
- cmd = "sys_mimo = sysgroup(sys2,sys3);";
- disp(cmd); eval(cmd);
- disp("The resulting state-space system (after changing signal names");
- disp("in sys2) is");
- cmd = "sysout(sys_mimo)";
- eval(cmd);
- disp("\nNotice that there are now 2 inputs and 2 outputs, and that it did");
- disp("not matter what form the two systems were in when they were grouped.");
- disp(["\nTo view the system\'s bode plots, execute the",...
- " following command:\n"])
- cmd = "bode(sys_mimo);";
- run_cmd;
- prompt
- disp("\nTo view the bode plots for selected channels, the command form changes:")
- cmd = "wrange = [];";
- disp(cmd)
- eval(cmd);
- cmd = "out = 1;";
- disp(cmd)
- eval(cmd);
- cmd = "in = 1;";
- disp(cmd)
- eval(cmd);
- cmd = "bode(sys_mimo,wrange,out,in);";
- run_cmd;
- disp("\nNotice that this bode plot is the same as the plot from example 2.");
- prompt
- closeplot
-
- elseif( k1 == 2 )
- disp("")
- clc
- disp("\nDiscrete system bode analysis\n");
- disp("Display bode plots of a discrete SISO system (dbode)\n")
- disp("Example #1, Consider the following discrete transfer");
- disp(" function:\n");
- cmd = "sys1 = tf2sys([0.00100502, -0.00099502], [1, -2, 1], 0.001);";
- disp(cmd);
- eval(cmd);
- cmd = "sysout(sys1)";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine open loop zeros and poles of the system,");
- disp("use the command:\n")
- cmd = "sysout(sys1,""zp"");";
- run_cmd;
- disp("\nTo view the system's bode plots, execute the following");
- disp("command:\n")
- cmd = "bode(sys1);";
- run_cmd;
- disp("\nNotice (1) the plot label uses exp(jwT) for its title axis. This")
- disp(" allows the user to determine what kind of system was")
- disp(" used to generate the bode plot");
- disp(" (2) the system poles are both at z=1, (break frequency at")
- disp(" jwT = 0); pure integrator poles like this are discarded")
- disp(" by Octave when computing the plot frequency range.")
-
- disp("\nIf magnitude, phase, and frequency data are also desired,");
- disp(" perform the following command instead:\n");
- disp("[M,P,w]=dbode(num,den,T,wrange).\n Where:");
- disp("M => Bode magnitude response data");
- disp("P => Bode phase response data");
- disp("w => frequencies that M and P were evaluated at");
- disp("sys1 => system data structure")
- disp("T => sample period")
- disp("wrange => optional vector of frequencies")
- disp(" if wrange is entered in the argument list, the");
- disp(" system will be evaluated at these specific");
- disp(" frequencies\n");
-
- prompt
- disp("")
- clc
- disp("Example #2, Consider the following set of discrete poles and");
- disp("zeros:\n")
- cmd = "sys2 = zp2sys([0.99258;0.99745],[0.99961;0.99242],1,0.001);";
- disp(cmd);
- eval(cmd);
- cmd = "sysout(sys2)";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system's bode plots, execute the following");
- disp("command:\n")
- cmd = "bode(sys2);";
- run_cmd;
- disp("Notice that the bode command is the same in both of the previous");
- disp("examples. The bode command is also the same for the continuous case.");
- disp("The function, dbode, is no longer used.");
-
- prompt
- disp("")
- clc
- disp("\nExample #3, Now consider the following state space system:\n");
- cmd = "sys3 = ss2sys([.857, .0011; 0, .99930],[1;1],[-.6318, .0057096],5.2, .001);";
- disp(cmd);
- eval(cmd);
- cmd = "sysout(sys3);";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system's bode plots, execute the following command:\n")
- cmd = "bode(sys3);";
- run_cmd;
- disp("\nAgain, notice that the bode command is the same regardless of the form");
- disp("of the system.");
- disp("\nSuppose the user is interested in the response of the system");
- disp("defined over the input frequency range of 1 - 1000 rad/s.\n");
- disp("First, a frequency vector is required. It can be created");
- disp("with the command:\n");
- cmd = "wrange = logspace(log10(1),log10(1000),100);";
- disp(cmd);
- eval(cmd);
- disp("\nThis creates a logrithmetically scaled frequency vector with");
- disp("100 values between 1 and 1000 rad/s\n");
- disp("Then, the bode command includes wrange in the input arguments");
- disp("like this:");
- cmd = "bode(sys3,wrange);";
- run_cmd;
- prompt;
-
- disp("")
- clc
- disp("\nExample #4, We will now examine a MIMO state-space system. Systems");
- disp("two and three will be grouped.");
- cmd = "[nn,nz] = sysdimensions(sys2);";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"out\",\"y_sys2\");";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"in\",\"u_sys2\");";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"st\",sysdefioname(nn+nz,\"x_sys2\"));";
- disp(cmd); eval(cmd);
- cmd = "sys_mimo = sysgroup(sys2,sys3);";
- disp(cmd); eval(cmd);
- cmd = "sysout(sys_mimo);";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system's bode plots, execute the following command:\n")
- cmd = "bode(sys_mimo);";
- run_cmd;
- prompt
-
- disp("\nThe bode plot of a single channel is viewed as follows:")
- cmd = "wrange = [];";
- disp(cmd)
- eval(cmd);
- cmd = "out = 1;";
- disp(cmd)
- eval(cmd);
- cmd = "in = 1;";
- disp(cmd)
- eval(cmd);
- cmd = "bode(sys_mimo,wrange,out,in);";
- run_cmd;
- disp("\nNotice that this bode plot is the same as the plot from example 2.");
- prompt
- closeplot
-
- elseif( k1 == 3 )
- help bode
- prompt
- endif
- endwhile
- elseif (j == 2)
- k2 = 0;
- disp("");
- while (k2 != 4)
- disp("\n");
- help nyquist
- prompt;
- disp("")
- clc;
-
- k2 = menu("Nyquist analysis (Nyquist)",...
- 'Continuous system nyquist analysis',...
- 'Discrete system nyquist analysis',...
- 'Mixed system nyquist analysis',...
- 'Return to frdemo menu');
-
- if( k2 == 1 )
- disp("")
- clc
- disp("\nContinuous system nyquist analysis\n");
- disp("Display Nyquist plots of a SISO system (nyquist)\n")
- disp("Example #1, Consider the following transfer function:\n")
- cmd = "sys1 = tf2sys(1, [1, 0.8, 1]);";
- disp(cmd);
- eval(cmd);
- disp("To examine the transfer function, use the command:");
- cmd = "sysout(sys1);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the open loop zeros and poles, use the command:");
- cmd = "sysout(sys1,""zp"");";
- run_cmd;
- disp("\nTo view the system""s nyquist plot, execute the following");
- disp("command:\n")
- cmd = "nyquist(sys1);";
- run_cmd;
- disp("\nIf the real and imaginary parts of the response are desired,");
- disp("use the following command:");
- disp("command: [R,I,w]=nyquist(sys1);\n");
- disp("If the user desires to evaluate the response in a certain");
- disp("frequency range, he may do so by entering the following:");
- disp("command: [M,P,w]=nyquist(num,den,wrange).\n")
- disp("wrange is a vector of frequencies that spans the desired");
- disp("viewing range.\n");
- disp("This will be illustrated in the third nyquist example.\n")
- disp("Variable Description:\n")
- disp("R => real part of response")
- disp("I => imaginary part of response")
- disp("w => frequencies that the transfer function was evaluated at")
- disp("sys1 => system data structure")
- disp("wrange => optional vector of frequencies")
- disp(" if wrange is entered in the argument list, the");
- disp(" system will be evaluated at these specific");
- disp(" frequencies\n")
- prompt
-
- disp("")
- clc
- disp("Example #2, Consider the following set of poles and zeros:\n")
- cmd = "sys2 = zp2sys([-1;-4],[-2+1.4142i;-2-1.4142i],1);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the poles and zeros, use the command:");
- cmd = "sysout(sys2)";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system""s nyquist plot, execute the following");
- disp("command:\n")
- cmd = "nyquist(sys2);";
- run_cmd;
- prompt
-
- disp("")
- clc
- disp("\nExample #3, Consider the following state space system:\n")
- cmd = "sys3 = ss2sys([0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; 0, 0, -20, -12],[0;0;0;1],[50, 100, 0, 0],0);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the state-space system, use the command:");
- cmd = "sysout(sys3)";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the poles and zeros, use the command:");
- cmd = "sysout(sys3,""zp"")";
- run_cmd;
- disp("\nTo view the system""s nyquist plot, execute the following");
- disp("commands:\n")
- cmd = "nyquist(sys3);";
- run_cmd;
- prompt
-
- disp("Example #3 (continued), If the user wishes to evaluate the");
- disp("system response over a desired frequency range, he must first");
- disp("create a frequency vector.\n")
- disp("For example, suppose the user is interested in the response");
- disp("of the system defined above over input frequency range of");
- disp("3 - 100 rad/s.\n")
- disp("A frequency vector can be created using the command:\n");
- cmd = "wrange = logspace(log10(3),log10(100),100);";
- disp(cmd);
- eval(cmd);
- disp("\nNyquist can be run again using the frequency vector as");
- disp("follows:\n")
- cmd = "nyquist(sys3,wrange);";
- run_cmd;
- prompt
-
- disp("")
- clc
- disp("Example #4, Nyquist can be used for MIMO systems if the system has");
- disp("an equal number of inputs and outputs. Otherwise, nyquist returns");
- disp("an error. To examine a MIMO system, systems 2 and 3 will be grouped");
- cmd = "[nn,nz] = sysdimensions(sys2);";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"out\",\"y_sys2\");";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"in\",\"u_sys2\");";
- disp(cmd); eval(cmd);
- cmd = "sys2 = syssetsignals(sys2,\"st\",sysdefioname(nn+nz,\"x_sys2\"));";
- disp(cmd); eval(cmd);
- cmd = "sys_mimo = sysgroup(sys2,sys3);";
- disp(cmd); eval(cmd);
- cmd = "sysout(sys_mimo);";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system's nyquist plot, execute the following command:\n")
- cmd = "nyquist(sys_mimo);";
- run_cmd;
- prompt
- disp("\nTo view the nyquist plots for selected channels, the command form changes:")
- cmd = "nyquist(sys_mimo,[],1,1);";
- run_cmd;
- disp("\nNotice that this bode plot is the same as the plot from example 2.");
- prompt
- closeplot
-
-
-
- elseif( k2 == 2 )
- disp("")
- clc
- disp("\nDiscrete system nyquist analysis\n");
- disp("Display Nyquist plots of a discrete SISO system (nyquist)\n")
- disp("We will first define a sampling time, T");
- cmd = "T = 0.01;";
- disp(cmd);
- eval(cmd);
- disp("\nExample #1, Consider the following transfer function:\n")
- cmd = "sys1 = tf2sys([2, -3.4, 1.5],[1, -1.6, 0.8],T);";
- disp(cmd);
- eval(cmd);
- disp("To examine the transfer function, use the command:");
- cmd = "sysout(sys1);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the open loop zeros and poles, use the command:");
- cmd = "sysout(sys1,""zp"")";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system""s nyquist plot, execute the following");
- disp("command:")
- cmd = "nyquist(sys1);";
- run_cmd;
- disp("To change the range used for the frequency, a frequency");
- disp("is needed. Suppose the user would like to examine the");
- disp("nyquist plot in the frequency range of 0.01 - 31.6 rad/s.");
- disp("\nThe frequency vector needed to do this is created with the");
- disp("command:");
- cmd = "wrange = logspace(-2,1.5,200);";
- disp(cmd);
- eval(cmd);
- disp("\nNyquist can be run again with this frequency vector");
- cmd = "nyquist(sys1,wrange);";
- run_cmd;
- disp("\nIf the real and imaginary parts of the response are desired,");
- disp("perform the following command:\n");
- disp("[R,I,w]=nyquist(sys,wrange)\n")
- disp("Variable Description:\n")
- disp("R => real part of response")
- disp("I => imaginary part of response")
- disp("w => frequencies that the transfer function was evaluated at")
- disp("sys => The system data structure");
- disp("wrange => optional vector of frequencies")
- disp(" if wrange is entered in the argument list, the");
- disp(" system will be evaluated at these specific");
- prompt
-
- disp("")
- clc
- disp("\nExample #2, Consider the following set of poles and zeros:\n")
- cmd = "sys2 = zp2sys([0.98025 + 0.01397i; 0.98025 - 0.01397i],[0.96079;0.99005],1,T);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the open loop zeros and poles, use the command:");
- cmd = "sysout(sys2)";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system's nyquist plot between the frequencies");
- disp("0.01 - 100 rad/s, execute the following commands:\n")
- cmd = "wrange = logspace(-2,2,100);";
- disp(cmd);
- eval(cmd);
- cmd = "nyquist(sys2,wrange);";
- run_cmd;
- prompt;
-
- disp("")
- clc
- disp("\nExample #3, Consider the following discrete state space");
- disp("system:\n");
- disp("This example will use the same system used in the third");
- disp("example in the continuous nyquist demo. First, that system");
- disp("will have to be re-entered useing the following commands:\n");
- cmd = "sys3 = ss2sys([0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; 0, 0, -20, -12],[0;0;0;1],[50, 100, 0, 0],0);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the state-space system, use the command:");
- cmd = "sysout(sys3)";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the poles and zeros, use the command:");
- cmd = "sysout(sys3,""zp"")";
- disp(cmd);
- eval(cmd);
- disp("\nTo convert the system to discrete time, we need a sampling");
- disp("time which can be entered like this:");
- cmd = "T = 0.01";
- disp(cmd);
- eval(cmd);
- disp("\nNow the command, c2d, is used to convert the system from");
- disp("continuous to discrete time, with the following command");
- cmd = "dsys3 = c2d(sys3,T);";
- run_cmd;
- disp("\nTo examine the new discrete state-space system, use the");
- disp("command");
- cmd = "sysout(dsys3);";
- disp(cmd);
- eval(cmd);
- disp("\nTo examine the new discrete poles and zeros, use the command:");
- cmd = "sysout(dsys3,""zp"")";
- disp(cmd);
- eval(cmd);
- disp("\nTo view the system's nyquist plot, execute the following");
- disp("commands:\n");
- cmd = "gset xrange [-4:2];";
- disp(cmd); eval(cmd);
- cmd = "gset yrange [-2.5:2.5];";
- disp(cmd); eval(cmd);
- cmd = "nyquist(dsys3);";
- run_cmd;
- disp("Notice that the asymptotes swamp out the behavior of the plot")
- disp("near the origin. You may use interactive nyquist plots")
- disp("to \"zoom in\" on a plot as follows:")
-
- cmd = "atol = 1;";
- disp(cmd)
- eval(cmd)
- cmd = "nyquist(dsys3,[],[],[],atol);";
- run_cmd
- prompt
-
-
- disp("")
- clc
- disp("MIMO SYSTEM: Nyquist cannot be used for discrete MIMO systems");
- disp("at this time.");
- ## cmd = "dsys_mimo = sysgroup(sys2,dsys3);";
- ## disp(cmd);
- ## eval(cmd);
- ## cmd = "sysout(dsys_mimo);";
- ## disp(cmd);
- ## eval(cmd);
- ## disp("\nTo view the system's nyquist plot, execute the following command:\n")
- ## cmd = "nyquist(dsys_mimo);";
- ## run_cmd;
- ## prompt
- ## disp("\nTo view the nyquist plots for selected channels, the command form changes:")
- ## cmd = "nyquist(dsys_mimo,[],1,1);";
- ## run_cmd;
- ## disp("\nNotice that this bode plot is the same as the plot from example 2.");
- prompt
- closeplot
-
-
- elseif( k2 == 3 )
- disp("\nMixed system nyquist analysis\n");
- disp("Nyquist exits with an error if it is passed a ""mixed"" system (one")
- disp("with both continuous and discrete states). Use c2d or d2c to")
- disp("convert the system to either pure digital or pure continuous form");
- endif
- endwhile
- elseif (j == 3)
- help nichols
- prompt
- endif
- endwhile
- endfunction
-