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);