home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ZP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.4 KB  |  72 lines
  1. function [zepo,Kg]=zp(th,ku,ky,thresh)
  2. %ZP    Computes zeros, poles and static gains associated with a model
  3. %
  4. %    [ZEPO,K] = zp(TH)
  5. %
  6. %    For a model defined by TH (in the format described by HELP THETA)
  7. %    ZEPO is returned as the zeros and poles and K as the static gains
  8. %    of the model. The first rows in these matrices contain integers with
  9. %    information about the column in question. For the zeros and the gains
  10. %    this integer is simply the input number, while for the poles it is
  11. %    20 + the input number. The noise source e is in this context counted
  12. %    as input number 0. For multi-output models, 1000*(output # - 1) is
  13. %    added to these numbers.
  14. %
  15. %    With  [ZEPO,K] = zp(TH,KU,KY) the zeros and poles associated with
  16. %    the input and output numbers given by the entries of the row vectors
  17. %    KU and KY, respectively, are computed.
  18. %    Default values: KU=[1:number of inputs], KY=[1:number of outputs].
  19. %    Noise input is counted as input # 0.
  20. %    The zeros and poles can be plotted by ZPPLOT. zpplot(zp(TH)) is a
  21. %    possible construction.
  22. %    See also th2zp.
  23.  
  24. %    L. Ljung 10-1-86, revised 7-3-87,10-2-90.
  25. %    Copyright (c) 1986-90 by the MathWorks, Inc.
  26. %    All Rights Reserved.
  27.  
  28. % *** Set up default values ***
  29. if nargin<4, thresh=[];end
  30. if nargin<3, ky=[];end
  31. if nargin<2, ku=[];end
  32. if isempty(thresh),thresh=100000;end
  33. if isthss(th),eval('[zepo,Kg]=zpss(th,ku,ky,thresh);'),return,end
  34. nu=th(1,3);
  35. if nu==0,kudef=0;else kudef=1:nu;end
  36. if isempty(ku), ku=kudef;end
  37. if any(ku==-1), ku(find(ku==-1))=0;end
  38. if max(ku)>nu | min(ku)<0, error('Input indices outside # of inputs in theta')
  39. return,end
  40.  
  41.  
  42.    [a,b,c,d,f]=polyform(th);
  43.  
  44.  
  45.   if length(find(ku==0))>0
  46.        b(nu+1,1:length(c))=c;
  47.        f(nu+1,1:length(d))=d;
  48.   end
  49.  
  50.   [n1,nb]=size(b); [n1,nf]=size(f);
  51.  
  52.   nzp=max(nb-1,nf+length(a)-2);
  53.  
  54.  zepo=zeros(nzp+1,2*length(ku))+inf;
  55.  
  56.       cc=1;
  57.     for k=ku
  58.       if k<=0,k1=501;k=nu+1;else k1=k;end
  59.       zepo(1,cc:cc+1)=[k1 20+k1];Kg(1,(cc+1)/2)=k1;
  60.       dg=conv(a,f(k,:));
  61.       zb=roots(b(k,:)); nzb=length(zb);
  62.             kind=find(abs(zb)>thresh);
  63.         if ~isempty(kind),zb(kind)=inf*ones(length(kind),1);end
  64.       if nzb>0, zepo(2:nzb+1,cc)=zb;end
  65.       zg=roots(dg); nzg=length(zg);
  66.       if nzg>0, zepo(2:nzg+1,cc+1)=zg;end
  67.       [tmp1,tmp2]=size(b(k,:));
  68.       [md,nd] = size(dg);
  69.       Kg(2,(cc+1)/2)=(b(k,:)*ones(tmp1,tmp2)')/(dg*ones(md,nd)');
  70.       cc=cc+2;
  71.     end
  72.