% RLS 算法 
   clear 
   close all
 M=100;
 N=1024;
 w1=zeros(M,N);
 w2=zeros(M,N);
 ee=zeros(M,N);
 Lamda1 = 0.95 ;
 w10=zeros(M,N);
 w20=zeros(M,N);
 ee0=zeros(M,N);
 Lamda2 = 0.6 ; 
 Delta =10000;
 C=Delta*eye(2,2);
  
 a1=1.558;
 a2=-0.81;
 for m=1:M
 xu=randn(N,1);
 a(1)=xu(1);
 a(2)=a1*a(1)+xu(2);
   for n=3:N
   a(n)=xu(n)+a1*a(n-1)+a2*a(n-2);
   end
   for n=1:N
   x(n,1)=a(n);  %列向量
   end
   for n=3:N
     X=[x(n-1);x(n-2)];
     W=[w1(m,n-1);w2(m,n-1)];
     e(n)=x(n,1)-X'*W;
     u=X'*C*X;
     g=C*X/( Lamda1+u);
     W=W+g*e(n);
     w1(m,n)=W'*[1;0];
     w2(m,n)=W'*[0;1];
     C=1/Lamda1*(C-g*X'*C);
     
     ee(m,n)=e(n)^2;
   end
end
we1=zeros(1,N);
we2=zeros(1,N);
   for n=1:N
       for m=1:M
           we1(n)=we1(n)+w1(m,n);
           we2(n)=we2(n)+w2(m,n);
       end
       wee1(n)=we1(n)/M;
       wee2(n)=we2(n)/M;
   end 
   ee2=zeros(1,N);
   for n=1:N
       for m=1:M
         ee2(n)=ee2(n)+ee(m,n);
       end
       eee(n)=ee2(n)/M;
   end 
   %不同的lamda值
   C=Delta*eye(2,2);
for m=1:M
  xu=randn(N,1);
  a(1)=xu(1);
  a(2)=a1*a(1)+xu(2);
   for n=3:N
   a(n)=xu(n)+a1*a(n-1)+a2*a(n-2);
   end
   for n=1:N
   x(n,1)=a(n);  %列向量
   end
   for n=3:N
     X=[x(n-1);x(n-2)];
     W=[w10(m,n-1);w20(m,n-1)];
     e0(n)=x(n,1)-X'*W;
     u=X'*C*X;
     g=C*X/( Lamda2+u);
     W=W+g*e0(n);
     w10(m,n)=W'*[1;0];
     w20(m,n)=W'*[0;1];
     C=1/Lamda2*(C-g*X'*C);
     
     ee0(m,n)=e0(n)^2;
   end
 end
we10=zeros(1,N);
we20=zeros(1,N);
   for n=1:N
       for m=1:M
           we10(n)=we10(n)+w10(m,n);
           we20(n)=we20(n)+w20(m,n);
       end
       wee10(n)=we10(n)/M;
       wee20(n)=we20(n)/M;
   end 
   ee20=zeros(1,N);
   for n=1:N
       for m=1:M
         ee20(n)=ee20(n)+ee0(m,n);
       end
       eee0(n)=ee20(n)/M;
   end 
end
figure(1)
plot(w1(1,:),'r')
hold on
plot(w2(1,:),'b')
plot(wee1,'m')
plot(wee2,'k')
legend('w1(n)','w2(n)','E[w1(n)]','E[w2(n)]')
hold off
axis([0,1024,-2,2])
title('lamda1時的權系數變化曲線')
xlabel('n');
ylabel('w(n)');

figure(2)
subplot(3,1,1)
plot(20*log10(eee),'r')
hold on
plot(20*log10(ee(1,:)),'b')
legend('平均100次','單次運算')
hold off
axis([0,1024,-40,20])
title('平方誤差的變化曲線')
xlabel('n');
ylabel('e^2');

figure(3)
plot(w1(1,:),'r')
hold on
plot(w2(1,:),'b')
plot(w10(1,:),'m')
plot(w20(1,:),'k')
legend('lamda=0.95時w1(n)','lamda=0.95時w2(n)','lamda=0.6時w1(n)','lamda=0.6時w2(n)')
hold off
axis([0,1024,-2,2])
title('不同lamda時的權系數變化曲線比較')
xlabel('n');
ylabel('w(n)');

figure(4)
subplot(3,1,1)
plot(20*log10(eee),'r')
hold on
plot(20*log10(eee0),'b')
hold off
axis([0,1024,0,20])
title('lamda不同時的學習曲線比較')
legend('lamda=0.95','lamda=0.6')
xlabel('n');
ylabel('e^2');