← Back to branch summary

~ubuntu-branches/ubuntu/quantal/octave-control/quantal

« back to all changes in this revision

Viewing changes to inst/__bodquist__.m

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2009-05-24 12:05:04 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20090524120504-seb43rd9wxtp17ka
Tags: 1.0.10-1
* New upstream release
* debian/patches/no-check-failure-acker-place.diff: Get around precision
  problems when testing functions acker and place
* debian/rules:
  + Include patchsys-quilt.mk
  + Remove unneeded documentation sources from /usr/share/doc
* debian/control:
  + (Build-Depends): Add quilt
  + (Standards-Version): Bump to 3.8.1 (no changes needed)
  + (Depends): Add ${misc:Depends}
  + (Vcs-Git, Vcs-Browser): Adjust to new Git repository
* debian/source.lintian-overrides: Add file for overriding Lintian
  warnings build-depends-without-arch-dep
* debian/copyright:
  + Use DEP5 URL in Format-Specification
  + Use separate License stanzas for instructing about the location of
    the different licenses used in the package
* debian/README.source: Add file explaining the quilt patch system, as
  required by the Policy

Show diffs side-by-side

added added

removed removed

Lines of Context:
inst/__bodquist__.m
 
1
## Copyright (C) 1996, 1998, 2000, 2002, 2003, 2004, 2005, 2006, 2007
 
2
##               Auburn University.  All rights reserved.
 
3
##
 
4
##
 
5
## This program is free software; you can redistribute it and/or modify it
 
6
## under the terms of the GNU General Public License as published by
 
7
## the Free Software Foundation; either version 3 of the License, or (at
 
8
## your option) any later version.
 
9
##
 
10
## This program is distributed in the hope that it will be useful, but
 
11
## WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
13
## General Public License for more details.
 
14
##
 
15
## You should have received a copy of the GNU General Public License
 
16
## along with this program; see the file COPYING.  If not, see
 
17
## <http://www.gnu.org/licenses/>.
 
18
 
 
19
## Undocumented internal function.
 
20
 
 
21
## -*- texinfo -*-
 
22
## @deftypefn {Function File} {[@var{f}, @var{w}, @var{rsys}] =} __bodquist__ (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx})
 
23
## Used internally by @command{bode}, @command{nyquist}; compute system frequency response.
 
24
##
 
25
## @strong{Inputs}
 
26
## @table @var
 
27
## @item sys
 
28
## input system structure
 
29
## @item w
 
30
## range of frequencies; empty if user wants default
 
31
## @item out_idx
 
32
## @itemx in_idx
 
33
## names or indices of output/input signal names; empty if user wants all
 
34
## @item rname
 
35
## name of routine that called __bodquist__ ("bode", "nyquist", or "nichols")
 
36
## @end table
 
37
## @strong{Outputs}
 
38
## @table @var
 
39
## @item w
 
40
## list of frequencies
 
41
## @item f
 
42
## frequency response of sys; @math{f(ii) = f(omega(ii))}
 
43
## @item rsys
 
44
## system with selected inputs and outputs
 
45
## @end table
 
46
##
 
47
## @code{bode}, @code{nichols}, and @code{nyquist} share the same 
 
48
## introduction, so the common parts are
 
49
## in __bodquist__.  It contains the part that finds the number of arguments,
 
50
## determines whether or not the system is @acronym{SISO}, and computes the frequency
 
51
## response.  Only the way the response is plotted is different between the
 
52
## these functions.
 
53
## @end deftypefn
 
54
 
 
55
function [f, w, rsys] = __bodquist__ (sys, w, outputs, inputs, rname)
 
56
 
 
57
  ## check number of input arguments given
 
58
  if (nargin != 5)
 
59
    print_usage ();
 
60
  endif
 
61
 
 
62
  ## check each argument to see if it's in the correct form
 
63
  if (! isstruct (sys))
 
64
    error ("sys must be a system data structure");
 
65
  endif
 
66
 
 
67
  ## let __freqresp__ determine w if it's not already given
 
68
  USEW = freqchkw (w);
 
69
 
 
70
  ## get initial dimensions (revised below if sysprune is called)
 
71
  [nn, nz, mm, pp] = sysdimensions (sys);
 
72
 
 
73
  ## check for an output vector and to see whether it`s correct
 
74
  if (! isempty (outputs))
 
