

%=======================================================================================
%    8 PAM simulation
%
%    TKK - Comlab 2007 - SL
%=======================================================================================
%
% Purpose: simulate a 8 level PAM communication system and
%   calculate bit and symbol error probabilities.
%   Simulation results are compared against theoretical
%   error probabilities.
%
% Channel: AWGN
%
%=======================================================================================

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

clear all

M = 8;                                      % modulation is 8 PAM 

b = 3;                                      % number of bits to represent 8 levels is 3

N = 3000;                                   % number of symbols to transmit 

d = 1;                                      % distance

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
%----------------------------------------

  % ** STEP 7 **------------------------- 
  % Repeat the simulation for different 
  % SNR. Uncomment the line below (for)
  % and 'end' on line 274
  % -------------------------------------  

 iSNR = 3;                                 % Simulation for only one SNR
  
%for iSNR = 1:length(SNR_bit_dB)



  % ** STEP 1 **------------------------- 
  % transform SNR [dB] stored in variable
  % SNR_bit_dB(iSNR) to SNR [abs]
  % -------------------------------------
  
  % SNR_bit_abs =  ....SNR_bit_dB(iSNR)....


  % ** STEP 2 **------------------------- 
  % complete standard deviation for
  % 8 PAM using SNR_bit_abs
  % -------------------------------------  
  
  % sigma = ....SNR_bit_abs.... 


  %%%%%%%%%%%%%%%
  % Transmitter %
  %%%%%%%%%%%%%%%

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


 
  % ** STEP 3 **------------------------- 
  % generate integer random number in the
  % interval [0,7] and verify the 
  % distribution with a histogram.
  % use randint() 
  % -------------------------------------  
  
  % I = randint(....);


  % observe the histogram of the random 
  % numbers and check that the distribution
  % of symbols is indeed in the on the
  % interval [0,7]

  % histogram_interval = -2:0.5:12;
  % hist(I,histogram_interval)



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

  S = bin2gray(I,'pam',M);                  % Gray-encode random numbers and assign to Symbols

  %
  % Symbols and signal constellation 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % Original symbol
  %
  %  0    1    2    3    4    5    6    7    
  %
  % Symbols Gray-encoded
  %
  %  0    1    3    2    6    7    5    4
  % 000  001  011  010  110  111  101  100
  %
  % Signal constellation 
  %
  % -7d  -5d  -3d   -d   d    3d   5d   7d
  %



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

  for ix = 1:N                              % Modulate all the symbols

    if     (S(ix) == 0), X(ix) = -7*d ;      
    elseif (S(ix) == 1), X(ix) = -5*d ;
    elseif (S(ix) == 3), X(ix) = -3*d ;
    elseif (S(ix) == 2), X(ix) = -1*d ;
    elseif (S(ix) == 6), X(ix) =  1*d ;
    elseif (S(ix) == 7), X(ix) =  3*d ;
    elseif (S(ix) == 5), X(ix) =  5*d ;
    elseif (S(ix) == 4), X(ix) =  7*d ;
    end;
  
  end
                                            %.....................................................
                                            %
                                            % % ** STEP  10 ** 
                                            % Optional.
                                            %
                                            % An alternative approach to modulete PAM is using
                                            % function pammod() from the Communications Toolbox.
                                            % Replace the lines above by 
                                            % X = pammod(S,M,0,'gray') .* d;
                                            % and observe that the modulation is the same.
                                            %.....................................................


  %---------------------------------------
  % average transmitted signal energy
  %---------------------------------------


  % ** STEP  9 **------------------------- 
  % Observe the average energy
  % -------------------------------------  

  % SNR_bit_dB(iSNR)
  % AvgEnergy = Y*Y'/N



  %%%%%%%%%%%
  % Channel %
  %%%%%%%%%%%

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


  % ** STEP 4 **------------------------- 
  % Add noise to symbols, use randn()
  % -------------------------------------  
  
  
  %Noise = randn(....).....
  %Y = X + Noise;                           % add noise to received symbol



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

  for ix = 1:N
  
    if     (Y(ix) < -6*d), S_hat(ix) = 0;
    elseif (Y(ix) < -4*d), S_hat(ix) = 1;
    elseif (Y(ix) < -2*d), S_hat(ix) = 3;
    elseif (Y(ix) <  0*d), S_hat(ix) = 2;
    elseif (Y(ix) <  2*d), S_hat(ix) = 6;  
    elseif (Y(ix) <  4*d), S_hat(ix) = 7;  
    elseif (Y(ix) <  6*d), S_hat(ix) = 5;  
    else                   S_hat(ix) = 4;    
    end;
  
  end


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

  I_hat = gray2bin(S_hat,'pam',M);
  
  %---------------------------------------
  % Symbol error counter
  %---------------------------------------
  
  symbol_errors = 0;
  
  for ix = 1:N
  
    if (I(ix) ~= I_hat(ix))
      symbol_errors = symbol_errors + 1;
    end;
  
  end
  

  % ** STEP 5 **------------------------- 
  % Calculate the symbol error probability
  % -------------------------------------  

  %P_symbol_err_sim(iSNR) = symbol_errors...


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

  % Gray coded symbols are used to calculate
  % the bits in error

  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

  % ** STEP 6 **------------------------- 
  % Calculate the bit error probability
  % -------------------------------------  

  % P_bit_err_sim(iSNR) = bit_errors.....


                                            %.....................................................
                                            % Comment: an alternative to dec2binvec() is de2bi()
                                            %.....................................................
  
