function [V,S,VS]=estorderls(y,u,nmax);
//function [V,S,VS]=estorderls(y,u,nmax);
//function for model order estimation. The ordre is estimated using least square methode,
//where the criterion has the following forme:
//                       1						      
//                 V(n)=---Min{|y(0) - R(n).theta|^2}  n=1,2,3,...nmax
//                       N
//and for n=0, V(n)=|y(0)^2|/N.
//inputs:
//y ... 	vector of acquired system outputs
//u ... 	vector of applied inputs
//nmax ...maximum order to be considered
//otputs:
//V...		normalized vector of the criterion (V=[V(0) V(1) V(2), ...]^T, where V(0)=1)
//S...		vector of penalty coefficients. Model complexity n is penalised.
//VS...	normalized vector of penalised criterion VS=V+S
//
//written by:  H. Prochazka, I.D. Landau
//7th june 2002

sz=size(y);
if(sz(2)~=1), y=y';end;//column
sz=size(u);
if(sz(2)~=1), u=u';end;

yz=[zeros(nmax-1,1);y];
uz=[zeros(nmax-1,1);u];

N=length(y)-nmax;
y0=y(nmax+1:N+nmax);   
V=sum(y0.^2)/N;
S=0;
VS=V;
for nn=1:nmax,
   //construction of R(nn)
   Rnny=[];
   Rnnu=[];
   for k=1:nn,
      yp=y(k:k+N-1);
      up=u(k:k+N-1);
      Rnny=[yp Rnny];
      Rnnu=[up Rnnu];
   end;
	y0=y(nn+1:N+nn);   
   Rnn=[Rnny Rnnu];
   theta_min=inv(Rnn'*Rnn)*Rnn'*y0;
   criter=sum((y0-Rnn*theta_min).^2)/N;
   V=[V;criter];
   penal=0.4*nn*log10(N)/N;
   S=[S;penal];
end;
//normalization
V=V/V(1);
VS=V+S;
VS=VS/VS(1);