function [Je,Jelms]=conerrorlms(rd0,rx,P,mu,w0,nite,nexp);
% Problema de predicción lineal con SD y LMS
% dibuja el contorno de la superficie de error y la evolucion del error con LMS, los coeficientes y la curva
% de aprendizaje
%
%	rd0 ..... potencia señal deseada
%  rx ...... función de autocorrelación señal datos
%  P ....... vector correlación cruzada
%  mu ...... paso de adaptación
%  w0 ...... vector inicial
%  nite .... numero de iteraciones
%  nexp .... numero de experimentos
%
%  Salidas
%
%  Je ...... curva de aprendizaje del SD
%  Jelms ... curva de aprendizaje del LMS
%
figure(1)
clf
Rx=toeplitz(rx);
[V,D]=eig(Rx);
D
V
fprintf(1,'Dispersion de autovalores %f\n',D(2,2)/D(1,1));
fprintf(1,'maximo mu ....\n');
fprintf(1,'         SD .. %f\n',2./D(2,2));
fprintf(1,'         LMS . %f\n',2./trace(Rx));
hopt=inv(Rx)*P
Jmin=rd0-P'*hopt;
d=filter(sqrt(Jmin),[1 -hopt'],randn(nite,1));
h1min=hopt(1)-2;
h1max=hopt(1)+2;
h2min=hopt(2)-2;
h2max=hopt(2)+2;
n1=1;
n2=1;
for h1=h1min:0.05:h1max
   for h2=h2min:0.05:h2max
      h=[h1 h2].';
      J(n1,n2)=Jmin+(h-hopt)'*Rx*(h-hopt);
      n2=n2+1;
   end
   n2=1;
   n1=n1+1;
end
h1=h1min:0.05:h1max;
h2=h2min:0.05:h2max;
contour(h1,h2,J,50);
hold on
plot(hopt(1),hopt(2),'*r');
%
%  Steepest Descent
%
want=w0;
w1(1)=w0(1);
w2(1)=w0(2);
w1lms(1)=w0(1);
w2lms(1)=w0(2);
plot(w1(1),w2(1),'og');
Je(1)=Jmin+(w0-hopt)'*Rx*(w0-hopt);
for n=2:nite
   w=want+mu*(P-Rx*want);
   w1(n)=w(1);
   w2(n)=w(2);
   plot(w(1),w(2),'og');
   want=w;
   Je(n)=Jmin+(w-hopt)'*Rx*(w-hopt);
end
%
% LMS con nexp realizaciones para calcular la curva de aprendizaje
%
Jelms=zeros(nite,1);
wlmsav=zeros(nite,2);
for i=1:nexp
   wantlms=w0;
   Jelms(1)=Jelms(1)+d(1).^2;
   x=filter([0 1],1,d);
	X=convm(x,2)';
	wlms=wantlms+mu*(X(:,1)*d(1));
   wantlms=wlms;
   wlmsav(1,:)=wlmsav(1,:)+wlms';
	for n=2:nite
   	e=d(n)-wlms'*X(:,n);
   	Jelms(n)=Jelms(n)+e.^2;
      wlms=wantlms+mu*(X(:,n)*e');
      wlmsav(n,:)=wlmsav(n,:)+wlms';
      if(i==1)
	   	w1lms(n)=wlms(1);
   	   w2lms(n)=wlms(2);
         plot(wlms(1),wlms(2),'om');
      end
      wantlms=wlms;
   end
   d=filter(1,[1 -hopt'],randn(nite,1));
   d=d./sqrt(mean(d.^2));
end
Jelms=Jelms./nexp;
wlmsav=wlmsav./nexp;
plot(w1,w2);
plot(w1lms,w2lms,'r');
plot(wlmsav(:,1),wlmsav(:,2),'c');
figure(2)
clf
plot(w1);
hold on
plot(w1lms,'r');
plot(ones(nite,1)*hopt(1));
plot(w2,'g');
plot(w2lms,'m');
plot(wlmsav(:,1),'y');
plot(wlmsav(:,2),'c');
plot(ones(nite,1)*hopt(2),'g');
grid on
figure(3)
clf
plot(Je);
hold on
plot(Jelms,'r');