function [pb,ps]=bper16QAMG3(snr_in_dB,N)
echo on
E=1;                                          %energy per symbol;
snr=10^(snr_in_dB/10);
numan=3;                                      %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/4)];                    % rate k=4
dsource2=[dsource((N/4+1):2*N/4)];
dsource3=[dsource((2*N/4+1):3*N/4)];
dsource4=[dsource((3*N/4+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/4));qam_sig((N/4+1):2.*N/4);qam_sig((2.*N/4+1):3.*N/4);qam_sig((3.*N/4+1):N)];
 dsourcebits1=dsourcebits((1:(N/4)),:);
 dsourcebits2=dsourcebits(((N/4+1):2.*N/4),:);
 dsourcebits3=dsourcebits(((2.*N/4+1):3.*N/4),:);
 dsourcebits4=dsourcebits(((3.*N/4+1):N),:);
 for k=1:(N/4),                                 %each frame 3antennas 8time slots 4symbols
     [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
     pathgain=[(a(1)+a(2).*i);
               (a(3)+a(4).*i);
               (a(5)+a(6).*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
     [n(9) n(10)]=gngauss(sgman);                %noise timeslot5
     [n(11) n(12)]=gngauss(sgman);                %noise timeslot6
     [n(13) n(14)]=gngauss(sgman);                %noise timeslot7
     [n(15) n(16)]=gngauss(sgman);                %noise timeslot8
     noise=[(n(1)+n(2).*i);
            (n(3)+n(4).*i);
            (n(5)+n(6).*i);
            (n(7)+n(8).*i);
            (n(9)+n(10).*i);
            (n(11)+n(12).*i);
            (n(13)+n(14).*i);
            (n(15)+n(16).*i)]; 
     C=[qam_sigkrate(1,k),qam_sigkrate(2,k),qam_sigkrate(3,k);
        -qam_sigkrate(2,k),qam_sigkrate(1,k),-qam_sigkrate(4,k);
        -qam_sigkrate(3,k),qam_sigkrate(4,k),qam_sigkrate(1,k);
        -qam_sigkrate(4,k),-qam_sigkrate(3,k),qam_sigkrate(2,k);
        conj(qam_sigkrate(1,k)),conj(qam_sigkrate(2,k)),conj(qam_sigkrate(3,k));
        -conj(qam_sigkrate(2,k)),conj(qam_sigkrate(1,k)),-conj(qam_sigkrate(4,k));
        -conj(qam_sigkrate(3,k)),conj(qam_sigkrate(4,k)),conj(qam_sigkrate(1,k));
        -conj(qam_sigkrate(4,k)),-conj(qam_sigkrate(3,k)),conj(qam_sigkrate(2,k))];
     r(:,k)=C*pathgain+noise;                  %received signals 
     for j=1:M,
         s1(k,j)=(abs((r(1,k)*conj(pathgain(1,1))+r(2,k)*conj(pathgain(2,1))+r(3,k)*conj(pathgain(3,1))+...
             conj(r(5,k))*pathgain(1,1)+conj(r(6,k))*pathgain(2,1)+conj(r(7,k))*pathgain(3,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2)).*(abs(mapping(j))).^2;
         s2(k,j)=(abs((r(1,k)*conj(pathgain(2,1))-r(2,k)*conj(pathgain(1,1))+r(4,k)*conj(pathgain(3,1))+...
             conj(r(5,k))*pathgain(2,1)-conj(r(6,k))*pathgain(1,1)+conj(r(8,k))*pathgain(3,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2)).*(abs(mapping(j))).^2;
         s3(k,j)=(abs((r(1,k)*conj(pathgain(3,1))-r(3,k)*conj(pathgain(1,1))-r(4,k)*conj(pathgain(2,1))+...
             conj(r(5,k))*pathgain(3,1)-conj(r(7,k))*pathgain(1,1)-conj(r(8,k))*pathgain(2,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2)).*(abs(mapping(j))).^2;
         s4(k,j)=(abs((-r(2,k)*conj(pathgain(3,1))+r(3,k)*conj(pathgain(2,1))-r(4,k)*conj(pathgain(1,1))-...
             conj(r(6,k))*pathgain(3,1)+conj(r(7,k))*pathgain(2,1)-conj(r(8,k))*pathgain(1,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain(1,1))).^2+(abs(pathgain(2,1))).^2+(abs(pathgain(3,1))).^2)).*(abs(mapping(j))).^2;
     end;
 end;
 numofsymbolerror=0;
 numofbiterror=0;
 for k=1:(N/4),
     [min_metric1 decis1]=min(s1(k,:));
     [min_metric2 decis2]=min(s2(k,:));
     [min_metric3 decis3]=min(s3(k,:));
     [min_metric4 decis4]=min(s4(k,:));
     decisbits1(k,:)=mappingbits(decis1,:);
     decisbits2(k,:)=mappingbits(decis2,:);
     decisbits3(k,:)=mappingbits(decis3,:);
     decisbits4(k,:)=mappingbits(decis4,:);
     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;
         if (decisbits4(k,l)~=dsourcebits4(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;
     if (decis4~=dsource4(k)),
         numofsymbolerror=numofsymbolerror+1;
     end;
 end;
 ps=numofsymbolerror/N;
 pb=numofbiterror/(L*N);