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