clear all;

k=0.5;
N=10;
kn=100*N;
m_carlo=50;
m_ignore=10;
dir_signal=0;
dir_gan=0.15;

snr_db=[0,10,20];
inr_db=10;

inr_fact=sqrt(10^(inr_db/10));
v = zeros(N,1);
u = zeros(N,1);
n = zeros(N,1);

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 
LNR=[-10 0 10];
plus=(10.^(LNR/10));

Rx=zeros(N,N);
Ru=zeros(N,N);
Rs=zeros(N,N);
Wx=zeros(N,1);

snr_vary1=zeros(1,kn-m_ignore);
snr_vary2=zeros(1,kn-m_ignore);
snr_vary3=zeros(1,kn-m_ignore);


snr_varyn=zeros(1,kn-m_ignore);
snr_l=zeros(1,kn-m_ignore);

for i_snr = 1:1:length(snr_db)
   snr_fact = sqrt(10^(snr_db(i_snr)/10));
   Rs = snr_fact*v*v'*snr_fact';
   Ru = inr_fact*u*u'*inr_fact'+eye(N);
   R = Rs+Ru;
   %for vary_temp = 1:1:kn
     for j_monte=1:1:m_carlo
       j_monte
       signal_power = snr_fact.*exp(j*2*pi*rand);
       gan_power = inr_fact.*exp(j*2*pi*rand);
       for i_temp = 1:1:kn 
         signal_power_w = signal_power.*((randn+sqrt(-1)*randn)/(sqrt(2)));
         gan_power_w = gan_power.*((randn+sqrt(-1)*randn)/(sqrt(2)));
         n=(randn(N,1)+sqrt(-1)*randn(N,1))/(sqrt(2)); 
         x_temp=signal_power_w*v+gan_power_w*u+n;
         u_temp=gan_power_w*u;
         %Ru_temp=u_temp*(u_temp')+eye(N);
         %Rx_temp=x_temp*(x_temp');
         %Rx=Rx+Rx_temp;
         %Ru=Ru+Ru_temp;
         if i_temp==1
           Rx=x_temp*(x_temp');   
         else
           Rx=((i_temp-1)*Rx+x_temp*x_temp')/i_temp;
         end   
         %______________________________________
         RxL1=Rx+plus(1)*eye(N);
         %Rs=signal_power_w*v*v'*signal_power_w;
           if i_temp>m_ignore
           diedai(i_temp-m_ignore)=i_temp-m_ignore;  
           Bx=v'*inv(RxL1)*v;
           Wx=inv(RxL1)*v*inv(Bx); 
           snr_temp=abs((Wx'*Rs*Wx)/(Wx'*Ru*Wx));       
           if j_monte==1
              snr_vary1(i_snr,i_temp-m_ignore)=10*log10(snr_temp);
           else
              snr_vary1(i_snr,i_temp-m_ignore)=((j_monte-1)*snr_vary1(i_snr,i_temp-m_ignore)+10*log10(snr_temp))/j_monte;
           end 
         end
         %______________________________________
         RxL2=Rx+plus(2)*eye(N);
         %Rs=signal_power_w*v*v'*signal_power_w;
           if i_temp>m_ignore
           diedai(i_temp-m_ignore)=i_temp-m_ignore;  
           Bx=v'*inv(RxL2)*v;
           Wx=inv(RxL2)*v*inv(Bx); 
           snr_temp=abs((Wx'*Rs*Wx)/(Wx'*Ru*Wx));       
           if j_monte==1
              snr_vary2(i_snr,i_temp-m_ignore)=10*log10(snr_temp);
           else
              snr_vary2(i_snr,i_temp-m_ignore)=((j_monte-1)*snr_vary2(i_snr,i_temp-m_ignore)+10*log10(snr_temp))/j_monte;
           end 
         end
         %______________________________________
         RxL3=Rx+plus(3)*eye(N);
         %Rs=signal_power_w*v*v'*signal_power_w;
           if i_temp>m_ignore
           diedai(i_temp-m_ignore)=i_temp-m_ignore;  
           Bx=v'*inv(RxL3)*v;
           Wx=inv(RxL3)*v*inv(Bx); 
           snr_temp=abs((Wx'*Rs*Wx)/(Wx'*Ru*Wx));       
           if j_monte==1
              snr_vary3(i_snr,i_temp-m_ignore)=10*log10(snr_temp);
           else
              snr_vary3(i_snr,i_temp-m_ignore)=((j_monte-1)*snr_vary3(i_snr,i_temp-m_ignore)+10*log10(snr_temp))/j_monte;
           end 
         end
        %-----------------------------------------
           if i_temp>m_ignore
           Bx=v'*inv(Rx)*v;
           Wx=inv(Rx)*v*inv(Bx);
           snr_temp=abs((Wx'*Rs*Wx)/(Wx'*Ru*Wx));
           if j_monte==1
             snr_varyn(i_snr,i_temp-m_ignore)=10*log10(snr_temp);
           else 
             snr_varyn(i_snr,i_temp-m_ignore)=((j_monte-1)*snr_varyn(i_snr,i_temp-m_ignore)+10*log10(snr_temp))/j_monte;
           end
         end
         %Rx=Rx-plus1*eye(N);
         
         if i_temp>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,i_temp-m_ignore) = 10*log10(snr_l_temp);
         end  
      end 
     end
end

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

% semilogx(10.^wy,snr_varyn(1,:),'k+');
% hold on;
semilogx(10.^wy,snr_vary3(1,:),'r=');
hold on;
semilogx(10.^wy,snr_vary3(2,:),'g/');
hold on;
semilogx(10.^wy,snr_vary3(3,:),'b^');
hold on;
semilogx(10.^wy,snr_l(1,:),'c*');
hold on;
semilogx(10.^wy,snr_l(2,:),'c+');
hold on;
semilogx(10.^wy,snr_l(3,:),'co');
hold off;

legend('SNR=0db dl ','SNR=10db dl','20db dl','optimum SNR=0db','optimum SNR=10db','doptimum SNR=20db'); 
xlabel('The iterative number');
ylabel('Average SINR(db)');
title('INR=10db,SNR=0,10,20db,LNR=0db');
% 
% p = [50,50,300,300];
% 
% set(1,'color','white');
% set(1,'position',p);
% plot(diedai,snr_varyn(1,:),'k');
% hold on;
% plot(diedai,snr_vary(1,:),'r');
% hold on;
% plot(diedai,snr_l(1,:),'g');
% hold off;
% title('smi mvpr','FontWeight','bold');
% xlabel('Iterative number','FontWeight','bold');
% ylabel('Snr of output','FontWeight','bold');
% legend('smi 10','smi coefficient','best'); 
 grid on;