close all;
clear
echo off
Tx_n = 4;
Rx_n = 4;
index = 2;
frame_length=1200;
SNRindB=0:2:20;
for i=1:length (SNRindB),
    SNR(i)=10^(SNRindB(i)/10);
    a=rand(1,frame_length);
    for L=1:frame_length
        if a(L)>.5
            a(L)=1;
        else  a(L)=0;
        end
    end

    a1=modulation(a,index);
    a2=reshape(a1,Tx_n,frame_length/index/Tx_n);
    sigma = 1/sqrt(2*Rx_n*SNR(i));
    AWGN_noise = sigma*(randn(Rx_n,frame_length/index/Tx_n)+j*randn(Rx_n,frame_length/index/Tx_n));
    dec0 = zeros(4,frame_length/index/Tx_n);
    dec1 = zeros(frame_length/index , index);
    dec2 = zeros(frame_length/index , index);
    dec3 = zeros(frame_length/index , index);
    dec4 = zeros(frame_length/index , index);
    for col_idx = 1:frame_length/index/Tx_n

        H=(randn(Rx_n,Tx_n)+j*randn(Rx_n,Tx_n))/sqrt(2);
        h=H;

        r=H*a2(:,col_idx)+AWGN_noise(:,col_idx);
        G=pinv(H);
        G1=inv(H'*H+sigma.^2*eye(Tx_n))*H';
        [gk1 p0]=min(sum(abs(G1).^2,2));
        [gk k0]=min(sum(abs(G).^2,2));
        r_ic = r;
         % ML
      
        count=0:1:15;
        count1=zeros(4,256);
        a4=zeros(4,256);
         for  n1=1:16
             for n2=1:16
       
              d=[-1,1];
              c1=de2bi(count(n2),4);
              c1=d(c1+1);
              c2=de2bi(count(n1),4);
              c2=d(c2+1);
              c=c1+j*c2;
              count1(:,(n1-1)*16+n2)=c.';
             a4(:,(n1-1)*16+n2)=r-H*count1(:,(n1-1)*16+n2);
             end;% n2
         end; % n1

             [gk k2]=min(sum(abs(a4).^2,1));
               dec0(:,col_idx)=count1(:,k2);
  
      
        %           % ZF
        det_zf = G*r;
        dec2((col_idx-1)*Tx_n+[1:Tx_n],:) = reshape(demodulation(det_zf.',index),2,Tx_n).';
        %           MMSE
         det_MMSE=G1*r;
         dec3((col_idx-1)*Tx_n+[1:Tx_n],:) = reshape(demodulation(det_MMSE.',index),2,Tx_n).';

        % SIC1
         for m=1:Tx_n
            k1(m)=k0;
            y=G(k1(m),:)*r_ic;
            dec1((col_idx-1)*Tx_n+k1(m),:) = demodulation(y,index);
            b=modulation(dec1((col_idx-1)*Tx_n+k1(m),:),index);
            r_ic = r_ic - b*H(:, k1(m));
            H(:, k1(m))=zeros(Rx_n,1);
            G=pinv(H);
            temp = sum(abs(G).^2,2);
            temp(k1(1:m)) = 1e10;
            [gk k0]=min(temp);
 
         end; %m
         %SIC2
           for n=1:Tx_n
            p1(n)=p0;
            y1=G1(p1(n),:)*r;
            dec4((col_idx-1)*Tx_n+p1(n),:) = demodulation(y1,index);
            b1=modulation(dec4((col_idx-1)*Tx_n+p1(n),:),index);
            r = r - b1*h(:, p1(n));
            h(:, p1(n))=zeros(Rx_n,1);
            G1=inv(h'*h+sigma.^2*eye(Tx_n))*h';
            temp1 = sum(abs(G1).^2,2);
            temp1(p1(1:n)) = 1e10;
            [gk1 p0]=min(temp1);
           end; %m
    end % col_idx
    
    dec0=reshape(dec0,1,frame_length/index);
    dec0=demodulation(dec0,index); 
    NumErr0 =sum(dec0~=a)
    NumErr1 =sum(abs(reshape(dec1.',1,frame_length)~=a))
    NumErr2= sum(abs(reshape(dec2.',1,frame_length)~=a))
    NumErr3= sum(abs(reshape(dec3.',1,frame_length)~=a))
    NumErr4= sum(abs(reshape(dec4.',1,frame_length)~=a))
    Ber0(i) = NumErr0/frame_length;
    Ber1(i) = NumErr1/frame_length;
    Ber2(i) = NumErr2/frame_length;
    Ber3(i) = NumErr3/frame_length;
    Ber4(i) = NumErr4/frame_length;
end; % Eb/N0
semilogy (SNRindB,Ber0,'-x')
hold on
semilogy (SNRindB,Ber1,'-*')
semilogy (SNRindB,Ber2,'-+')
semilogy (SNRindB,Ber3,'-o')
semilogy (SNRindB,Ber4,'-.')


grid on;
xlabel('Eb/No in dB');
ylabel('error probability');