%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%發射天線為5，接收天線為1的1/2速率的STBC編解碼仿真
%建立人：費禮
%建立時間：12月2日
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [FER,FER_uncoded,SER,SER_uncoded, BER, BER_uncoded]=stbc51(channel_model,K,Num,no_tx_antennas,no_rx_antennas,modulation)

%SNR upto 20 dBs
EbNo=[0:2:20];

%N, M: number of transmit and receive antennas

N=no_tx_antennas;
M=no_rx_antennas;

%initialize packet_count
idx=1;

h=waitbar(0,'Percentage Completed');
set(h,'Position',[230 60 275.25 56.25]);
set(h,'name','Please wait...');
wb=100/(length(EbNo));

for SNR=EbNo
%     sigma=0.5/(10^(SNR/10));
       for packet_count=1:Num
% we are interested in transmitting 'K' SYMBOLS not bits. Hence, K*2 for QPSK
% etc.
        switch (modulation)
                case 'BPSK '
                    data=randint(K,8);
                    BIT=1;
                case 'QPSK '
                    data=randint(K*2,8);
                    BIT=2;
                case '8QAM '
                    data=randint(K*3,8);
                    BIT=3;
                case '16QAM'
                    data=randint(K*4,8);
                    BIT=4;
                case '32QAM'
                    data=randint(K*5,8);
                    BIT=5;
                case '64QAM'
                    data=randint(K*6,8);
                    BIT=6;
                otherwise
                    disp('No Modulation')
        end
sigma_un=0.5/((10^(SNR/10))*BIT);%未編碼下的Eb/No與SNR之間的轉換關系
sigma=0.5/((10^(SNR/10))*BIT*0.5);;%Eb/No與SNR之間的轉換關系
tx_bits=data.';
temp1=[];
temp=[];
for i=1:8 %N is actually 5, but we require 8 tx antenna coding to be transmitted from 5 antennas.
    [temp1 s P]=tx_modulate(tx_bits(i,:),modulation);
    temp=[temp; temp1];
    temp1=0;
end

%ready to transmit symbols of length 'K'
X=temp.';
fr_length=length(X);
 
%block signals in the lth time slot-Block coding for G5 STBC 
 x0=X(:,1);% required to verify a 1x1 system
 x1=X(:,1:5);
 x2(:,1)=-X(:,2); x2(:,2)=X(:,1); x2(:,3)=X(:,4);x2(:,4)=-X(:,3);x2(:,5)=X(:,6);
 x3(:,1)=-X(:,3); x3(:,2)=-X(:,4); x3(:,3)=X(:,1);x3(:,4)=X(:,2);x3(:,5)=X(:,7);
 x4(:,1)=-X(:,4); x4(:,2)=X(:,3); x4(:,3)=-X(:,2);x4(:,4)=X(:,1);x4(:,5)=X(:,8);
 x5(:,1)=-X(:,5); x5(:,2)=-X(:,6); x5(:,3)=-X(:,7);x5(:,4)=-X(:,8);x5(:,5)=X(:,1);
 x6(:,1)=-X(:,6); x6(:,2)=X(:,5); x6(:,3)=-X(:,8);x6(:,4)=X(:,7);x6(:,5)=-X(:,2);
 x7(:,1)=-X(:,7); x7(:,2)=X(:,8); x7(:,3)=X(:,5);x7(:,4)=-X(:,6);x7(:,5)=-X(:,3);
 x8(:,1)=-X(:,8); x8(:,2)=-X(:,7); x8(:,3)=X(:,6);x8(:,4)=X(:,5);x8(:,5)=-X(:,4);
 x9=conj(x1);
 x10=conj(x2);
 x11=conj(x3);
 x12=conj(x4);
 x13=conj(x5);
 x14=conj(x6);
 x15=conj(x7);
 x16=conj(x8);

 % form the channel matrix
 for n=1:N %這里也許應該是M
     if channel_model=='AWGN    '
         Hr(n,:,:)=ones(fr_length,N);
     else
         Hr(n,:,:)=(randn(fr_length,N)+j*randn(fr_length,N))/sqrt(2);
     end
 end

 
 for n=1:M

