%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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 : 14th January 2009 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Script for computing BER with Binary Convolutional Code % and Viterbi decoding. % 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 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 ]; 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 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; end % trace back unit currHardState = 1; currSoftState = 1; ipHatHard_v = zeros(1,length(y)/2); ipHatSoft_v = zeros(1,length(y)/2); for jj = length(y)/2:-1:1 prevHardState = svHard_v(currHardState,jj); ipHatHard_v(jj) = ipLUT(currHardState,prevHardState); currHardState = prevHardState; prevSoftState = svSoft_v(currSoftState,jj); ipHatSoft_v(jj) = ipLUT(currSoftState,prevSoftState); currSoftState = prevSoftState; 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 for BPSK in AWGN');