75
    if (isempty (inputs))
 
76
      inputs = 1:mm;                    # use all inputs
 
77
      warning ("%s: outputs specified but not inputs", rname);
 
78
    elseif (is_signal_list (inputs) || ischar (inputs))
 
79
      inputs = sysidx (sys, "in", inputs);
 
80
    endif
 
81
    if (is_signal_list (outputs) || ischar (outputs))
 
82
      outputs = sysidx (sys, "out", outputs);
 
83
    endif
 
84
    sys = sysprune (sys, outputs, inputs);
 
85
    [nn, nz, mm, pp] = sysdimensions (sys);
 
86
  endif
 
87
 
 
88
  ## for speed in computation, convert local copy of
 
89
  ## SISO state space systems to zero-pole  form
 
90
  if (is_siso (sys) && strcmp (sysgettype (sys), "ss"))
 
91
    [zer, pol, k, tsam, inname, outname] = sys2zp (sys);
 
92
    sys = zp (zer, pol, k, tsam, inname, outname);
 
93
  endif
 
94
 
 
95
  ## get system frequency response
 
96
  [f, w] = __freqresp__ (sys, USEW, w);
 
97
 
 
98
  phase = arg(f)*180.0/pi;
 
99
 
 
100
  if (! USEW)
 
101
    ## smooth plots
 
102
    pcnt = 5;           # max number of refinement steps
 
103
    dphase = 5;         # desired max change in phase
 
104
    dmag = 0.2;         # desired max change in magnitude
 
105
    while (pcnt)
 
106
      pd = abs (diff (phase));                    # phase variation
 
107
      pdbig = find (pd > dphase);
 
108
 
 
109
      ## relative variation
 
110
      lp = length (f);
 
111
      lp1 = lp-1;
 
112
 
 
113
      fd = abs (diff (f));
 
114
      fm = max (abs ([f(1:lp1); f(2:lp)]));
 
115
      fdbig = find (fd > fm/10);
 
116
 
 
117
      bigpts = union (fdbig, pdbig);
 
118
 
 
119
      if (isempty (bigpts))
 
120
        pcnt = 0;
 
121
      else
 
122
        pcnt = pcnt - 1;
 
123
        wnew = [];
 
124
        crossover_points = find (phase(1:lp1).*phase(2:lp) < 0);
 
125
        pd(crossover_points) = abs (359.99+dphase - pd(crossover_points));
 
126
        np_pts = max (3, ceil(pd/dphase)+2);              # phase points
 
127
        nm_pts = max (3, ceil(log(fd./fm)/log(dmag))+2);  # magnitude points
 
128
        npts = min (5, max(np_pts, nm_pts));
 
129
 
 
130
        w1 = log10 (w(1:lp1));
 
131
        w2 = log10 (w(2:lp));
 
132
        for ii = bigpts
 
133
          if (npts(ii))
 
134
            wtmp = logspace (w1(ii), w2(ii), npts(ii));
 
135
            wseg(ii,1:(npts(ii)-2)) = wtmp(2:(npts(ii)-1));
 
136
          endif
 
137
        endfor
 
138
        wnew = vec(wseg)'; # make a row vector
 
139
        wnew = wnew(find (wnew != 0));
 
140
        wnew = sort (wnew);
 
141
        wnew = create_set (wnew);
 
142
        if (isempty (wnew))   # all small crossovers
 
143
          pcnt = 0;
 
144
        else
 
145
          ## get new freq resp points, combine with old, and sort.
 
146
          [fnew, wnew] = __freqresp__ (sys, 1, wnew);
 
147
          w = [w, wnew];
 
148
          f = [f, fnew];
 
149
          [w, idx] = sort (w);
 
150
          f = f (idx);
 
151
          phase = arg(f)*180.0/pi;
 
152
        endif
 
153
      endif
 
154
    endwhile
 
155
  endif
 
156
 
 
157
  ## ensure unique frequency values
 
158
  [w, idx] = sort (w);
 
159
  f = f(idx);
 
160
 
 
161
  w_diff = diff (w);
 
162
  w_dup = find (w_diff == 0);
 
163
  w_idx = complement (w_dup, 1:length(w));
 
164
  w = w(w_idx);
 
165
  f = f(w_idx);
 
166
 
 
167
  ## set rsys to pruned system
 
168
  rsys = sys;
 
169
 
 
170
endfunction