function th = canstart(z,pobs,nu,oe)
%CANSTART Produces a multivariable theta-model with initial parameter estimates
%	THETA = canstart(Z,PSOBS,NU,OE)
%	
%	THETA: The resulting model, parametrized according to the pseudoobser-
%              vability indices and with initial parameter estimates calculated
%	       from data. 
%	Z :    The output-input data [y u], with each input/output component
%	       as a column vector.
%	PSOBS: The pseudo-observability indices that define the multivariable
%	       state-space structure. This is a row vector with as many entries
%	       as there are outputs. Basically, PSOBS(k) gives the number of
%	       delays of output # k that are used. The order of the system
%	       is the sum of the psobs-indices. See also HELP CANFORM.
%	NU:    The number of inputs.
%	OE:    If OE='oe', then the K-matrix is THETA is fixed to zero, so
%	       that an output-error model structure is obtained. Otherwise
%	       (default) K is initialized at a value that gives a stable
%	       predictor.
%
%	THETA is just an initial estimate of the system, and should be further
%	      improved by TH = pem(z,THETA);

%	L. Ljung 21-9-90,11-2-92
%	Copyright (c) 1990-92 by the MathWorks, Inc.
%	All Rights Reserved.

if nargin<4,oe='no';end
if nu==0, error('Sorry, this routine does not work for time series!'),end
ny=length(pobs);
[Ncap,nz]=size(z);
if nz~=nu+ny, error(['The number of inputs and outputs in the data vector is not'...
'\nconsistent with the number of observability indices and the number of inputs!']),end
mp=max(pobs);
r=cumsum(pobs);
stateind=[]; newstate=[];
for kl=1:ny
    for ni=0:pobs(kl)-1
        stateind=[stateind, kl+ni*ny];
    end	
    newstate=[newstate,kl+(ni+1)*ny];
end
nx=length(stateind);
th = iv4mv(z,[mp*ones(ny,ny),mp*ones(ny,nu),ones(ny,nu)]);
[a,b,c,d]=th2ss(th);
[a,b,c]=minreal(a,b,c,d);
[nxa,ndum]=size(a);
co=ctrb(a,b);ob=obsv(a,c); ob=[ob;ob((nxa-1)*ny+1:nxa*ny,:)*a];
hm=ob*co;
apar=(hm(newstate,:)/hm(stateind,:))';
bpar=(hm(stateind,1:nu))';
ms=canform(pobs,nu);
th=ms2th(ms,'d',[apar(:)' bpar(:)' zeros(1,ny*nx)]);
if oe=='oe',
    kind=[];
    for kx=1:nx,for ky=1:ny
        kind=[kind;kx,ky];
    end,end
    th=fixpar(th,'K',kind);
else	
    [a,b,c]=th2ss(th);
    k=real(dlqe(a,eye(nx),c,eye(nx),eye(ny)));
    k=(a*k)';
    th=ms2th(ms,'d',[apar(:)' bpar(:)' k(:)']);
end
th(2,7)=24;


