clear all;
k=0.5;
pi=3.1415;
N=10;
kn=100*N;
m_carlo=50;
m_ignore=10;
zao_white=1;
u_forget=[0.95,0.99,0.999];
dir_signal=0;
dir_gan=0.15;
%snr_db=[0];
snr_db=[20];
inr_db=10;

v = zeros(N,1);
u = zeros(N,1);
n = zeros(N,1);
Rs = zeros(N,N);
Ru = zeros(N,N);

inr_fact=sqrt(10^(inr_db/10));

for i = 1:N
   v(i) = exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*k*sin(dir_signal));
   u(i) = exp(sqrt(-1)*(-(N+1)/2+i)*2*pi*k*sin(dir_gan));
end 

for i_snr=1:1:length(snr_db)
   snr_fact=sqrt(10^(snr_db(i_snr)/10));
   signal_power=snr_fact.*exp(j*2*pi*rand);
   gan_power=inr_fact.*exp(j*2*pi*rand);
   Rs=signal_power*v*v'*signal_power';
   Ru=gan_power*u*u'*gan_power'+eye(N);
   R=Rs+Ru;
   
   for j_monte=1:1:m_carlo
     j_monte 
     for u_f=1:length(u_forget)
       P_rls = zeros(N,N);     
       P_rls = (1/zao_white)*eye(N);
       san = zeros(1,kn+1);
       san(1,1) = 1/(v'*P_rls*v);
       W_rls = v/N;
       g_rls = zeros(N,1);
       for temp2=1:1:kn
     
          signal_rls = snr_fact.*exp(j*2*pi*rand);
          gan_rls = inr_fact.*exp(j*2*pi*rand);
          n = (randn(N,1)+j*randn(N,1))/(sqrt(2));
          x = signal_rls*v+gan_rls*u+n;
          g_temp = 1+(1/u_forget(u_f))*x'*P_rls*x;
          g_rls = (1/u_forget(u_f))*P_rls*x/g_temp;
          P_rls = (1/u_forget(u_f))*P_rls-(1/u_forget(u_f))*g_rls*x'*P_rls;
          san(1,temp2+1) = 1/(v'*P_rls*v);
          W_rls = san(1,temp2+1)/(u_forget(u_f)*san(1,temp2))*(eye(N)-g_rls*x')*W_rls;
          if temp2>m_ignore
             diedai(temp2-m_ignore)=temp2-m_ignore;
             snr_temp = abs((W_rls'*Rs*W_rls)/(W_rls'*Ru*W_rls));
             if u_f==1
               if j_monte==1      
                  snr_vary(i_snr,temp2-m_ignore) = 10*log10(snr_temp);
               else 
                  snr_vary(i_snr,temp2-m_ignore) = ((j_monte-1)*snr_vary(i_snr,temp2-m_ignore)+10*log10(snr_temp))/j_monte;
               end
           elseif u_f==2 
               if j_monte==1
                  snr_varyn(i_snr,temp2-m_ignore) = 10*log10(snr_temp);
               else
                  snr_varyn(i_snr,temp2-m_ignore) = ((j_monte-1)*snr_varyn(i_snr,temp2-m_ignore)+10*log10(snr_temp))/j_monte;
              end
          else 
               
               if j_monte==1
                  snr_varynx(i_snr,temp2-m_ignore) = 10*log10(snr_temp);
               else
                  snr_varynx(i_snr,temp2-m_ignore) = ((j_monte-1)*snr_varynx(i_snr,temp2-m_ignore)+10*log10(snr_temp))/j_monte;
               end
             end  
          end
      
          if temp2>m_ignore
             W_l = inv(R)*v*inv(v'*inv(R)*v);
             snr_l_temp = abs((W_l'*Rs*W_l)/(W_l'*Ru*W_l));
             snr_l(i_snr,temp2-m_ignore) = 10*log10(snr_l_temp);
          end  

       end
      end
   end
end

figure(105);
wyy=10:length(diedai)+9;
wy=log10(wyy);

semilogx(10.^wy,snr_vary(1,:),'r-');
hold on;
semilogx(10.^wy,snr_varyn(1,:),'g:');
hold on;
semilogx(10.^wy,snr_varynx(1,:),'b-.');
hold on;
semilogx(10.^wy,snr_l(1,:),'c*');
hold on;
hold off;
title('SNR=20db,INR=10db');
xlabel('The iterative number');
ylabel('Average SINR(db)');
legend('yiwangyinzi 0.95','yiwangyinzi 0.99','yiwangyinzi 0.999','optimum'); 
grid on;