function [pb,ps]=bper4pskG4(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=4;                                          %4PSK
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 1;
             1 0;
             1 1];
mapping=[cos(pi/4)+sin(pi/4).*i;
     cos(3*pi/4)+sin(3*pi/4).*i;
     cos(7*pi/4)+sin(7*pi/4).*i;
     cos(5*pi/4)+sin(5*pi/4).*i];               %4PSK
 for k=1:N,                                      %mapping symbol and bits
     fpsk_sig(k)=mapping(dsource(k));
     dsourcebits(k,:)=mappingbits(dsource(k),:);
 end;
 fpsk_sigkrate=[fpsk_sig(1:(N/4));fpsk_sig((N/4+1):2.*N/4);fpsk_sig((2.*N/4+1):3.*N/4);fpsk_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  receive antenna1
     [a(3) a(4)]=gngauss(sgmaa);                %path gain antenna2  receive antenna1
     [a(5) a(6)]=gngauss(sgmaa);                %path gain antenna3  receive antenna1
     [a(7) a(8)]=gngauss(sgmaa);                %path gain antenna4  receive antenna1
     pathgain1=[(a(1)+a(2).*i);
               (a(3)+a(4).*i);
               (a(5)+a(6).*i);
               (a(7)+a(8).*i)]; 
     [a(9) a(10)]=gngauss(sgmaa);               %path gain antenna1  receive antenna2
     [a(11) a(12)]=gngauss(sgmaa);              %path gain antenna2  receive antenna2
     [a(13) a(14)]=gngauss(sgmaa);              %path gain antenna3  receive antenna2
     [a(15) a(16)]=gngauss(sgmaa);              %path gain antenna4  receive antenna2
     pathgain2=[(a(9)+a(10).*i);                %path gain receive antenna2
               (a(11)+a(12).*i);
               (a(13)+a(14).*i);
               (a(15)+a(16).*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
     noise1=[(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)];
     [n(17) n(18)]=gngauss(sgman);              %noise timeslot1  receive antenna2
     [n(19) n(20)]=gngauss(sgman);              %noise timeslot2  receive antenna2
     [n(21) n(22)]=gngauss(sgman);              %noise timeslot3  receive antenna2
     [n(23) n(24)]=gngauss(sgman);              %noise timeslot4  receive antenna2
     [n(25) n(26)]=gngauss(sgman);              %noise timeslot5  receive antenna2
     [n(27) n(28)]=gngauss(sgman);              %noise timeslot6  receive antenna2
     [n(29) n(30)]=gngauss(sgman);              %noise timeslot7  receive antenna2
     [n(31) n(32)]=gngauss(sgman);              %noise timeslot8  receive antenna2
     noise2=[(n(17)+n(18).*i);                  %noise receive antenna2
            (n(19)+n(20).*i);
            (n(21)+n(22).*i);
            (n(23)+n(24).*i);
            (n(25)+n(26).*i);
            (n(27)+n(28).*i);
            (n(29)+n(30).*i);
            (n(31)+n(32).*i)]; 
     C=[fpsk_sigkrate(1,k),fpsk_sigkrate(2,k),fpsk_sigkrate(3,k),fpsk_sigkrate(4,k);
        -fpsk_sigkrate(2,k),fpsk_sigkrate(1,k),-fpsk_sigkrate(4,k),fpsk_sigkrate(3,k);
        -fpsk_sigkrate(3,k),fpsk_sigkrate(4,k),fpsk_sigkrate(1,k),-fpsk_sigkrate(2,k);
        -fpsk_sigkrate(4,k),-fpsk_sigkrate(3,k),fpsk_sigkrate(2,k),fpsk_sigkrate(1,k);
        conj(fpsk_sigkrate(1,k)),conj(fpsk_sigkrate(2,k)),conj(fpsk_sigkrate(3,k)),conj(fpsk_sigkrate(4,k));
        -conj(fpsk_sigkrate(2,k)),conj(fpsk_sigkrate(1,k)),-conj(fpsk_sigkrate(4,k)),conj(fpsk_sigkrate(3,k));
        -conj(fpsk_sigkrate(3,k)),conj(fpsk_sigkrate(4,k)),conj(fpsk_sigkrate(1,k)),-conj(fpsk_sigkrate(2,k));
        -conj(fpsk_sigkrate(4,k)),-conj(fpsk_sigkrate(3,k)),conj(fpsk_sigkrate(2,k)),conj(fpsk_sigkrate(1,k))];
     r1(:,k)=C*pathgain1+noise1;                  %received signals  receive antenna1
     r2(:,k)=C*pathgain2+noise2;                  %received signals  receive antenna2
     for j=1:M,
         s1(k,j)=(abs((r1(1,k)*conj(pathgain1(1,1))+r1(2,k)*conj(pathgain1(2,1))+r1(3,k)*conj(pathgain1(3,1))+r1(4,k)*conj(pathgain1(4,1))+...
             conj(r1(5,k))*pathgain1(1,1)+conj(r1(6,k))*pathgain1(2,1)+conj(r1(7,k))*pathgain1(3,1)+conj(r1(8,k))*pathgain1(4,1)+...
             r2(1,k)*conj(pathgain2(1,1))+r2(2,k)*conj(pathgain2(2,1))+r2(3,k)*conj(pathgain2(3,1))+r2(4,k)*conj(pathgain2(4,1))+...
             conj(r2(5,k))*pathgain2(1,1)+conj(r2(6,k))*pathgain2(2,1)+conj(r2(7,k))*pathgain2(3,1)+conj(r2(8,k))*pathgain2(4,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain1(1,1))).^2+(abs(pathgain1(2,1))).^2+(abs(pathgain1(3,1))).^2+(abs(pathgain1(4,1))).^2+...
             (abs(pathgain2(1,1))).^2+(abs(pathgain2(2,1))).^2+(abs(pathgain2(3,1))).^2+(abs(pathgain2(4,1))).^2)).*(abs(mapping(j))).^2;
         s2(k,j)=(abs((r1(1,k)*conj(pathgain1(2,1))-r1(2,k)*conj(pathgain1(1,1))-r1(3,k)*conj(pathgain1(4,1))+r1(4,k)*conj(pathgain1(3,1))+...
             conj(r1(5,k))*pathgain1(2,1)-conj(r1(6,k))*pathgain1(1,1)-conj(r1(7,k))*pathgain1(4,1)+conj(r1(8,k))*pathgain1(3,1)+...
             r2(1,k)*conj(pathgain2(2,1))-r2(2,k)*conj(pathgain2(1,1))-r2(3,k)*conj(pathgain2(4,1))+r2(4,k)*conj(pathgain2(3,1))+...
             conj(r2(5,k))*pathgain2(2,1)-conj(r2(6,k))*pathgain2(1,1)-conj(r2(7,k))*pathgain2(4,1)+conj(r2(8,k))*pathgain2(3,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain1(1,1))).^2+(abs(pathgain1(2,1))).^2+(abs(pathgain1(3,1))).^2+(abs(pathgain1(4,1))).^2+...
             (abs(pathgain2(1,1))).^2+(abs(pathgain2(2,1))).^2+(abs(pathgain2(3,1))).^2+(abs(pathgain2(4,1))).^2)).*(abs(mapping(j))).^2;
         s3(k,j)=(abs((r1(1,k)*conj(pathgain1(3,1))+r1(2,k)*conj(pathgain1(4,1))-r1(3,k)*conj(pathgain1(1,1))-r1(4,k)*conj(pathgain1(2,1))+...
             conj(r1(5,k))*pathgain1(3,1)+conj(r1(6,k))*pathgain1(4,1)-conj(r1(7,k))*pathgain1(1,1)-conj(r1(8,k))*pathgain1(2,1)+...
             r2(1,k)*conj(pathgain2(3,1))+r2(2,k)*conj(pathgain2(4,1))-r2(3,k)*conj(pathgain2(1,1))-r2(4,k)*conj(pathgain2(2,1))+...
             conj(r2(5,k))*pathgain2(3,1)+conj(r2(6,k))*pathgain2(4,1)-conj(r2(7,k))*pathgain2(1,1)-conj(r2(8,k))*pathgain2(2,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain1(1,1))).^2+(abs(pathgain1(2,1))).^2+(abs(pathgain1(3,1))).^2+(abs(pathgain1(4,1))).^2+...
             (abs(pathgain2(1,1))).^2+(abs(pathgain2(2,1))).^2+(abs(pathgain2(3,1))).^2+(abs(pathgain2(4,1))).^2)).*(abs(mapping(j))).^2;
         s4(k,j)=(abs((r1(1,k)*conj(pathgain1(4,1))-r1(2,k)*conj(pathgain1(3,1))+r1(3,k)*conj(pathgain1(2,1))-r1(4,k)*conj(pathgain1(1,1))+...
             conj(r1(5,k))*pathgain1(4,1)-conj(r1(6,k))*pathgain1(3,1)+conj(r1(7,k))*pathgain1(2,1)-conj(r1(8,k))*pathgain1(1,1)+...
             r2(1,k)*conj(pathgain2(4,1))-r2(2,k)*conj(pathgain2(3,1))+r2(3,k)*conj(pathgain2(2,1))-r2(4,k)*conj(pathgain2(1,1))+...
             conj(r2(5,k))*pathgain2(4,1)-conj(r2(6,k))*pathgain2(3,1)+conj(r2(7,k))*pathgain2(2,1)-conj(r2(8,k))*pathgain2(1,1))-mapping(j))).^2+...
             (-1+2.*((abs(pathgain1(1,1))).^2+(abs(pathgain1(2,1))).^2+(abs(pathgain1(3,1))).^2+(abs(pathgain1(4,1))).^2+...
             (abs(pathgain2(1,1))).^2+(abs(pathgain2(2,1))).^2+(abs(pathgain2(3,1))).^2+(abs(pathgain2(4,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);