// [theta] = LeastSquares(na,nb,nk,u,y,mu)
// Batch Least Square Identification for ARX models
// A(z)Y(z) = z^-nk B(z)U(z) + zita(z)
//
// input:
// 
// na - order of A(z)
// nb - order of B(z)
// nk - input-output delay
// u  - input vector
// y  - output vector
// mu - delay factor (1 for LS)
//
// output:
//
// theta - parameters estimation of dimension na+nb
//
// Feb 2008
function [theta] = LeastSquares(na,nb,nk,u,y,mu)

	// u,y should be column vector
	[d1,d2]=size(u);if d2>d1, u=u'; end
	[d1,d2]=size(y);if d2>d1, y=y'; end
	
	n = length(y);          // number of samples
	m = max(na,nb+nk-1);    // m is the first non-zero output of the predictor
	
	PHI = [];               // stacked regressor
	Mu = [];                // stacked weight vector
	for i=(m+1):n
		phi = [-y(i-1:-1:i-na)', u(i-nk:-1:i-nk-nb+1)'];   // i-sample regressor
		PHI = [PHI; phi];
		Mu = [Mu; mu^(n-i)];
	end
	W = diag(Mu);           // stacked weight matrix
	if rank(PHI'*W*PHI)<(na+nb)
		disp(['singular PHI''PHI matrix'])
		theta=zeros(na+nb,1);
		abort
	end
	theta = inv(PHI'*W*PHI)*PHI'*W*y(m+1:n); // pseudoinverse
	
	s = svd(PHI);
	disp([' regressor condition number : ' string(s(1)/s(length(s)))])
	disp([' norm of PHI''PHI           : ' string(norm(PHI'*PHI))])
	
	// display of the identified theta
	if mu==1
		disp([' theta (LS) = ' string(theta')])
	else
		disp([' theta (WLS) = ' string(theta')])
	end
	
	disp(' poles of the transfer function : '); 
	disp(roots([1; theta(1:na)]))

endfunction