home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / SS2ZP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.8 KB  |  112 lines
  1. function [z,p,k] = ss2zp(a,b,c,d,iu)
  2. %SS2ZP    State-space to zero-pole conversion.
  3. %    [Z,P,K] = SS2ZP(A,B,C,D,IU)  calculates the transfer function in
  4. %    factored form:
  5. %
  6. %                 -1           (s-z1)(s-z2)...(s-zn)
  7. %        H(s) = C(sI-A) B + D =  k ---------------------
  8. %                              (s-p1)(s-p2)...(s-pn)
  9. %    of the system:
  10. %        .
  11. %        x = Ax + Bu
  12. %        y = Cx + Du
  13. %
  14. %    from the single input IU.  The vector P contains the pole 
  15. %    locations of the denominator of the transfer function.  The 
  16. %    numerator zeros are returned in the columns of matrix Z with as 
  17. %    many columns as there are outputs y.  The gains for each numerator
  18. %    transfer function are returned in column vector K.
  19. %
  20. %    See also: ZP2SS,PZMAP,TZERO, and EIG.
  21.  
  22. %     J.N. Little 7-17-85
  23. %    Revised 3-12-87 JNL, 8-10-90 CLT, 1-18-91 ACWG.
  24. %    Copyright (c) 1985-92 by the MathWorks, Inc.
  25.  
  26. error(nargchk(4,5,nargin));
  27. error(abcdchk(a,b,c,d));
  28.  
  29. [nx,ns] = size(a);
  30.  
  31. if nargin==4,
  32.     if nx>0,
  33.         [nb,nu] = size(b);
  34.     else
  35.         [ny,nu] = size(d);
  36.     end
  37.     if (nu<=1), 
  38.         iu = 1;
  39.     else
  40.         error('IU must be specified for systems with more than one input.');
  41.     end
  42. end
  43.  
  44. % Remove relevant input:
  45. if ~isempty(b), b = b(:,iu); end
  46. if ~isempty(d), d = d(:,iu); end
  47.  
  48. % Trap gain-only models
  49. if nx==0&~isempty(d), z = []; p = []; k = d; return, end
  50.  
  51. % Do poles first, they're easy:
  52. p = eig(a);
  53.  
  54. k = []; 
  55. % Compute zeros and gains using transmission zero calculation
  56. % Check if Control Toolbox is on path, in which case use the 
  57. % more accurate tzero method. 
  58. if exist('tzreduce')
  59. % New code if new tzero code exists:
  60.     [ny,nu] = size(d);
  61.     z = [];
  62.     k = [];
  63.     for i=1:ny
  64.       [zi,gi] = eval('tzero(a,b,c(i,:),d(i,:))');
  65.       [mz,nz] = size(z);
  66.       nzi = length(zi);
  67.       if i==1,
  68.         z = zi;
  69.       else
  70.          z = [[z; inf*ones(max(0,nzi-mz),max(nz,1))], ...
  71.              [zi;inf*ones(max(0,mz-nzi),1)]];
  72.       end
  73.       k = [k;gi];
  74.     end
  75.  
  76. else
  77. % If Control Toolbox not on path then use generalized eigenvalues:
  78.  
  79. % Now try zeros, they're harder:
  80.     [no,ns] = size(c);
  81.     z = zeros(ns,no) + inf;         % Set whole Z matrix to infinities
  82. % Loop through outputs, finding zeros
  83.     for i=1:no
  84.         s1 = [a b;c(i,:) d(i)];
  85.         s2 = diag([ones(1,ns) 0]);
  86.         zv = eig(s1,s2);
  87.         % Now put NS valid zeros into Z. There will always be at least one
  88.         % NaN or infinity.
  89.         zv = zv(~isnan(zv) & finite(zv));
  90.         if ~isempty(zv)
  91.             z(1:length(zv),i) = zv;
  92.         end
  93.     end
  94. end
  95. % Now finish up by finding gains using Markov parameters
  96. if isempty(k)
  97.     k = d;  CAn = c; iter = 0;
  98.     while any(k==0) % do until all k's are finished
  99.         i = find(k==0);
  100.         if iter > ns
  101.             % Note: iter count is needed because when B is all zeros
  102.             % it does not otherwise converge.
  103.             k(i) = zeros(max(size(i)),1);
  104.             break
  105.         end
  106.         iter = iter + 1;
  107.         markov = CAn*b;
  108.         k(i) = markov(i);
  109.         CAn = CAn*a;
  110.     end
  111. end
  112.