Global Index (short | long) | Local contents | Local Index (short | long)
th=thss2th(eta,iy,noises)
THSS2TH Converts an internal TH(SS)-format to the standard multi-input-
single-output THETA-format. Auxiliary routine to th2ff and th2zp.
TH = thss2th(TH_SS,IY)
TH_SS: The model defined in the THETA(SS) format (See help thss)
IY: The output number in TH_SS to be considered as the
output in TH (The THETA-format handles only single-
output models). Default IY=1
TH: The corresponding model in the THETA-format including
the transformed covariance matrix (See help theta).
By default, a diagonal noise matrix (H) is assumed. To obtain models
for all noise sources, use TH=thss2th(TH_SS,IY,'noises');
Then the noise sources are treated and counted as additional inputs.
The tranformation of the covariance matrix is carried out using
the Gauss' approximation formula with numerical differentiation.
The stepsize in the differentiation is decided by the
m-file nuderst.
| This function calls | This function is called by |
|---|---|
function th=thss2th(eta,iy,noises)
% L. Ljung 10-2-90,11-11-93
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.4 $ $Date: 1997/12/02 03:42:47 $
if nargin<3, noises='no';end
if nargin<2, iy=1;end
[nr,nc]=size(eta);
T=eta(1,2);nu=eta(1,3);ny=eta(1,4);
sspm=getmfth(eta);
nd=eta(1,5);
[etapar,p,lambda]=th2par(eta);if length(etapar)==0,etapar=0;end
dt=nuderst(etapar);
arg=getargth(eta);
if any(eta(2,8)==[2 3 4]) & T<0,flag=1;else flag=0;end
if eta(2,8)==3 & T<0,flaga=1;else flaga=0;end
if any(eta(2,8)==[2 3]),Tmod=-1;else Tmod=T;end
[a,b,c,d,k,x0]=feval(sspm,etapar,Tmod,arg);
[ny,nx]=size(c);
bk=[b k];de=[d,eye(ny)];
if flag, [a,bk,c,d]=contth(a,bk,c,d,abs(T),flaga);
b=bk(:,1:nu);k=bk(:,nu+1:nu+ny);end
c=c(iy,:);
de=de(iy,:);
%if nu>0,d=d(iy,:);else d=[];end
Apol=poly(a);
[nx,nu]=size(b);
for kk=1:nu+ny
Bpol(kk,:)=poly(a-bk(:,kk)*c)+(de(kk)-1)*Apol;
bs=(abs(Bpol(kk,:))>eps);
Bpol(kk,:)=bs.*Bpol(kk,:);
end
if strcmp(noises,'no')
Cpol=Bpol(nu+iy,:);
Bpol=Bpol(1:nu,:);
lambb=lambda(iy,iy);
if T>0,ia=max(find(abs(Cpol)>eps));Cpol=Cpol(1:ia);end
else
Cpol=1;
lambb=0;
end
if T>0,ia=max(find(abs(Apol)>eps));Apol=Apol(1:ia);end
Ncap=getncap(eta);
th=poly2th(Apol,Bpol,Cpol,1,[],lambb,T,1);
if ~isempty(Ncap),
ndtemp=length(th2par(th));th(2,1)=lambb*(1+ndtemp/Ncap)/(1-ndtemp/Ncap);
else th(2,1)=lambb;
end
nu1=th(1,3);
if nu1>0
na=th(1,4);nb=th(1,5:4+nu1);nc=th(1,5+nu1);
ndd=th(1,6+nu1);nf=th(1,7+nu1:6+2*nu1); nk=th(1,7+2*nu1:6+3*nu1);
else na=th(1,4);nb=0;nc=th(1,5);ndd=th(1,6);nf=0;nk=0;
end
if ~isempty(p), if norm(p)>0
for kl=1:nd
drawnow
th1=etapar;th1(kl)=th1(kl)+dt(kl)/2;
th2=etapar;th2(kl)=th2(kl)-dt(kl)/2;
[A1,B1,C1,D1,K1,X1]=feval(sspm,th1,Tmod,arg);
[A2,B2,C2,D2,K2,X2]=feval(sspm,th2,Tmod,arg);
BK1=[B1 K1];BK2=[B2 K2];DE1=[D1 eye(ny)];DE2=[D2 eye(ny)];
if flag [A1,BK1,C1,D1]=contth(A1,[B1 K1],C1,D1,abs(T),flaga);
[A2,BK2,C2,D2]=contth(A2,[B2 K2],C2,D2,abs(T),flaga);
if nu>0,B1=BK1(:,1:nu);B2=BK2(:,1:nu);end
K1=BK1(:,nu+1:nu+ny);K2=BK2(:,nu+1:nu+ny);
end
Apol1=poly(A1);Apol2=poly(A2);C1=C1(iy,:);C2=C2(iy,:);
DE1=DE1(iy,:);DE2=DE2(iy,:);
% if nu>0,D1=D1(iy,:);D2=D2(iy,:);end
for kk=1:nu+ny
B1pol(kk,:)=poly(A1-BK1(:,kk)*C1)+(DE1(kk)-1)*Apol1;
B2pol(kk,:)=poly(A2-BK2(:,kk)*C2)+(DE2(kk)-1)*Apol2;
end
if strcmp(noises,'no')
C1pol=B1pol(nu+iy,:);
B1pol=B1pol(1:nu,:);
C2pol=B2pol(nu+iy,:);
B2pol=B2pol(1:nu,:);nu1=nu;
else
C1pol=1;C2pol=1;nu1=nu+ny;
end
% C1pol=poly(A1-K1(:,iy)*C1);C2pol=poly(A2-K2(:,iy)*C2);
if na>0, diff(kl,1:na)=(Apol1(2:na+1)-Apol2(2:na+1))/dt(kl);end
if nc>0, diff(kl,na+sum(nb)+1:na+nc+sum(nb))=...
(C1pol(2:nc+1)-C2pol(2:nc+1))/dt(kl);end
sb=na;sf=na+sum(nb)+nc+ndd;
for k=1:nu1
if nb(k)>0,diff(kl,sb+1:sb+nb(k))=...
(B1pol(k,nk(k)+1:nk(k)+nb(k))-B2pol(k,nk(k)+1:nk(k)+nb(k)))/dt(kl);
end
sb=sb+nb(k);
end
end
PTH=diff'*p*diff;
[d,nd]=size(diff');
if nd>0,th(4:d+3,1:d)=PTH;end
end,end