function [B,A]=nrls(y,u,na,nb,d);
//function [B,A]=nrls(y,u,na,nb,d);
//linear digital model identification using non-recursive least square methode. The least square
//problem is resolved by QR-triangularization. The identified model is digital transfer 
//function of the following forme:
//                              z^-d (b_1z^-1 + b_2z^-2 + ... + b_nbz^-nb)
//                     G(z^-1)= ------------------------------------------
//                               1 + a_1z^-1 + a_2z^-2 + ... + a_naz^-na
//
//where na is the order of denominator, nb is the order of numerator and d is system's delay.
//inputs:
//y ...	vector of system outputs
//u ...	vector of system inputs
//na ...	order of model denominator
//nb ...	order of model numerator
//d ...	model's delay
//outputs:
//B ...	vector of numerator coefficients B=[0 0 ... 0 b_1 b_2 ... b_nb] (zeros is d+1)
//A ...	vector of denominator coefficients A=[1 a_1 a_2 ... a_na]
//
//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;

nmax=max(na,nb+d);
ntheta=na+nb;
N=length(y)-nmax-1;

//y vector
Rnny=[];
for k=2:na+1,
   yp=y(k:k+N-1);
   Rnny=[yp Rnny];
end;
//u vector
Rnnu=[];
for k=d+2:nb+d+1,
   up=u(k:k+N-1);
   Rnnu=[up Rnnu];
end;
y0=y(1:N);   
Rnn=[Rnny Rnnu];
//QR-triangularization
[Q,Res]=qr(Rnn);
Y0=Q*y0;
Res1=Res(1:2*ntheta,:);
Y01=Y0(1:2*ntheta);
theta = Res\(Res'\(Rnn'*y0));
//refinement
ther=y0 - Rnn*theta;
er=Res\(Res'\(Rnn'*ther));
theta=theta + er;

A=[1 theta(1:na)'];
B=[zeros(1,d) 0 theta(na+1:na+nb)'];
