%***********************************************************************
% cm_sm32.m
% cm_sm32  finds the probability of bit error and symbol error for the 
% given value of snr_in_dB, signal-to-noise ratio in dB.
% Programmed by Liu Dejian
% Dtat:2008-3-8
%***********************************************************************
function[pb,ps]=cm_sm32(snr_in_dB)
%   [pb,ps]=cm_sm32(snr_in_dB)
%		CM_SM32  finds the probability of bit error and symbol error for the 
%   		given value of snr_in_dB, signal-to-noise ratio in dB.

N=10000;                  % 在不同信噪比Eb/N0下發送10000符號進行蒙特卡羅仿真
E=1;                      % energy per symbol（這里E=Es=2Eb）
snr=10^(snr_in_dB/10);    % signal-to-noise ratio
sgma=sqrt(E/(3*2*snr));	  % noise variance
% the signal mapping
s000=[1 0];
s001=[cos(pi/4) sin(pi/4)];
s011=[0 1];
s010=[cos(3*pi/4) sin(3*pi/4)];
s110=[-1 0];
s111=[cos(5*pi/4) sin(5*pi/4)];
s101=[0 -1];
s100=[cos(7*pi/4) sin(7*pi/4)];
% generation of the data source
for i=1:N,
   temp=rand;			  	% a uniform random variable between 0 and 1  
   if(temp<0.125),		  	% With probability 1/8, source output is "000."
      dsource1(i)=0;
      dsource2(i)=0;
      dsource3(i)=0;
   elseif(temp<0.25),		% With probability 1/8, source output is "001."
      dsource1(i)=0;
      dsource2(i)=0;
      dsource3(i)=1;
   elseif(temp<0.375),		% With probability 1/8, source output is "010."
      dsource1(i)=0;
      dsource2(i)=1;
      dsource3(i)=0;
   elseif(temp<0.5),		% With probability 1/8, source output is "011."
      dsource1(i)=0;
      dsource2(i)=1;
      dsource3(i)=1;
   elseif(temp<0.625),		% With probability 1/8, source output is "100."      
      dsource1(i)=1;
      dsource2(i)=0;
      dsource3(i)=0;
   elseif(temp<0.75),		% With probability 1/8, source output is "101."
      dsource1(i)=1;
      dsource2(i)=0;
      dsource3(i)=1;
   elseif(temp<0.875),		% With probability 1/8, source output is "110."
      dsource1(i)=1;
      dsource2(i)=1;
      dsource3(i)=0;
   else             		% With probability 1/8, source output is "111."
      dsource1(i)=1;
      dsource2(i)=1;
      dsource3(i)=1;
   end;
end;
% detection and the probability of error calculation（判決誤碼率的計算）
numofsymbolerror=0;
numofbiterror=0;
for i=1:N,
   % The received signal at the detector, for the ith symbol, is:
   % （在判決器處的接收信號，對于第i個符號為：）
   n(1)=gngauss(sgma);
   n(2)=gngauss(sgma);
   
   if((dsource1(i)==0)&(dsource2(i)==0)&(dsource3(i)==0)),
      r=s000+n;
   elseif((dsource1(i)==0)&(dsource2(i)==0)&(dsource3(i)==1)),
      r=s001+n;
   elseif((dsource1(i)==0)&(dsource2(i)==1)&(dsource3(i)==0)),
      r=s010+n;
   elseif((dsource1(i)==0)&(dsource2(i)==1)&(dsource3(i)==1)),
      r=s011+n;
   elseif((dsource1(i)==1)&(dsource2(i)==0)&(dsource3(i)==0)),
      r=s100+n;
   elseif((dsource1(i)==1)&(dsource2(i)==0)&(dsource3(i)==1)),
      r=s101+n;
   elseif((dsource1(i)==1)&(dsource2(i)==1)&(dsource3(i)==0)),
      r=s110+n;
   else
      r=s111+n;
   end;
   % The correlation metrics are computed below.
   % （以下為計算互相關量度）
   c000=dot(r,s000);
   c001=dot(r,s001);
   c010=dot(r,s010);
   c011=dot(r,s011);
   c100=dot(r,s100);
   c101=dot(r,s101);
   c110=dot(r,s110);
   c111=dot(r,s111);
   % The decision on the ith symbol is made next.
   % （第i個符號的判決按如下方式進行）
   c_max=max([c000 c001 c010 c011 c100 c101 c110 c111]);
   
   if(c000==c_max),
      decis1=0;decis2=0;decis3=0;
   elseif(c001==c_max),
      decis1=0;decis2=0;decis3=1;
   elseif(c010==c_max),
      decis1=0;decis2=1;decis3=0;
   elseif(c011==c_max),
      decis1=0;decis2=1;decis3=1;
   elseif(c100==c_max),
      decis1=1;decis2=0;decis3=0;
   elseif(c101==c_max),
      decis1=1;decis2=0;decis3=1;
   elseif(c110==c_max),
      decis1=1;decis2=1;decis3=0;
   else
      decis1=1;decis2=1;decis3=1;
   end;
   % Increment the error counter, if the decision is not correct.   
   symbolerror=0;
   if(decis1~=dsource1(i)),
      numofbiterror=numofbiterror+1;
      symbolerror=1;
   end;
   
   if(decis2~=dsource2(i)),
      numofbiterror=numofbiterror+1;
      symbolerror=1;
   end;
   
   if(decis3~=dsource3(i)),
      numofbiterror=numofbiterror+1;
      symbolerror=1;
   end;
   
   if(symbolerror==1),
      numofsymbolerror=numofsymbolerror+1;
   end;
end;

ps=numofsymbolerror/N;  % since there are totally N symbols
pb=numofbiterror/(3*N); % since 3N bits are transmitted