function [theta,mserror,yhat]=orthobases(ze,zv,psi)
%--------------------------------------------------------------------------
% Usage:  [theta,mserror,yhat]=orthobases(ze,zv,psi)
%--------------------------------------------------------------------------
% Author: Juan Carlos Gomez
% Last Modified: May 02, 2002
%--------------------------------------------------------------------------
% This function computes the Least Squares Estimates of the parameter vector
% 'theta' for the following model structure
%      
%       y(n)= G(q,theta)u(n) + v(n)
%
% where G(q,theta) is represented as a (truncated) series expansion
%
%		  G(q,theta)= SUM_(k=0)^(p-1) theta_k B_k (q)
%
% where B_k (q) are the Orthonormal Bases with Fixed Poles at the locations 
% given in vector 'psi'. The orthonormal bases are defined as
%
%                  sqrt(1 - |psi_0|^2)
%       B_0 (q) = ---------------------
%                       q - psi_0
%
%                  sqrt(1 - |psi_k|^2)                     1 - q psi_i^*
%       B_k (q) = --------------------- PROD_[i=0]^[k-1] ----------------
%                       q - psi_k                          q - psi_i
%       for k >= 1
%
%       psi=[psi_0,psi_1, ..., psi_{p-1)]     real stable poles
% 		  ze=[ye,ue]                            estimation data
%       zv=[yv,uv]                            validation data
%       theta                                 parameter estimation
%       mserror                               output mean square error
%       yhat                                  estimated output
%-----------------------------------------------------------------------
ue=ze(:,2);
uv=zv(:,2);
ye=ze(:,1);
yv=zv(:,1);
p=length(psi);
num1=sqrt(1-psi(1)^2);
den1=[1 -psi(1)];
num=num1;
den=den1;
phit=dlsim(num,den,ue);
for k=2:length(psi)
   numk=conv(num,[-psi(k-1) 1])*sqrt(1-psi(k)^2)/sqrt(1-psi(k-1)^2);
   denk=conv(den,[1 -psi(k)]);
   num=numk;
   den=denk;
   phik=dlsim(num,den,ue);
   phit=[phit,phik];
end
theta=phit\ye;
num=num1;
den=den1;
phitv=dlsim(num,den,uv);
for k=2:length(psi)
   numk=conv(num,[-psi(k-1) 1])*sqrt(1-psi(k)^2)/sqrt(1-psi(k-1)^2);
   denk=conv(den,[1 -psi(k)]);
   num=numk;
   den=denk;
   phikv=dlsim(num,den,uv);
   phitv=[phitv,phikv];
end
yhat=phitv*theta;
phi=phit.';
mserror=norm(yv-yhat)^2/length(yv);
plot([yv, yhat])
title('True (blue) and Estimated (green) Output - Validation data')
grid