%The BER and throughput comparison among chase combined scheme ,our
%proposed HARQ scheme and STBC scheme

 SNR=-5:5:25;
 Nt=2;
 Nr=2;
 M=4;
SymTime=10;
packetNum=1000;
ERROR_ARQ=[];
ERROR_HARQ=[];
ERROR_STBC=[];
RATE_ARQ=[];
RATE_HARQ=[];
RATE_STBC=[];
%The symbol constellation
for m=1:M
    s(m)=qammod(m-1,M);
end
 d1=[];
d2=[];
  
for SNRindex=1:length(SNR)

       error_arq=0; 
       error_harq=0;
       rate_arq=0;
       rate_harq=0;
       error_stbc=0;
        rate_stbc=0;
    for packetindex=1:packetNum
    H=(randn(Nr,Nt)+j*randn(Nr,Nt))/sqrt(2);
      Out=[];
      data=[];
 for t=1:SymTime 
     data=floor(rand(Nt,1)*M);
     a=qammod(data,M);
     for condition=1:2
         switch condition
             case 1
    r1=awgn(H*a,SNR(SNRindex),'measured');
    Out=blastdetect(H,Nt,r1,M); 
   [number ratio]=biterr(data,Out);
if ratio==0
     rate_harq=rate_harq+2;
      rate_arq=rate_arq+2;
else

   for cond=1:2
       switch cond
       case 1        %repeat the same symbol vector

             r1_1=awgn(H*a,SNR(SNRindex),'measured');
             combined_r1=(r1+r1_1)/2;
             Out=blastdetect(H,Nt,combined_r1,M); 
             [number ratio]=biterr(data,Out);
              if ratio==0
                   rate_arq=rate_arq+1;
              end
         error_arq=error_arq+ratio;  
       case 2       %repeat the stbc code vector
            x=[];
             x(1,1)=-conj(a(2));
             x(2,1)=conj(a(1));
             r2=awgn(H*x,SNR(SNRindex),'measured');
             r=[r1 r2];
           Out=decode22(H,r,M);
       [number ratio]=biterr(data,Out);
        if ratio==0
        rate_harq=rate_harq+1;
        end
        error_harq=error_harq+ratio;
       end
   end
end

       case 2
           code=encode22(a);
           r=awgn(H*code,SNR(SNRindex),'measured');
          Out=decode22(H,r,M);
       [number ratio]=biterr(data,Out);
        if ratio==0
        rate_stbc=rate_stbc+1;
        end
        error_stbc=error_stbc+ratio;
           end
    end
 end
end

ERROR_ARQ(SNRindex)= error_arq/SymTime/packetNum
ERROR_HARQ(SNRindex)= error_harq/SymTime/packetNum
ERROR_STBC(SNRindex)= error_stbc/SymTime/packetNum
RATE_ARQ(SNRindex)=rate_arq./SymTime/packetNum
RATE_HARQ(SNRindex)=rate_harq./SymTime/packetNum
RATE_STBC(SNRindex)=rate_stbc./SymTime/packetNum
end
subplot(1,2,1)
semilogy(SNR,ERROR_ARQ,'-bo',SNR,ERROR_HARQ,'-rs',SNR,ERROR_STBC,'-g+');
title('average BER')
axis([-5 25 10^-5 10^-0])
xlabel('SNR')
ylabel('BER')
grid on
legend('uncoded with "osic" detection','proposed scheme','STBC')

subplot(1,2,2)
plot(SNR,RATE_ARQ,'-bo',SNR,RATE_HARQ,'-rs',SNR,RATE_STBC,'-g+');
title('average transmission rate')
axis([-5 25 0.1 2])
xlabel('SNR')
ylabel('average rate')
grid on
legend('uncoded with "osic" detection','proposed scheme','STBC')