function [model, fit, yhat]=orthobasesid(datae,datav,poles)
%--------------------------------------------------------------------------
% Usage:  [model, fit, yhat]=orthobasesid(datae,datav,poles)
%         [model, fit, yhat]=orthobasesid(datae,datav)
%--------------------------------------------------------------------------
% Author: Juan Carlos Gomez
% Last Modified: April 07, 2006
%--------------------------------------------------------------------------
% This function returns an 'idmodel' obtained by computing 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 'poles'. 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
%
%       poles=[psi_0,psi_1, ..., psi_{p-1)]     real stable poles
% 		datae                                   estimation data (iddata)
%       datav                                   validation data (iddata)
%       model                                   estimated model (idmodel)
%       fit                                     model fit (validation data)
%       yhat                                    simulated model output
% When the third input argument is ommitted ('poles'), initial guesses for the poles
% are computed using the 'n4sid' algorithm.
%-----------------------------------------------------------------------
ue=datae.U;
uv=datav.U;
ye=datae.Y;
yv=datav.Y;
Ts=datae.Ts;
if nargin==3
p=length(poles);
num1=sqrt(1-poles(1)^2);
den1=[1 -poles(1)];
num=num1;
den=den1;
phit=dlsim(num,den,ue);
for k=2:length(poles)
   numk=conv(num,[-poles(k-1) 1])*sqrt(1-poles(k)^2)/sqrt(1-poles(k-1)^2);
   denk=conv(den,[1 -poles(k)]);
   num=numk;
   den=denk;
   phik=dlsim(num,den,ue);
   phit=[phit,phik];
end
theta=phit\ye;
[A,B,C,D] = ob2ss(poles,theta);
model=idss(A,B,C,D,zeros(size(A,1),1),zeros(size(A,1),1),Ts);
[yy, fit]=compare(datav,model);
compare(datav,model)
yh=yy{1};
yhat=yh.Y;
else 
modelo=n4sid(datae,'best','focus','prediction','n4weight','cva');
[AA,BB,CC,DD]=ssdata(modelo);
%poles=abs(eig(AA)).*sign(real(eig(AA)));
poles=real(eig(AA));
p=length(poles);
num1=sqrt(1-poles(1)^2);
den1=[1 -poles(1)];
num=num1;
den=den1;
phit=dlsim(num,den,ue);
for k=2:length(poles)
   numk=conv(num,[-poles(k-1) 1])*sqrt(1-poles(k)^2)/sqrt(1-poles(k-1)^2);
   denk=conv(den,[1 -poles(k)]);
   num=numk;
   den=denk;
   phik=dlsim(num,den,ue);
   phit=[phit,phik];
end
theta=phit\ye;
[A,B,C,D] = ob2ss(poles,theta);
model=idss(A,B,C,D,zeros(size(A,1),1),zeros(size(A,1),1),Ts);
[yy, fit]=compare(datav,model);
compare(datav,model)
yh=yy{1};
yhat=yh.Y;
end

