% 本程序在給定信噪比的前提下，求出比特誤碼率和符號誤碼率；
% 本程序的基本過程如下：(1)發送前的信號映射；(2)序列產生；
% (3)在信道中加入衰落和噪聲；(4)在接受端做相關檢測；
% (5)擇大判決;(6)統計并計算誤碼率;

function [pb,ps]=cm_sm32(snr_in_dB)
% [pb,ps]=cm_sm32(snr_in_dB)
% CM_SM3 finds the probability of bit error and
% symbol error for 
% the given value of snr_in_dB, signal to noise ratio in dB.

N=100;
% 產生瑞利分布的因子；
n_factor=0.8;
% energy per symbol
% 符號能量；
E=1;
numofsymbolerror=0;
numofbiterror=0;
counter=0;

% signal to noise ratio
% 信噪比；
snr=10^(snr_in_dB/10); 
% 噪聲方差,由調用函數時輸入的信噪比決定；
sgma=sqrt(E/snr)/2; % noise variance 

% signal mapping
% 信號映射，實際比特流是如下所示的0101的代碼，但再發送的時候
% 要做映射，使他們的平均值為零；
s00=[1 0]; s01=[0 1]; s11=[-1 0]; s10=[0 -1];

% generation of the data source
while(numofbiterror<100)
    % 以下相當于產生100個4進制信號，兩列100個2進制信號；
for i=1:N
    % a uniform random variable between 0 and 1
    temp=rand;  
    % with probability 1/4, source output is "00"
    % 由于是從小到大順序判斷，不會有誤判的可能；
    if (temp<0.25)
        dsource1(i)=0; dsource2(i)=0;
         % with probability 1/4, source output is "01"
    elseif (temp<0.5)
        dsource1(i)=0; dsource2(i)=1;
        % with probability 1/4, source output is "10"
    elseif (temp<0.75)
        dsource1(i)=1; dsource2(i)=0;
        % with probability 1/4, source output is "11"
    else
        dsource1(i)=1; dsource2(i)=1;
    end;
end;
% 下面是檢測部分；
% detection and the probability of error calculation
% r表示接受信號，信號在傳輸的過程中受到信道衰，減加入噪聲；
% 上面的命令表示產生100個發送端符號，下面的命令是逐個發送，相關
% 判決；
for i=1:N
    
    % 這里信道的衰減因子是一定的,
    % 這里把衰落因子拿到程序前面去設定;
    ray=raylrnd(n_factor);
    % 2 normal distributed r.v with 0, variance sgma
    n=sgma*randn(1,2);                          
    if ((dsource1(i)==0) & (dsource2(i)==0))
        r=ray*s00+n;
    elseif ((dsource1(i)==0) & (dsource2(i)==1))
        r=ray*s01+n;
    elseif ((dsource1(i)==1) & (dsource2(i)==0))
        r=s10*ray+n;
    else
        r=s11*ray+n;
    end;
    % The correlation metrics are computed below
    % 接受端作相關檢測；
    c00=dot(r,s00); c01=dot(r,s01);
    c10=dot(r,s10); c11=dot(r,s11);
    % The decision on the ith symbol is made next
    % 判決，選擇相關值最大的；
    c_max=max([c00,c01,c10,c11]);
    if (c00==c_max)
        decis1=0; decis2=0;
    elseif (c01==c_max)
        decis1=0; decis2=1;
    elseif (c10==c_max)
        decis1=1; decis2=0;
    else
        decis1=1; decis2=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 (symbolerror==1)
        numofsymbolerror=numofsymbolerror+1;
    end;
end;
counter=counter+1;
end;

% 基本過程是這樣的: 產生多次循環,每次循環產生100個符號,對這100個
% 符號統計比特誤碼個數,如果誤比特個數小于100,再進行下一次循環;

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