frdemo.m

similer program for matlab

## Copyright (C) 1996, 1998, 2000, 2004, 2005, 2007
##               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 3 of the License, 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, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} frdemo ()
## Octave Control Toolbox demo: Frequency Response demo.
## @end deftypefn

## Author: David Clem
## Created: August 15, 1994
## a s hodel: updated to match new order of ss2zp outputs
## J Ingram:  updated for system data structure format August 1996

function frdemo ()

  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=tf([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=zp(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=ss([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 = tf([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 = zp([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 = ss([.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