

%=======================================================================================
%    Biorthogonal simulation
%
%    TKK - Comlab 2007 - SL
%=======================================================================================
%
% Purpose: simulate a communication system using biorthogonal
%          waveforms and calculate bit and symbol error probabilities.
%
% Channel: AWGN
%
%=======================================================================================

%--------------------------------------
% Settings
%--------------------------------------


clear all

%M = 4;


  % ** STEP 1 **------------------------- 
  % Set M = 2 
  %
  % -------------------------------------  


%M = ...; 

%>> SOLUTION
M = 2; 
%<<


b = log2(M);                                      % 2 bits per symbol

N = 1000;                                   % number of symbols to transmit 

E = 1;                                      % energy


SNR_bit_dB_MAX = 16;                        % maximum SNR in dB to simulate

SNR_bit_dB = 0:1:SNR_bit_dB_MAX;            % vector with the bit SNR in dB = Eb/No


%----------------------------------------
% Repeat the simulation for different SNR
%----------------------------------------



for iSNR = 1:length(SNR_bit_dB)



  SNR_bit_abs = 10^(SNR_bit_dB(iSNR)/10);   % transform SNR in dB to abs values

%  sigma = sqrt(E/(4*SNR_bit_abs));       % standard deviation 


  % ** STEP 2 **------------------------- 
  % Set standard deviation
  %
  % -------------------------------------  

  %  sigma = sqrt(E/(...*SNR_bit_abs));

  %>> SOLUTION
  sigma = sqrt(E/(2*SNR_bit_abs));
  %<<


  %---------------------------------------
  % Generate N random symbols
  %---------------------------------------


  I = randint(1,N,M);                       % Uniform random number generator interval [0,M-1]


  %---------------------------------------
  % Gray-encoded the N symbols
  %---------------------------------------

  S = bin2gray(I,'pam',M);                  % Gray-encode random numbers and assign to Symbols
                                            % Gray coding for PAM is suitable for this example


  %
  % Symbols and signal constellation 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % Symbols Gray-encoded
  %
  %  0    1    3    2  
  % 00   01   11   10 
  %



  %---------------------------------------
  % Modulation
  %---------------------------------------


  % ** STEP 3 **------------------------- 
  % Modify modulation
  % Hint: comment out 2 lines
  % -------------------------------------  



  for ix = 1:N                              % Modulate all the symbols

    if     (S(ix) == 0), X(:,ix) = [ 1  0] ;      
    elseif (S(ix) == 1), X(:,ix) = [ 0  1] ;

%    elseif (S(ix) == 3), X(:,ix) = [-1  0] ;
%    elseif (S(ix) == 2), X(:,ix) = [ 0 -1] ;
    end;
  
  end


  %>> SOLUTION
  % Comment out the following lines above
  %    elseif (S(ix) == 3), X(:,ix) = [-1  0] ;
  %    elseif (S(ix) == 2), X(:,ix) = [ 0 -1] ;
  %<<



  %---------------------------------------
  % Add noise to symbols
  %---------------------------------------


  Noise = randn(2,N) * sigma;               %  Normally distributed noise

  Y = X + Noise;                            % add noise to received symbol




  %%%%%%%%%%%%
  % Receiver %
  %%%%%%%%%%%%
    
  %---------------------------------------
  % Detector
  %---------------------------------------

  for ix = 1:N
  
%    if     (Y(1,ix) > abs(Y(2,ix))), S_hat(ix) = 0;
%    elseif (Y(2,ix) > abs(Y(1,ix))), S_hat(ix) = 1;
%    elseif (Y(1,ix) <-abs(Y(2,ix))), S_hat(ix) = 3;
%    else                             S_hat(ix) = 2;    
%    end;
    

    % ** STEP 4 **------------------------- 
    % Comment out the decision areas above 
    % and use the following conditions
    % -------------------------------------  

    if ( (Y(2,ix) > abs(Y(1,ix))) || (Y(1,ix) <-abs(Y(2,ix))) ), S_hat(ix) = 1;
    else                                                         S_hat(ix) = 0;    
    end;


  
  end



  %---------------------------------------
  % Symbol error counter
  %---------------------------------------
 
  symbol_errors = 0;
  
  for ix = 1:N
  
    if (S(ix) ~= S_hat(ix))
      symbol_errors = symbol_errors + 1;
    end;
  
  end

  P_symbol_err_sim(iSNR) = symbol_errors / N;


  %---------------------------------------
  % Bit error counter
  %---------------------------------------

  bit_errors = 0;
  
  for ix = 1:N

    s1 = dec2binvec(S(ix),3);               % convert symbols expressed as decimal numbers
    s2 = dec2binvec(S_hat(ix),3);           %  to binary vector with 3 elements  

    bit_err_n = sum(xor(s1,s2));            % count how many bits are in error

    bit_errors = bit_errors  + bit_err_n;
  
  end

  P_bit_err_sim(iSNR) = bit_errors / (b*N);



  %---------------------------------------
  % Decode gray coded symbols
  %---------------------------------------

  I_hat = gray2bin(S_hat,'pam',M);


%---------------------------------------
% end loop SNR
%---------------------------------------

end


%---------------------------------------
% Plot symbol error prob.got by simulation
%---------------------------------------
  
figure(1)                                   % plot on figure 1
semilogy(SNR_bit_dB,P_symbol_err_sim,'r:*');% Y axis in the plot is logarithmic
grid on                                     % show grid



%---------------------------------------
% Plot bit error prob.got by simulation
%---------------------------------------

hold on
semilogy(SNR_bit_dB,P_bit_err_sim,'b:s'); 
hold off

legend('SER','BER')
title('Symbol and bit error probabilities.')
xlabel('Eb/No')
ylabel('Symbol/Bit Error Probability')



% ** STEP 5 **------------------------- 
% Observe symbol and bit error rates
%
% -------------------------------------  


% ** STEP 6 **------------------------- 
% Plot theoretical bit error probabilities
% orthogonal and biorthogonal
% Comment out the code below
% Observe 3dB difference.
% -------------------------------------  


%---------------------------------------
% Theoretical bit error probability 
%---------------------------------------

SNR_bit_dB_t = 0:0.2:12;

for iSNR = 1:length(SNR_bit_dB_t)

  SNR_bit_abs = 10^(SNR_bit_dB_t(iSNR)/10);

    % Symbol Error Probability for BPSK = Bit Error Probability
    P_symbol_err_theo_orthogonal(iSNR) = qfunc(sqrt(SNR_bit_abs));
    P_symbol_err_theo_biorthogonal(iSNR) = qfunc(sqrt(2*SNR_bit_abs));

end

hold on
semilogy(SNR_bit_dB_t,P_symbol_err_theo_orthogonal,'k'); 
semilogy(SNR_bit_dB_t,P_symbol_err_theo_biorthogonal,'k--'); 
hold off

legend('SER','BER','BER Theoretical Orthogonal','BER Theoretical Biorthogonal')