%---------------------------------------
% end loop SNR
%---------------------------------------

end;



% ** STEP  8 **------------------------- 
% Plot using the code below.
% Observe the error probabilities
% Observe constellation and histogram of
% the last samples in memory.
% -------------------------------------  


%---------------------------------------
% Theoretical symbol error probability
%---------------------------------------

SNR_bit_dB_t = 0:0.2:SNR_bit_dB_MAX;

for iSNR = 1:length(SNR_bit_dB_t)

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

  P_symbol_err_theo(iSNR) = (2*(M-1)/M)*qfunc(sqrt((6*log2(M)/(M^2-1))*SNR_bit_abs));

end


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

SNR_bit_dB_t = 0:0.2:SNR_bit_dB_MAX;

for iSNR = 1:length(SNR_bit_dB_t)

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

  P_bit_err_theo(iSNR) = (1/3) * (2*(M-1)/M)*qfunc(sqrt((6*log2(M)/(M^2-1))*SNR_bit_abs));

end



%---------------------------------------
% Plot Theoretical symbol error prob.
%---------------------------------------

figure(1)                                   % plot on figure 1
semilogy(SNR_bit_dB_t,P_symbol_err_theo); 
grid on                                     % show grid
title('Symbol error probabilities. 8 PAM')
xlabel('Eb/No')
ylabel('Symbol Error Probability')


%---------------------------------------
% Plot symbol error prob.got by simulation
%---------------------------------------
  
figure(1)                                   % plot on figure 1
hold on                                     % keep previous plot on the same figure
semilogy(SNR_bit_dB,P_symbol_err_sim,'r*'); % Y axis in the plot is logarithmic
hold off
legend('theoretical','simulation')


%---------------------------------------
% Plot Theoretical bit error prob.
%---------------------------------------

figure(2)                                   % plot on figure 2
semilogy(SNR_bit_dB_t,P_bit_err_theo);
grid on
title('Bit error probabilities. 8 PAM')
xlabel('Eb/No')
ylabel('Bit Error Probability')
 
%---------------------------------------
% Plot bit error prob.got by simulation
%---------------------------------------

figure(2)                                   % plot on figure 2
hold on
semilogy(SNR_bit_dB,P_bit_err_sim,'r*'); 
hold off
legend('theoretical','simulation')



%---------------------------------------
% Plot Constellation of the last samples
%---------------------------------------

figure(3)
plot(Y,1,'k*')
grid
title('Signal Constellation')

%---------------------------------------
% Plot Histogram of the last samples
%---------------------------------------

figure(4)
histogram_interval = -8:0.1:8;
hist(Y,histogram_interval);
title('Histogram')

%[n,xout] = hist(Y,histogram_interval);
%bar(xout,n./N)
 