%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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 % Convolutional code of Rate-1/2, Generator polynomial - [7,5] octal % Hard decision and soft decison decoding is used. clear N = 10^6 ;% number of bits or symbols tic Eb_N0_dB = [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 = 15; decodeDel = 50; for yy = 1:length(Eb_N0_dB) % 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(yy)/20)*n; % additive white gaussian noise % receiver cipHard = real(y)>0; % hard decision cipSoft = real(y); % soft decision % Viterbi decoding pmHard = zeros(4,1); % hard path metric svHard_v = zeros(4,length(y)/2); % hard survivor path pmSoft = zeros(4,1); % soft path metric svSoft_v = zeros(4,length(y)/2); % soft survivor path count=1; ipHatHard_v = zeros(1,length(y)/2); ipHatSoft_v = zeros(1,length(y)/2); svHard_v = zeros(4,length(y)/2); svSoft_v = zeros(4,length(y)/2); for ii = 1:length(y)/2 rHard = cipHard(2*ii-1:2*ii); % taking 2 hard bits rSoft = cipSoft(2*ii-1:2*ii); % taking 2 soft bits % computing the Hamming distance and euclidean distance rHardv = kron(ones(4,1),rHard); rSoftv = kron(ones(4,1),rSoft); hammingDist = sum(xor(rHardv,refHard),2); euclideanDist = sum(rSoftv.*refSoft,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; bm1Soft = pmSoft(1,1) + euclideanDist(1); pmSoft_n(1,1) = bm1Soft; svSoft(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; bm1Soft = pmSoft(3,1) + euclideanDist(3); pmSoft_n(2,1) = bm1Soft; svSoft(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; bm1Soft = pmSoft(1,1) + euclideanDist(4); pmSoft_n(3,1) = bm1Soft; svSoft(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; bm1Soft = pmSoft(3,1) + euclideanDist(2); pmSoft_n(4,1) = bm1Soft; svSoft(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; bm1Soft = pmSoft(1,1) + euclideanDist(1); bm2Soft = pmSoft(2,1) + euclideanDist(4); [pmSoft_n(1,1) idx] = min([bm1Soft,bm2Soft]); svSoft(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; bm1Soft = pmSoft(3,1) + euclideanDist(3); bm2Soft = pmSoft(4,1) + euclideanDist(2); [pmSoft_n(2,1) idx] = min([bm1Soft,bm2Soft]); svSoft(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; bm1Soft = pmSoft(1,1) + euclideanDist(4); bm2Soft = pmSoft(2,1) + euclideanDist(1); [pmSoft_n(3,1) idx] = min([bm1Soft,bm2Soft]); svSoft(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; bm1Soft = pmSoft(3,1) + euclideanDist(2); bm2Soft = pmSoft(4,1) + euclideanDist(3); [pmSoft_n(4,1) idx] = min([bm1Soft,bm2Soft]); svSoft(4,1) = idx+2; end pmHard = pmHard_n; svHard_v(:,ii) = svHard; pmSoft = pmSoft_n; svSoft_v(:,ii) = svSoft; % trace back unit if ( (mod(ii,count*decodeDel+TB) == 0) || (ii==(length(y)/2)) ) if (ii == length(y)/2) % starting with state 00 when the packet end is reached currHardState = 1; currSoftState = 1; stopTB = startTB+1; startTB = length(y)/2; stopDecode = startDecode+1; startDecode = startTB; else % finding the minimum pathmetric to start the traceback [ dum1 currHardState ] = min(pmHard); [ dum2 currSoftState ] = min(pmSoft); startTB = count*decodeDel+TB; stopTB = startTB-TB+1; startDecode = stopTB-1; stopDecode = startDecode - decodeDel+1; end if (ii ~= (length(y)/2)) % start and stop of traceback for jj = startTB:-1:stopTB prevHardState = svHard_v(currHardState,jj); currHardState = prevHardState; prevSoftState = svSoft_v(currSoftState,jj); currSoftState = prevSoftState; end end % start and stop of decoding + traceback for tt = startDecode:-1:stopDecode prevHardState = svHard_v(currHardState,tt); ipHatHard_v(tt) = ipLUT(currHardState,prevHardState); currHardState = prevHardState; prevSoftState = svSoft_v(currSoftState,tt); ipHatSoft_v(tt) = ipLUT(currSoftState,prevSoftState); currSoftState = prevSoftState; end count = count+1; end end % counting the errors nErrHardViterbi(yy) = size(find([ip - ipHatHard_v(1:N)]),2); nErrSoftViterbi(yy) = size(find([ip - ipHatSoft_v(1:N)]),2); end simBer_HardViterbi = nErrHardViterbi/N; % simulated ber - hard decision Viterbi decoding BER simBer_SoftViterbi = nErrSoftViterbi/N; % simulated ber - soft decision Viterbi decoding BER theoryBer = 0.5*erfc(sqrt(10.^(Eb_N0_dB/10))); % theoretical ber uncoded AWGN close all figure semilogy(Eb_N0_dB,theoryBer,'bd-','LineWidth',2); hold on semilogy(Eb_N0_dB,simBer_HardViterbi,'mp-','LineWidth',2); semilogy(Eb_N0_dB,simBer_SoftViterbi,'cd-','LineWidth',2); axis([0 10 10^-5 0.5]) grid on legend('theory - uncoded', 'simulation - hard Viterbi', 'simulation - soft Viterbi'); xlabel('Eb/No, dB'); ylabel('Bit Error Rate'); title('BER for BCC with Viterbi decoding (TB=15) for BPSK in AWGN'); print('ber_plot_bpsk_viterbi_decode_finite_traceback_depth.png','-dpng','-S448,336') toc