Global Index (short | long) | Local contents | Local Index (short | long)
[zepo,Kg]=th2zp(th,ku,ky,thresh)
TH2ZP Computes zeros, poles, static gains and their standard deviations. [ZEPO,K] = TH2ZP(TH) For a model defined by TH (in the theta format. See also THETA.) ZEPO is returned as the zeros and poles and their standard deviations. Its first row contains integers with information about the column in question. See the manual. The rows of K contain in order: the input/output number, the static gain, and its standard deviation. Both discrete and continuous time models are handled by TH2ZP With [ZEPO,K] = TH2ZP(TH,KU,KY) the zeros and poles associated with the input and output numbers given by the entries of the row vectors KU and KY, respectively, are computed. Default values: KU=[1:number of inputs], KY=[1:number of outputs]. The noise e is then regarded as input # 0. The information is best displayed by ZPPLOT. ZPPLOT(TH2ZP(TH),sd) is a possible construction. See also TH2FF, TH2POLY, TH2TF, TH2SS, and ZP.
| This function calls | This function is called by |
|---|---|
function [zepo,Kg]=th2zp(th,ku,ky,thresh)
% L. Ljung 7-2-87,10-12-93
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.4 $ $Date: 1998/07/15 19:32:10 $
if nargin<1
disp('Usage: ZEPO = TH2ZP(TH)')
disp(' [ZEPO,K] = TH2ZP(TH,INPUTS,OUTPUTS)')
return
end
pars=th2par(th);
% *** Set up default values ***
if nargin<4,thresh=[];end
if nargin<3,ky=[];end
if nargin<2,ku=[];end
[dum,P,dum]=th2par(th);
if any(abs(imag(pars))>eps)
disp('This model has complex valued parameters.')
disp('The calculation of confidence regions is then not supported.')
disp('ZP is used instead.')
eval('[zepo,Kg]=zp(th,ku,ky,thresh);')
return
end
if isempty(P)|norm(P)==0
eval('[zepo,Kg]=zp(th,ku,ky,thresh);')
return
end
if isthss(th),
if length(th2par(th))>10
disp(str2mat('This may take a while.',...
'ZP will be faster, but does not evaluate standard deviations.'))
end
eval('[zepo,Kg]=zpsssd(th,ku,ky,thresh);'),return
end
nu=th(1,3);T=gett(th);
if isempty(thresh),thresh=100000;end
if nu==0, kudef=0;else kudef=1:nu;end
if isempty(ku), ku=kudef;end
if any(ku==-1),ku(find(ku==-1))=0;end
if max(ku)>nu | min(ku)<-1,error('Input indices outside # of inputs in theta')
return,end
[nt,ntt]=size(th);
Novar=0;
% Sort out the orders of the polynomials
na=th(1,4); nb=th(1,5:nu+4);nc=th(1,nu+5);nd=th(1,nu+6);
nf=th(1,nu+7:2*nu+6);
ns=2*nu+3;
Ncum=cumsum(th(1,4:3+ns));
[par,coP]=th2par(th);
if isempty(coP) | norm(coP)==0,Novar=1;end
% set up orders for the noise case
kzero=find(ku<=0);kku=ku;if length(kzero)>0,kku(kzero)=nu+1;end
nb(nu+1)=nc;nf(nu+1)=nd;
nnb=max(nb(kku));nnf=max(nf(kku));if nnb==nc,nnb=nnb+1;end
nzp=max(nnb-1,nnf+na);
zepo=zeros(nzp+1,4*length(ku))+inf;
ll=1:na;a=[1 th(3,ll)];aa=0;
if na>0,
if Novar, tempP=[];else tempP=coP(ll,ll);end
Apoles=rotvar(a,tempP);
[aa,sl]=size(Apoles);
end
cc=1;
for k=ku
if k>0,k1=k;else k1=501;end
zepo(1,cc:cc+3)=[k1 60+k1 20+k1 80+k1];Kg(1,(cc+3)/4)=k1;
if na>0,zepo(2:aa+1,cc+2:cc+3)=Apoles;end
if k>0,
llb=Ncum(k)+1:Ncum(k+1);
llf=Ncum(nu+2+k)+1:Ncum(nu+3+k);
b=th(3,llb);nbb=nb(k);nff=nf(k);
else
llb=Ncum(nu+1)+1:Ncum(nu+2);
llf=Ncum(nu+2)+1:Ncum(nu+3);
b=[1 th(3,llb)];nbb=nc;nff=nd;
end
f=[1 th(3,llf)];
if length(b)>1,
if Novar,tempP=[];else tempP=coP(llb,llb);end
Bzeros=rotvar(b,tempP);%Corr 89-07-30
[bb,sl]=size(Bzeros);
for kinf=1:sl
kind=find(abs(Bzeros(:,kinf))>thresh);
if ~isempty(kind),Bzeros(kind,kinf)=inf*ones(length(kind),1);end
end
zepo(2:bb+1,cc:cc+1)=Bzeros;
end
if length(llf)>0,
if Novar,tempP=[];else tempP=coP(llf,llf);end
Fpoles=rotvar(f,tempP);
[ff,sl]=size(Fpoles);zepo(aa+2:aa+1+ff,cc+2:cc+3)=Fpoles;
end
if T>0
bst=sum(b);ast=sum(a);fst=sum(f);
if abs(ast*fst)>eps,
kgg=bst/ast/fst;
dKdth=[-ones(1,na)*kgg/ast ones(1,nbb)/ast/fst -ones(1,nff)*kgg/fst];
else kgg=inf;
end
else
if abs(a(length(a))*f(length(f)))>eps
kgg=b(length(b))/a(length(a))/f(length(f));
if na>0,dna=[zeros(1,na-1) 1];else dna=[];end
if nbb>0,dnb=[zeros(1,nbb-1) 1];else dnb=[];end
if nff>0,dnf=[zeros(1,nff-1) 1];else dnf=[];end
dKdth=[-dna*kgg/a(length(a)) dnb/a(length(a))/f(length(f)) -dnf*kgg/f(length(f))];
else kgg=inf;
end
end
Kg(2,(cc+3)/4)=kgg;
if kgg==inf,Kg(3,(cc+3)/4)=inf;
else
if Novar,
Kg(3,(cc+3)/4)=0;
else
Pk=dKdth*coP([ll llb llf],[ll llb llf])*dKdth';
Kg(3,(cc+3)/4)=sqrt(Pk);
end
end
cc=cc+4;
end
zepo(1,:)=sign(T)*zepo(1,:);