home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / TH2ZP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.0 KB  |  126 lines
  1. function [zepo,Kg]=th2zp(th,ku,ky,thresh)
  2. %TH2ZP    Computes zeros, poles, static gains and their standard deviations
  3. %
  4. %    [ZEPO,K] = th2zp(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 their standard deviations.
  8. %    Its first row contains integers with information about the column in
  9. %    question. See the manual.
  10. %
  11. %    The rows of K contain in order: the input/output number, the static 
  12. %     gain, and its standard deviation.
  13. %
  14. %    Both discrete and continuous time models are handled by TH2ZP
  15. %
  16. %    With  [ZEPO,K] = th2zp(TH,KU,KY) the zeros and poles associated with
  17. %    the input and output numbers given by the entries of the row vectors
  18. %    KU and KY, respectively, are computed.
  19. %    Default values: KU=[1:number of inputs], KY=[1:number of outputs].
  20. %    The noise e is then regarded as input # 0.
  21. %    The information is best displayed by ZPPLOT. zpplot(th2zp(TH),sd) is a
  22. %    possible construction.
  23.  
  24. %    L. Ljung 7-2-87,10-2-90
  25. %    Copyright (c) 1987-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 isthss(th),eval('[zepo,Kg]=zpsssd(th,ku,ky,thresh);'),return,end
  33. nu=th(1,3);T=gett(th);
  34. if isempty(thresh),thresh=100000;end
  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)<-1,error('Input indices outside # of inputs in theta')
  39. return,end
  40.  
  41. [nt,ntt]=size(th);
  42. Novar=0;
  43. % Sort out the orders of the polynomials
  44. na=th(1,4); nb=th(1,5:nu+4);nc=th(1,nu+5);nd=th(1,nu+6);
  45. nf=th(1,nu+7:2*nu+6);
  46. ns=2*nu+3;
  47. Ncum=cumsum(th(1,4:3+ns));
  48. [par,coP]=th2par(th);
  49. if isempty(coP) | norm(coP)==0,Novar=1;end
  50.  
  51. % set up orders for the noise case
  52. kzero=find(ku<=0);kku=ku;if length(kzero)>0,kku(kzero)=nu+1;end
  53. nb(nu+1)=nc;nf(nu+1)=nd;
  54. nnb=max(nb(kku));nnf=max(nf(kku));if nnb==nc,nnb=nnb+1;end
  55. nzp=max(nnb-1,nnf+na);
  56.  
  57. zepo=zeros(nzp+1,4*length(ku))+inf;
  58. ll=1:na;a=[1 th(3,ll)];aa=0;
  59. if na>0,
  60.    if Novar, tempP=[];else tempP=coP(ll,ll);end
  61.    Apoles=rotvar(a,tempP);
  62.    [aa,sl]=size(Apoles);
  63. end
  64. cc=1;
  65. for k=ku
  66.     if k>0,k1=k;else k1=501;end
  67.     zepo(1,cc:cc+3)=[k1 60+k1 20+k1 80+k1];Kg(1,(cc+3)/4)=k1;
  68.  
  69.     if na>0,zepo(2:aa+1,cc+2:cc+3)=Apoles;end
  70.  
  71.     if k>0,
  72.         llb=Ncum(k)+1:Ncum(k+1);
  73.     llf=Ncum(nu+2+k)+1:Ncum(nu+3+k);
  74.     b=th(3,llb);nbb=nb(k);nff=nf(k);
  75.     else
  76.         llb=Ncum(nu+1)+1:Ncum(nu+2);
  77.     llf=Ncum(nu+2)+1:Ncum(nu+3);
  78.     b=[1 th(3,llb)];nbb=nc;nff=nd;
  79.     end
  80.     f=[1 th(3,llf)];
  81.     if length(b)>1,
  82.     if Novar,tempP=[];else tempP=coP(llb,llb);end
  83.         Bzeros=rotvar(b,tempP);%Corr 89-07-30
  84.         [bb,sl]=size(Bzeros);
  85.     for kinf=1:sl
  86.             kind=find(abs(Bzeros(:,kinf))>thresh);
  87.         if ~isempty(kind),Bzeros(kind,kinf)=inf*ones(length(kind),1);end
  88.         end
  89.     zepo(2:bb+1,cc:cc+1)=Bzeros;
  90.     end
  91.     if length(llf)>0,
  92.     if Novar,tempP=[];else tempP=coP(llf,llf);end
  93.     Fpoles=rotvar(f,tempP);
  94.         [ff,sl]=size(Fpoles);zepo(aa+2:aa+1+ff,cc+2:cc+3)=Fpoles; 
  95.     end
  96.     if T>0
  97.         bst=sum(b);ast=sum(a);fst=sum(f);
  98.         if abs(ast*fst)>eps,
  99.            kgg=bst/ast/fst;
  100.            dKdth=[-ones(1,na)*kgg/ast ones(1,nbb)/ast/fst -ones(1,nff)*kgg/fst];
  101.         else kgg=inf;
  102.         end
  103.     else
  104.         if abs(a(length(a))*f(length(f)))>eps
  105.            kgg=b(length(b))/a(length(a))/f(length(f));
  106.            if na>0,dna=[zeros(1,na-1) 1];else dna=[];end
  107.            if nbb>0,dnb=[zeros(1,nbb-1) 1];else dnb=[];end
  108.            if nff>0,dnf=[zeros(1,nff-1) 1];else dnf=[];end
  109.            dKdth=[-dna*kgg/a(length(a)) dnb/a(length(a))/f(length(f)) -dnf*kgg/f(length(f))];
  110.         else kgg=inf;
  111.         end
  112.     end
  113.     Kg(2,(cc+3)/4)=kgg;
  114.     if kgg==inf,Kg(3,(cc+3)/4)=inf; 
  115.     else
  116.     if Novar, 
  117.        Kg(3,(cc+3)/4)=0;
  118.     else
  119.        Pk=dKdth*coP([ll llb llf],[ll llb llf])*dKdth';
  120.        Kg(3,(cc+3)/4)=sqrt(Pk);
  121.     end
  122.     end
  123.     cc=cc+4;
  124. end
  125. zepo(1,:)=sign(T)*zepo(1,:);
  126.