%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All rights reserved by Krishna Pillai, http://www.dsplog.com
% The file may not be re-distributed without explicit authorization
% from Krishna Pillai.
% Checked for proper operation with Octave Version 3.0.0
% Author        : Krishna Pillai
% Email         : krishna@dsplog.com
% Version       : 1.0
% Date          : 26th July 2009
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Script for computing BER with Binary Convolutional Code 
% and Viterbi decoding with finite survivor state memory and
% variable traceback depth 
% Convolutional code of Rate-1/2, Generator polynomial - [7,5] octal 
% Hard decision decoding is used. 

clear
N = 10^4 ;% number of bits or symbols
tic
Eb_N0_dB = 4; % [0:1:10]; % multiple Eb/N0 values
Ec_N0_dB = Eb_N0_dB - 10*log10(2);

refHard = [0 0 ; 0 1; 1 0  ; 1 1 ];
refSoft = -1*[-1 -1; -1 1 ;1 -1; 1 1 ];

ipLUT = [ 0   0   0   0;...
	  0   0   0   0;...
          1   1   0   0;...
	  0   0   1   1 ];
TB_v = [1:2:21] ;
decodeDel = 1;

for yy = 1:length(TB_v)
  TB = TB_v(yy);

   % Transmitter
   ip = rand(1,N)>0.5; % generating 0,1 with equal probability

   % convolutional coding, rate - 1/2, generator polynomial - [7,5] octal
   cip1 = mod(conv(ip,[1 1 1 ]),2);
   cip2 = mod(conv(ip,[1 0 1 ]),2);
   cip = [cip1;cip2];
   cip = cip(:).';

   s = 2*cip-1; % BPSK modulation 0 -> -1; 1 -> 0 

   n = 1/sqrt(2)*[randn(size(cip)) + j*randn(size(cip))]; % white gaussian noise, 0dB variance 

   % Noise addition
   y = s + 10^(-Ec_N0_dB/20)*n; % additive white gaussian noise

   % receiver 
   cipHard = real(y)>0; % hard decision

   % Viterbi decoding
   pmHard  = zeros(4,1);  % hard path metric
   svHard_v  = zeros(4,length(y)/2); % hard survivor path
   
   count=1; 
   ipHatHard_v = zeros(1,length(y)/2);
   svHard_v    = zeros(4,length(y)/2); 
   for ii = 1:length(y)/2
      rHard = cipHard(2*ii-1:2*ii); % taking 2 hard bits
      
      % computing the Hamming distance and euclidean distance
      rHardv = kron(ones(4,1),rHard);
      hammingDist = sum(xor(rHardv,refHard),2);
 
      if (ii == 1) || (ii == 2) 

         % branch metric and path metric for state 0
         bm1Hard = pmHard(1,1) + hammingDist(1);
         pmHard_n(1,1)  = bm1Hard; 
         svHard(1,1)  = 1; 
 
         % branch metric and path metric for state 1
         bm1Hard = pmHard(3,1) + hammingDist(3);
         pmHard_n(2,1) = bm1Hard;
         svHard(2,1)  = 3; 
         

         % branch metric and path metric for state 2
         bm1Hard = pmHard(1,1) + hammingDist(4);
         pmHard_n(3,1) = bm1Hard;
         svHard(3,1)  = 1; 

         % branch metric and path metric for state 3
         bm1Hard = pmHard(3,1) + hammingDist(2);
         pmHard_n(4,1) = bm1Hard;
         svHard(4,1)  = 3; 

      else
         % branch metric and path metric for state 0
         bm1Hard = pmHard(1,1) + hammingDist(1);
         bm2Hard = pmHard(2,1) + hammingDist(4);
         [pmHard_n(1,1) idx] = min([bm1Hard,bm2Hard]);
         svHard(1,1)  = idx; 
 
         % branch metric and path metric for state 1
         bm1Hard = pmHard(3,1) + hammingDist(3);
         bm2Hard = pmHard(4,1) + hammingDist(2);
         [pmHard_n(2,1) idx] = min([bm1Hard,bm2Hard]);
         svHard(2,1)  = idx+2; 

         % branch metric and path metric for state 2
         bm1Hard = pmHard(1,1) + hammingDist(4);
         bm2Hard = pmHard(2,1) + hammingDist(1);
         [pmHard_n(3,1) idx] = min([bm1Hard,bm2Hard]);
         svHard(3,1)  = idx; 

         % branch metric and path metric for state 3
         bm1Hard = pmHard(3,1) + hammingDist(2);
         bm2Hard = pmHard(4,1) + hammingDist(3);
         [pmHard_n(4,1) idx] = min([bm1Hard,bm2Hard]);
         svHard(4,1)  = idx+2; 

      end

   pmHard = pmHard_n; 
   svHard_v(:,ii) = svHard;
   
   % trace back unit
   if ( (mod(ii,count*decodeDel+TB) == 0) || (ii==(length(y)/2)) )
   
      if (ii == length(y)/2)
         currHardState0 = 1;
         currHardState_min = 1;
	 stopTB  = startTB+1;
         startTB = length(y)/2;
         stopDecode  = startDecode+1;
         startDecode = startTB;
      else
         [ dum1 currHardState_min ] = min(pmHard);
         currHardState0 = 1;
         startTB = count*decodeDel+TB;
         stopTB  = startTB-TB+1;
         startDecode = stopTB-1;
         stopDecode  = startDecode - decodeDel+1;
      end

    if (ii ~= (length(y)/2))
       for jj = startTB:-1:stopTB
           prevHardState0   = svHard_v(currHardState0,jj); 
           currHardState0   = prevHardState0;

           prevHardState_min   = svHard_v(currHardState_min,jj); 
           currHardState_min   = prevHardState_min;

        end
    end
     
     for tt = startDecode:-1:stopDecode
        prevHardState0   = svHard_v(currHardState0,tt); 
        ipHatHard_v0(tt) = ipLUT(currHardState0,prevHardState0);
        currHardState0   = prevHardState0;

        prevHardState_min   = svHard_v(currHardState_min,tt); 
        ipHatHard_v_min(tt) = ipLUT(currHardState_min,prevHardState_min);
        currHardState_min   = prevHardState_min;
     end

     count = count+1;
    end
  
  end
   
   % counting the errors
   nErrHardViterbi0(yy) = size(find([ip - ipHatHard_v0(1:N)]),2);
   nErrHardViterbi_min(yy) = size(find([ip - ipHatHard_v_min(1:N)]),2);

end

simBer_HardViterbi0 = nErrHardViterbi0/N;
simBer_HardViterbi_min = nErrHardViterbi_min/N;

close all
figure
semilogy(TB_v,simBer_HardViterbi0,'sb-','LineWidth',2);
hold on
semilogy(TB_v,simBer_HardViterbi_min,'mp-','LineWidth',2);
grid on
axis([0 25 10^-3 10^0])
xlabel('traceback depth, TB')
ylabel('Bit Error Rate');
title('BER with Viterbi decoding vs traceback at Eb/N0=4dB')
legend('start = state_{00}', 'start = min_{path metric}')
toc