

%=======================================================================================
%    QPSK/BPSK simulation
%
%    TKK - Comlab 2007 - SL
%=======================================================================================
%
% Purpose: simulate a communication system using QPSK/BPSK modulation
%          and calculate bit and symbol error probabilities.
%
% Channel: AWGN
%
%=======================================================================================

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

clear all


% ** STEP 1 **------------------------- 
% Set modulation to BPSK
% 
% -------------------------------------

%Modulation = 'B';                           % Modulation is BPSK


% ** STEP 5 **------------------------- 
% Change modulation to QPSK
% 
% -------------------------------------

%Modulation = 'Q';                           % Modulation is QPSK



if     (Modulation == 'B'), M = 2;          % BPSK case
elseif (Modulation == 'Q')  M = 4;          % QPSK case
end 


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

N = 1000;                                  % number of symbols to transmit 

E = 1;                                      % energy


SNR_bit_dB_MAX = 10;                        % 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)

% iSNR = 3;                                 % Simulation for only one SNR


  SNR_bit_abs = 10^(SNR_bit_dB(iSNR)/10);   % transform SNR in dB to abs values
    
  if     (Modulation == 'B')
    sigma = sqrt(E/(2*SNR_bit_abs));        % standard deviation BPSK case    
  elseif (Modulation == 'Q')  
    sigma = sqrt(E/(4*SNR_bit_abs));        % standard deviation QPSK case
  end 



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


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

  % ** STEP 2 **------------------------- 
  % Gray-encode symbols, use bin2gray()
  % and modulation psk
  % -------------------------------------

  %S = bin2gray(....);                 


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


  if     (Modulation == 'B')

    for ix = 1:N                            % Modulate all the symbols BPSK
      if     (S(ix) == 0), X(ix) =  1;      
      elseif (S(ix) == 1), X(ix) = -1;      
      end;
    end

  elseif (Modulation == 'Q')  

    for ix = 1:N                            % Modulate all the symbols QPSK
      if     (S(ix) == 0), X(ix) =  1+j*0;      
      elseif (S(ix) == 1), X(ix) =  0+j*1;      
      elseif (S(ix) == 3), X(ix) = -1+j*0;      
      elseif (S(ix) == 2), X(ix) =  0-j*1;      
      end;
    end

  end 


                                            %.....................................................
                                            %
                                            % % ** STEP 9 ** 
                                            % Optional.
                                            %
                                            % An alternative approach to modulete PSK is using
                                            % function pskmod() from the Communications Toolbox.
                                            % For QPSK case replace the lines above by 
                                            % X = pskmod(S,4,0,'gray');
                                            % and observe that the modulation is the same.
                                            %.....................................................

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


  %---------------------------------------
  % Rotating constellation in the channel
  %---------------------------------------


  if     (Modulation == 'Q')

    X_orig = X;                               % keep a copy of the received symbols
  
  
    R=5+i*8;                                  % rotation of constellation
    R=R/abs(R);



    % ** STEP 6 **------------------------- 
    % Introduce a rotation in the channel
    % 
    % -------------------------------------


    % X = ....


  end

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

  if     (Modulation == 'B')
    Noise = sigma*(randn(1,N));             % BPSK case
  elseif (Modulation == 'Q')  
    Noise = sigma*(randn(1,N)+j*randn(1,N));% QPSK case      
  end 

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

  Y_noise = Y;
    
  %%%%%%%%%%%%
  % Receiver %
  %%%%%%%%%%%%


  %---------------------------------------
  % Equalize the channel
  %---------------------------------------

  % ** STEP 7 **------------------------- 
  % Equalize the channel by multiplying
  % the received symbols with the channel 
  % complex conjugate
  % -------------------------------------

  if     (Modulation == 'Q')

    %Y = ...*Y;
  

  end  
    
  %---------------------------------------
  % Demodulator/Detector
  %---------------------------------------


  % ** STEP 3 **------------------------- 
  % Demodulate PSK by using pskdemod()
  % 
  % -------------------------------------


  %S_hat = pskdemod(....);




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

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

    
  %---------------------------------------
  % Symbol error counter
  %---------------------------------------

  symbol_errors = 0;
  
  for ix = 1:N

    if (I(ix) ~= I_hat(ix))
      symbol_errors = symbol_errors + 1;
    end;
  
  end

  P_symbol_err_sim(iSNR) = symbol_errors/N;

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

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

  b1 = de2bi(S,2);                          % Convert decimal numbers to binary vectors
  b2 = de2bi(S_hat,2);                      % Convert decimal numbers to binary vectors  
 
  bit_errors = sum(sum(b1~=b2));


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


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

end;



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


  if     (Modulation == 'B')

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

  elseif (Modulation == 'Q')  

    % approximation of Symbol Error Probability for M > 4
    P_symbol_err_theo(iSNR) = 2*qfunc(sin(pi/M)*sqrt(2*k*SNR_bit_abs));

  end 

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) = qfunc(sqrt(2*SNR_bit_abs));

end




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



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

figure(1)                                   % plot on figure 1
semilogy(SNR_bit_dB_t,P_symbol_err_theo,'b'); 
grid on                                     % show grid

%if     (Modulation == 'B')
%  title('Symbol/Bit error probabilities. BPSK ')
%elseif (Modulation == 'Q')  
%  title('Symbol/Bit error probabilities. QPSK ')
%end 


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



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

figure(1)                                   % plot on figure 2
hold on                                     % keep previous plot on the same figure
semilogy(SNR_bit_dB_t,P_bit_err_theo,'b--');
hold off
grid on

if     (Modulation == 'B')
  title('Symbol/Bit error probabilities. BPSK')
elseif (Modulation == 'Q')  
  title('Symbol/Bit error probabilities. QPSK')
end 

xlabel('Eb/No')
ylabel('Bit Error Probability')
 
%---------------------------------------
% Plot bit error prob.got by simulation
%---------------------------------------

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


legend('SER theoretical','SER simulation','BER theoretical','BER simulation')
    

% ** STEP 8 **------------------------- 
% Observe the rotation and noise 
% in the the constellation
% for the last received symbols
% -------------------------------------

if     (Modulation == 'Q')


  figure(3)
  plot(real(X_orig),imag(X_orig),'k*')
  grid
  title('Signal Constellation')
  xlabel('real')
  ylabel('imag')

  figure(4)
  plot(real(X),imag(X),'k*')
  grid
  title('Rotated Signal Constellation')
  xlabel('real')
  ylabel('imag')

  figure(5)
  plot(real(Y_noise ),imag(Y_noise),'k*')
  grid
  title('Rotated Signal Constellation with Noise')
  xlabel('real')
  ylabel('imag')

  figure(6)
  plot(real(Y),imag(Y),'k*')
  grid
  title('Signal Constellation after Equalization')
  xlabel('real')
  ylabel('imag')

end 




