%The delay and throughput comparison between chase combinded scheme and
%our proposed HARQ scheme under the limited retransmission number condition
SNR=-10:5:20;
 Nt=2;
 Nr=2 ;
 M=4;
symtime=10;
packetNum=100;
pureharq_averdelay=[];
stbc22harq_averdelay=[];
pureharq_averber=[];
stbc22harq_averber=[];
pureharq_rate=[];
stbc22harq_rate=[];

for SNRindex=1:length(SNR)
    pureharq_delay=0;stbc22harq_delay=0;
    ber1=0; ber2=0; pureharq_rate=0;stbc22harq_rate=0;
    for packetindex=1:packetNum
           H=(randn(Nr,Nt)+j*randn(Nr,Nt))/sqrt(2);
        for t=1:symtime
          
            data=floor(rand(Nt,1)*M);
           a=qammod(data,M);
             transTime=1;transTime1=1;
          
                        r1=[];
                    Out=[];x=[];
                    x=a;
                   r1=awgn(H*x,SNR(SNRindex),'measured');
                  record_r1=r1;
            
                  Out=blastdetect(H,Nt,r1,M);

                     [number ratio1]=symerr(data,Out);
                     for cond=1:2
                         switch cond
                         case 1
                             ratio=ratio1;
                     r1(:,1)=record_r1;
                     while (ratio~=0&transTime<=30)
                     transTime=transTime+1;
                      r1(:,transTime)=awgn(H*x,SNR(SNRindex),'measured');
                    combined_r1=sum(r1.').'/transTime;
                   Out=blastdetect(H,Nt,combined_r1,M);
                    [number ratio]=biterr(data,Out);
                     end
                    pureharq_delay=pureharq_delay+transTime-1;
                    ber1=ber1+ratio;
              if ratio==0
              pureharq_rate=pureharq_rate+2/transTime;
              end
                        case 2
                            ratio=ratio1;
                            r2=[];code=encode22(a);
                            r2(:,1)=record_r1;
                            Out1=[];
                            while (ratio~=0&transTime1<=10)
                                transTime1=transTime1+1;
                                N1=rem(transTime1-1,2)+1;
                                 if transTime1<=2
                                x=code(:,transTime1);
                                r2(:,transTime1)=awgn(H*x,SNR(SNRindex),'measured');
                             Out=decode22(H,r2,M);
                               elseif transTime1>2   
                       x=code(:,N1);
                      r21=awgn(H*x,SNR(SNRindex),'measured');
                     r2(:,N1)=(r2(:,N1)+r21)/2;
                  Out=decode22(H,r2,M);
                  end
              [number ratio]=biterr(data,Out);
                            end
              stbc22harq_delay=stbc22harq_delay+transTime1-1;
                ber2=ber2+ratio;
          if ratio==0
         stbc22harq_rate=stbc22harq_rate+2/transTime1;
         end
                     end
                end
                            
                   
            end
        end

   pureharq_averdelay(SNRindex)=pureharq_delay/symtime/packetNum
    stbc22harq_averdelay(SNRindex)=stbc22harq_delay/symtime/packetNum
    pureharq_averber(SNRindex)=ber1/symtime/packetNum
    stbc22harq_averber(SNRindex)=ber2/symtime/packetNum
   pureharq_averrate(SNRindex)=pureharq_rate/symtime/packetNum
    stbc22harq_averrate(SNRindex)=stbc22harq_rate/symtime/packetNum
end
   subplot(1,3,1)
   plot(SNR, pureharq_averdelay,'-bo',SNR,stbc22harq_averdelay,'-rs')
   grid on
   title('average delay comparision')
   axis([-10 20 0 10])
   legend('pureharq','orthogonal stbc ')
   subplot(1,3,2)
   semilogy(SNR, pureharq_averber,'-bo',SNR,stbc22harq_averber,'-rs')
   grid on
   title('average BER comparision')
   axis([-10 20 10^-5 10^0])
   legend('pureharq','orthogonal stbc ')
  subplot(1,3,3)
   plot(SNR, pureharq_averrate,'-bo',SNR,stbc22harq_averrate,'-rs')
   grid on
   title('average rate comparision')
   axis([-10 20 0 2])
   legend('pureharq','orthogonal stbc ')