% Show the convergence of the parameter for fixed pole.
%
%       INPUT 
%   y: output of the system 
%   u: input of the system; white noise.
%   P: pole location
%
%       INTERNAL VARIABLES
%   N = number of blocks 
%   S = number of samples
%   T = number of true systems
%

N = size(P); N=N(1,1)/2;  
T = size(y); S=T(1,1);  T=T(1,2);

for n = 1:N,              % find c and b for the complex poles
 c(n) = - real(P(n*2))^2 - imag(P(n*2))^2;    
 b(n) = - 2*real(P(n*2))/(c(n) - 1);
end

n = 1;
den1 = [1 b(n)*(c(n)-1) -c(n)];
L(:,1) = dlsim([0 1 -b(n)],den1,u); 
 
num1 = [0 0 1];
L(:,2) = dlsim(num1,den1,u); 
   
num1 = 1;
if N > 1,
 for n=2:N,
   num2 = conv(num1,[-c(n-1) b(n-1)*(c(n-1)-1) 1]);
   den2 = conv(den1,[1 b(n)*(c(n)-1) -c(n)]);
   L(:,n*2-1)=dlsim(conv([0 1 -b(n)],num2),den2,u); 
   L(:,n*2)=dlsim(conv([0 0 1],num2),den2,u); 
   num1 = num2;
   den1 = den2;
 end
end

theta_hat = zeros(2*N,1);
size_L = size(L)
disp('Recursive estimation')
for s = 1:S,
 theta1(s) = theta_hat(1,1);
 error(s) = y(s) - L(s,:) * theta_hat;
 theta_hat = theta_hat + .03*L(s,:)'*error(s)/(.001+L(s,:)*L(s,:)');
end
plot(log(error.*error))