%transmission matrix   
H=reshape(Hr(n,:,:),fr_length,N);
Habs(:,n)=sum(abs(H).^2,2);

 %received signal per receiver antenna   
 r1(:,n)=sum(H.*x1,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r2(:,n)=sum(H.*x2,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r3(:,n)=sum(H.*x3,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r4(:,n)=sum(H.*x4,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r5(:,n)=sum(H.*x5,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r6(:,n)=sum(H.*x6,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r7(:,n)=sum(H.*x7,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r8(:,n)=sum(H.*x8,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r9(:,n)=sum(H.*x9,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r10(:,n)=sum(H.*x10,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r11(:,n)=sum(H.*x11,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r12(:,n)=sum(H.*x12,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r13(:,n)=sum(H.*x13,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r14(:,n)=sum(H.*x14,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 
 r15(:,n)=sum(H.*x15,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));
 r16(:,n)=sum(H.*x16,2)/sqrt(N)+sqrt(sigma)*(randn(fr_length,1)+j*randn(fr_length,1));

 % demodulate the received signals 
  z1_1(:,n)=r1(:,n).*conj(H(:,1))+r2(:,n).*conj(H(:,2))+r3(:,n).*conj(H(:,3))+r4(:,n).*conj(H(:,4))+r5(:,n).*conj(H(:,5));
  z1_2(:,n)=conj(r9(:,n)).*H(:,1)+conj(r10(:,n)).*H(:,2)+conj(r11(:,n)).*H(:,3)+conj(r12(:,n)).*H(:,4)+conj(r13(:,n)).*H(:,5);
  z1(:,n)=z1_1(:,n)+z1_2(:,n);
  
  z2_1(:,n)=r1(:,n).*conj(H(:,2))-r2(:,n).*conj(H(:,1))+r3(:,n).*conj(H(:,4))-r4(:,n).*conj(H(:,3))-r6(:,n).*conj(H(:,5));
  z2_2(:,n)=conj(r9(:,n)).*H(:,2)-conj(r10(:,n)).*H(:,1)+conj(r11(:,n)).*H(:,4)-conj(r12(:,n)).*H(:,3)-conj(r14(:,n)).*H(:,5);
  z2(:,n)=z2_1(:,n)+z2_2(:,n);
  
  z3_1(:,n)=r1(:,n).*conj(H(:,3))-r2(:,n).*conj(H(:,4))-r3(:,n).*conj(H(:,1))+r4(:,n).*conj(H(:,2))-r7(:,n).*conj(H(:,5));
  z3_2(:,n)=conj(r9(:,n)).*H(:,3)-conj(r10(:,n)).*H(:,4)-conj(r11(:,n)).*H(:,1)+conj(r12(:,n)).*H(:,2)-conj(r15(:,n)).*H(:,5);
  z3(:,n)=z3_1(:,n)+z3_2(:,n);
  
  z4_1(:,n)=r1(:,n).*conj(H(:,4))+r2(:,n).*conj(H(:,3))-r3(:,n).*conj(H(:,2))-r4(:,n).*conj(H(:,1))-r8(:,n).*conj(H(:,5));
  z4_2(:,n)=conj(r9(:,n)).*H(:,4)+conj(r10(:,n)).*H(:,3)-conj(r11(:,n)).*H(:,2)-conj(r12(:,n)).*H(:,1)-conj(r16(:,n)).*H(:,5);
  z4(:,n)=z4_1(:,n)+z4_2(:,n);
  
  z5_1(:,n)=r1(:,n).*conj(H(:,5))-r5(:,n).*conj(H(:,1))+r6(:,n).*conj(H(:,2))+r7(:,n).*conj(H(:,3))+r8(:,n).*conj(H(:,4));
  z5_2(:,n)=conj(r9(:,n)).*H(:,5)-conj(r13(:,n)).*H(:,1)+conj(r14(:,n)).*H(:,2)+conj(r15(:,n)).*H(:,3)+conj(r16(:,n)).*H(:,4);
  z5(:,n)=z5_1(:,n)+z5_2(:,n);
  
  z6_1(:,n)=r2(:,n).*conj(H(:,5))-r5(:,n).*conj(H(:,2))-r6(:,n).*conj(H(:,1))-r7(:,n).*conj(H(:,4))+r8(:,n).*conj(H(:,3));
  z6_2(:,n)=conj(r10(:,n)).*H(:,5)-conj(r13(:,n)).*H(:,2)-conj(r14(:,n)).*H(:,1)-conj(r15(:,n)).*H(:,4)+conj(r16(:,n)).*H(:,3);
  z6(:,n)=z6_1(:,n)+z6_2(:,n);
  
  z7_1(:,n)=r3(:,n).*conj(H(:,5))-r5(:,n).*conj(H(:,3))+r6(:,n).*conj(H(:,4))-r7(:,n).*conj(H(:,1))-r8(:,n).*conj(H(:,2));
  z7_2(:,n)=conj(r11(:,n)).*H(:,5)-conj(r13(:,n)).*H(:,3)+conj(r14(:,n)).*H(:,4)-conj(r15(:,n)).*H(:,1)-conj(r16(:,n)).*H(:,2);
  z7(:,n)=z7_1(:,n)+z7_2(:,n);
  
  z8_1(:,n)=r4(:,n).*conj(H(:,5))-r5(:,n).*conj(H(:,4))-r6(:,n).*conj(H(:,3))+r7(:,n).*conj(H(:,2))-r8(:,n).*conj(H(:,1));
  z8_2(:,n)=conj(r12(:,n)).*H(:,5)-conj(r13(:,n)).*H(:,4)-conj(r14(:,n)).*H(:,3)+conj(r15(:,n)).*H(:,2)-conj(r16(:,n)).*H(:,1);
  z8(:,n)=z8_1(:,n)+z8_2(:,n);
 end

%uncoded(1,1)
r01=H(:,1).*x0+sqrt(sigma_un)*(randn(fr_length,1)+j*randn(fr_length,1));
     
%form estimates
for m=1:P
   d01(:,m)=abs(r01-H(:,1)*s(m)).^2; %uncoded signal
%coded signals   
   d1(:,m)=abs(sum(z1,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d2(:,m)=abs(sum(z2,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d3(:,m)=abs(sum(z3,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d4(:,m)=abs(sum(z4,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d5(:,m)=abs(sum(z5,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d6(:,m)=abs(sum(z6,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d7(:,m)=abs(sum(z7,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
   d8(:,m)=abs(sum(z8,2)-s(m)).^2+(-1+2*sum(Habs,2))*abs(s(m))^2;
end
% determine the minimum of estimates      

%decision for detecting uncoded
[y0,i0]=min((d01),[],2);
      s0d=s(i0).';
      
      clear d01 
%decision for detecting s1     
      [y1,i1]=min(d1,[],2);
      s1d=s(i1).';
      
      clear d1
%decision for detecting s2     
      [y2,i2]=min(d2,[],2);
      s2d=s(i2).';
      
      clear d2
%decision for detecting s3     
      [y3,i3]=min(d3,[],2);
      s3d=s(i3).';
      
      clear d3
%decision for detecting s4     
      [y4,i4]=min(d4,[],2);
      s4d=s(i4).';
      
      clear d4
      
%decision for detecting s5     
      [y5,i5]=min(d5,[],2);
      s5d=s(i5).';
      
      clear d5
%decision for detecting s6     
      [y6,i6]=min(d6,[],2);
      s6d=s(i6).';
      
      clear d6
%decision for detecting s7     
      [y7,i7]=min(d7,[],2);
      s7d=s(i7).';
      
      clear d7
%decision for detecting s8     
      [y8,i8]=min(d8,[],2);
      s8d=s(i8).';
      
      clear d8

% form received symbols    
    Xd=[s1d s2d s3d s4d s5d s6d s7d s8d];
      
%determine symbol errors   
error_un(packet_count)=sum(X(:,1)~=s0d);% for uncoded
temp1=X>0;
temp2=Xd>0;
error(packet_count)=sum(sum(temp1~=temp2));% for coded
       end% end of FOR loop for "packet_count"
%calculate FER, SER and BER for current idx

%for uncoded signal
  SER_uncoded(idx)=sum(error_un)/(Num*K);
  %采用格雷碼映射時應除以BIT/2 但BPSK除外
  if BIT == 1
  BER_uncoded(idx)=SER_uncoded(idx)/BIT;
  else
  BER_uncoded(idx)=SER_uncoded(idx)/(BIT/2);
  end
  FER_uncoded(idx)=SER_uncoded(idx)*K;

  %for coded signal
  SER(idx)=sum(error)/(Num*K);
  %采用格雷碼映射時應除以BIT/2 但BPSK除外
  if BIT == 1
  BER(idx)=SER(idx)/BIT;
  else
  BER(idx)=SER(idx)/(BIT/2);
  end
  FER(idx)=SER(idx)*K;
  
%increment idx  
  idx=idx + 1;
  
  
str_bar=[num2str(wb) '% Completed'];
waitbar(wb/100,h,str_bar);
wb=wb+100/(length(EbNo));
  
end% end of FOR loop for SNR
close(h);