function [pb,ps]=bper16QAMH3(snr_in_dB,N)
echo on
E=1;                                           %energy per symbol;
snr=10^(snr_in_dB/10);
numan=4;                                       %the number of transmit antennas
M=16;                                          %16QAM
L=log(M)/log(2);
d=sqrt(1/10);
sgman=sqrt(E*numan/(2*snr));                 %noise variance
sgmaa=sqrt(0.5);
for k=1:N,
    temp=rand;
    dsource(k)=1+floor(M*temp);                % a integer number between 1 and M, uniform       
end;
dsource1=[dsource(1:N/3)];                     % rate k=3
dsource2=[dsource((N/3+1):2*N/3)];
dsource3=[dsource((2*N/3+1):N)];
mappingbits=[0 0 0 0;0 0 0 1;0 0 1 0;0 0 1 1;
             0 1 0 0;0 1 0 1;0 1 1 0;0 1 1 1;
             1 0 0 0;1 0 0 1;1 0 1 0;1 0 1 1;
             1 1 0 0;1 1 0 1;1 1 1 0;1 1 1 1];
mapping=d.*[1+i;  1+3.*i;  3+i;  3+3.*i;
         1-i;  1-3.*i;  3-i;  3-3.*i;
         -1+i; -1+3.*i; -3+i; -3+3.*i;
         -1-i; -1-3.*i; -3-i; -3-3.*i];          %16QAM
 for k=1:N,                                      %mapping symbol and bits
     qam_sig(k)=mapping(dsource(k));
     dsourcebits(k,:)=mappingbits(dsource(k),:);
 end;
 qam_sigkrate=[qam_sig(1:(N/3));qam_sig((N/3+1):2.*N/3);qam_sig((2.*N/3+1):N)];
 dsourcebits1=dsourcebits((1:(N/3)),:);
 dsourcebits2=dsourcebits(((N/3+1):2.*N/3),:);
 dsourcebits3=dsourcebits(((2.*N/3+1):N),:);
 for k=1:(N/3),                                 %each frame 4antennas 4time slots 3symbols
     [a(1) a(2)]=gngauss(sgmaa);                %path gain antenna1
     [a(3) a(4)]=gngauss(sgmaa);                %path gain antenna2
     [a(5) a(6)]=gngauss(sgmaa);                %path gain antenna3
     [a(7) a(8)]=gngauss(sgmaa);                %path gain antenna4
     pathgain=[(a(1)+a(2).*i);
               (a(3)+a(4).*i);
               (a(5)+a(6).*i);
               (a(7)+a(8).*i)]; 
     [n(1) n(2)]=gngauss(sgman);                %noise timeslot1
     [n(3) n(4)]=gngauss(sgman);                %noise timeslot2
     [n(5) n(6)]=gngauss(sgman);                %noise timeslot3
     [n(7) n(8)]=gngauss(sgman);                %noise timeslot4
     noise=[(n(1)+n(2).*i);
            (n(3)+n(4).*i);
            (n(5)+n(6).*i);
            (n(7)+n(8).*i)]; 
     C=[qam_sigkrate(1,k),qam_sigkrate(2,k),(qam_sigkrate(3,k)./sqrt(2)),(qam_sigkrate(3,k)./sqrt(2));
        -conj(qam_sigkrate(2,k)),conj(qam_sigkrate(1,k)),(qam_sigkrate(3,k)./sqrt(2)),(-qam_sigkrate(3,k)./sqrt(2));
        (conj(qam_sigkrate(3,k))./sqrt(2)),(conj(qam_sigkrate(3,k))./sqrt(2)),...
                     ((-qam_sigkrate(1,k)-conj(qam_sigkrate(1,k))+qam_sigkrate(2,k)-conj(qam_sigkrate(2,k)))/2),...
                     ((-qam_sigkrate(2,k)-conj(qam_sigkrate(2,k))+qam_sigkrate(1,k)-conj(qam_sigkrate(1,k)))/2);
        (conj(qam_sigkrate(3,k))./sqrt(2)),(-conj(qam_sigkrate(3,k))./sqrt(2)),...
                     ((qam_sigkrate(2,k)+conj(qam_sigkrate(2,k))+qam_sigkrate(1,k)-conj(qam_sigkrate(1,k)))/2),...
                     (-(qam_sigkrate(1,k)+conj(qam_sigkrate(1,k))+qam_sigkrate(2,k)-conj(qam_sigkrate(2,k)))/2)];
     r(:,k)=C*pathgain+noise;                  %received signals 
     for j=1:M,
         s1(k,j)=(abs((r(1,k)*conj(pathgain(1,1))+conj(r(2,k))*pathgain(2,1)+...
             (r(4,k)-r(3,k))*(conj(pathgain(3,1))-conj(pathgain(4,1)))/2-conj(r(3,k)+r(4,k))*(pathgain(3,1)+pathgain(4,1))/2)-mapping(j))).^2+...
             (-1+(abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2+(abs(pathgain(4,1))).^2).*(abs(mapping(j))).^2;
         s2(k,j)=(abs((r(1,k)*conj(pathgain(2,1))-conj(r(2,k))*pathgain(1,1)+...
             (r(4,k)+r(3,k))*(conj(pathgain(3,1))-conj(pathgain(4,1)))/2+conj(-r(3,k)+r(4,k))*(pathgain(3,1)+pathgain(4,1))/2)-mapping(j))).^2+...
             (-1+(abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2+(abs(pathgain(4,1))).^2).*(abs(mapping(j))).^2;
         s3(k,j)=(abs(((r(1,k)+r(2,k))*conj(pathgain(3,1))/(sqrt(2))+(r(1,k)-r(2,k))*conj(pathgain(4,1))/(sqrt(2))+conj(r(3,k))*(pathgain(1,1)+pathgain(2,1))/(sqrt(2))+...
             conj(r(4,k))*(pathgain(1,1)-pathgain(2,1))/(sqrt(2)))-mapping(j))).^2+...
             (-1+(abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2+(abs(pathgain(4,1))).^2).*(abs(mapping(j))).^2;
     end;
 end;
 numofsymbolerror=0;
 numofbiterror=0;
 for k=1:(N/3),
     [min_metric1 decis1]=min(s1(k,:));
     [min_metric2 decis2]=min(s2(k,:));
     [min_metric3 decis3]=min(s3(k,:));
     decisbits1(k,:)=mappingbits(decis1,:);
     decisbits2(k,:)=mappingbits(decis2,:);
     decisbits3(k,:)=mappingbits(decis3,:);
     for l=1:L,
         if (decisbits1(k,l)~=dsourcebits1(k,l)),
             numofbiterror=numofbiterror+1;
         end;
         if (decisbits2(k,l)~=dsourcebits2(k,l)),
             numofbiterror=numofbiterror+1;
         end;
         if (decisbits3(k,l)~=dsourcebits3(k,l)),
             numofbiterror=numofbiterror+1;
         end;
     end;
     if (decis1~=dsource1(k)),
         numofsymbolerror=numofsymbolerror+1;
     end;
     if (decis2~=dsource2(k)),
         numofsymbolerror=numofsymbolerror+1;
     end;
     if (decis3~=dsource3(k)),
         numofsymbolerror=numofsymbolerror+1;
     end;
 end;
 ps=numofsymbolerror/N;
 pb=numofbiterror/(L*N);