home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ZPSS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.8 KB  |  57 lines
  1. function [zepo,k] = zpss(th,nnu,nny,thresh)
  2. %ZPSS    An auxiliary function to ZP
  3. %
  4. %    [zepo,k] = zpss(th,nnu,nny)
  5. %
  6.  
  7. %    L. Ljung 10-2-90,11-11-91
  8. %    Copyright (c) 1990 by the MathWorks, Inc.
  9. %    All Rights Reserved.
  10.  
  11. if ~exist('ss2zp'),error('This routine requires the SIGNAL PROCESSING TOOLBOX'),end
  12. [a,b,c,d,km]=eta2ss(th); [nx,nu]=size(b);[ny,nx]=size(c);
  13. if any(nnu>nu),error('There are not that many inputs in the model!'),end
  14. if any(nnu<-ny),error('There are not that many noise sources in the model!'),end
  15. if any(nny>ny) | any(nny<1)
  16.    error('There are not that many outputs in the model!'),end
  17. T=gett(th);
  18. if isempty(nnu),nnu=1:nu;end
  19. if isempty(nny),nny=1:ny;end
  20. lny=length(nny);
  21. if any(nnu==0),
  22.     nnu=[-ny:-1,nnu(find(nnu>0))];
  23. end
  24. if length(nnu)>0
  25. zepo=[];k=[];
  26. for ku=nnu
  27.     G1=zeros(1+nx,2*lny)+inf;ka1=zeros(2,lny);
  28.     if ku<0, btemp=km;dtemp=eye(ny);ku1=abs(ku);ku=500+ku1;
  29.     else btemp=b;dtemp=d;ku1=ku;end
  30.     [z1,p1]=ss2zp(a,btemp,c(nny,:),dtemp(nny,:),ku1);
  31.     [nzr,nzc]=size(z1);
  32.     if nzc<lny,  %This fix is due to a bug in ss2zp
  33.        z1=inf*ones(nx,lny);kcount=1;
  34.        for kky=nny 
  35.            z1t=ss2zp(a,btemp,c(kky,:),dtemp(kky,:),ku1);
  36.            if ~isempty(z1t),z1(1:length(z1t),kcount)=z1t;end,kcount=kcount+1;
  37.        end
  38.     end
  39.     [k1,ph1]=trfsaux(a,btemp,c(nny,:),dtemp(nny,:),ku1,eps,T);
  40.     k1=k1.*(1-2*abs(round(rem(ph1/180,2))));
  41.    [rz1,cz1]=size(z1);
  42.    for kinf=1:cz1
  43.       kind=find(abs(z1(:,kinf))>thresh);
  44.       if ~isempty(kind),z1(kind,kinf)=inf*ones(length(kind),1);end
  45.    end
  46.    G1(1,1:2:2*lny)=1000*(nny-1)+ku;
  47.    G1(1,2:2:2*lny)=1000*(nny-1)+20+ku;ka1(1,1:lny)=1000*(nny-1)+ku;
  48.    [nzr,nzc]=size(z1);
  49.    G1(2:nzr+1,1:2:2*lny)=z1;
  50.    G1(2:length(p1)+1,2:2:2*lny)=p1*ones(1,lny);
  51.    ka1(2,1:lny)=k1;
  52.    k=[k ka1];
  53.    zepo=[zepo G1];
  54. end,end
  55. zepo(1,:)=sign(T)*zepo(1,:);
  56.  
  57.