function [FER,FER_uncoded,SER,SER_uncoded, BER, BER_uncoded]=stbc31(channel_model,K,Num,no_tx_antennas,no_rx_antennas,modulation)%SNR upto 20 dBsEbNo=[0:2:20];%N, M: number of transmit and receive antennasN=no_tx_antennas;M=no_rx_antennas;%initialize packet_countidx=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,4);                    BIT=1;                case 'QPSK '                    data=randint(K*2,4);                    BIT=2;                case '8QAM '                    data=randint(K*3,4);                    BIT=3;                case '16QAM'                    data=randint(K*4,4);                    BIT=4;                case '32QAM'                    data=randint(K*5,4);                    BIT=5;                case '64QAM'                    data=randint(K*6,4);                    BIT=6;                otherwise                    disp('No Modulation')        endsigma_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:4 %N is actually 3, but we require 4 tx antenna coding to be transmitted from 3 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 G3 STBC  x0=X(:,1);% required to verify a 1x1 system x1=X(:,1:3); x2(:,1)=-X(:,2); x2(:,2)=X(:,1); x2(:,3)=-X(:,4); x3(:,1)=-X(:,3); x3(:,2)=X(:,4); x3(:,3)=X(:,1); x4(:,1)=-X(:,4); x4(:,2)=-X(:,3); x4(:,3)=X(:,2); x5=conj(x1); x6=conj(x2); x7=conj(x3); x8=conj(x4);  % 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));  % demodulate the received signals  z1_1(:,n)=r1(:,n).*conj(H(:,1))+r2(:,n).*conj(H(:,2))+r3(:,n).*conj(H(:,3));  z1_2(:,n)=conj(r5(:,n)).*H(:,1)+conj(r6(:,n)).*H(:,2)+conj(r7(:,n)).*H(:,3);  z1(:,n)=z1_1(:,n)+z1_2(:,n);    z2_1(:,n)=r1(:,n).*conj(H(:,2))-r2(:,n).*conj(H(:,1))+r4(:,n).*conj(H(:,3));  z2_2(:,n)=conj(r5(:,n)).*H(:,2)-conj(r6(:,n)).*H(:,1)+conj(r8(:,n)).*H(:,3);  z2(:,n)=z2_1(:,n)+z2_2(:,n);    z3_1(:,n)=r1(:,n).*conj(H(:,3))-r3(:,n).*conj(H(:,1))-r4(:,n).*conj(H(:,2));  z3_2(:,n)=conj(r5(:,n)).*H(:,3)-conj(r7(:,n)).*H(:,1)-conj(r8(:,n)).*H(:,2);  z3(:,n)=z3_1(:,n)+z3_2(:,n);    z4_1(:,n)=-r2(:,n).*conj(H(:,3))+r3(:,n).*conj(H(:,2))-r4(:,n).*conj(H(:,1));  z4_2(:,n)=-conj(r6(:,n)).*H(:,3)+conj(r7(:,n)).*H(:,2)-conj(r8(:,n)).*H(:,1);  z4(:,n)=z4_1(:,n)+z4_2(:,n);end%uncoded(1,1)r01=H(:,1).*x0+sqrt(sigma_un)*(randn(fr_length,1)+j*randn(fr_length,1));  %form estimatesfor m=1:P   d01(:,m)=abs(r01-H(:,1)*s(m)).^2; %uncoded signal%coded signals 以下的解碼式似乎缺少系數2         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;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% form received symbols        Xd=[s1d s2d s3d s4d];%determine symbol errors   error_un(packet_count)=sum(X(:,1)~=s0d);% for uncodedtemp1=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);