% File: AWGN_VD_SPC.m
% Description:
% Simulation of Viterbi decoding of a (5,4,2) SPC code based on its
% syndrome trellis:
%
% 0 --- 0 --- 0 --- 0 --- 0 --- 0
% \ \ / \ / \ / /
% \ x x x /
% \ / \ / \ / \ /
% 1 --- 1 --- 1 --- 1
%
% Copyright (c) 2007. Robert Morelos-Zaragoza. All rights reserved.
clear all
n=5; % Code length
k=4; % Code dimension
rate=k/n; % Code rate
G = [ 1 0 0 0 1 % Generator matrix (for encoding)
0 1 0 0 1
0 0 1 0 1
0 0 0 1 1 ];
W = [ 1 0 10 0 5 0]; % Weight distribution
nsim = 50000; % Number of codewords per SNR value
init = 0; % Inital value of Eb/No (dB)
inc = 0.5; % Increment in Eb/No (dB)
final = 6; % Final value of Eb/No (dB)
num=0;
seed=123445678;
rand('state',seed);
randn('state',seed);
fprintf('Eb/No(dB) \tBER_HD \tErr_HD \t BER_ML \tErr_ML \n');
for SNRdB = init:inc:final
error_HD = 0;
error_ML = 0;
TotalNumError = 0;
TotalN = 0;
SNRdBs = SNRdB + 10*log(rate)/log(10);
No = (10.^(-SNRdBs/10));
Var = No/2;
for Ns = 1:nsim
msg = randint(1,k,[0,1]); % Random binary message vector
codeword = mod(msg*G,2); % Compute the codeword in C
x = (-1).^codeword; % BPSK mapping
rec = x + sqrt(Var)*randn(1,n); % Add AWGN to obtain received vector
% ------------------------------------------------------------
% Viterbi decoding with HARD DECISIONS
% ------------------------------------------------------------
rh = (1-sign(rec))/2; % Hard decision vector
% First trellis stage
M0(1) = rh(1); % State 0
M1(1) = 1-rh(1); % State 1
% ======== Trellis stage 2 ========
BM0 = rh(2); BM1 = 1 - BM0; % Branch metrics (Hamming distance)
% State 0
if M0(1)+BM0 < M1(1)+BM1
M0(2) = M0(1)+BM0; v0(1)=0; v0(2)=0;
else
M0(2) = M1(1)+BM1; v0(1)=1; v0(2)=1;
end
% State 1
if M0(1)+BM1 < M1(1)+BM0
M1(2) = M0(1)+BM1; v1(1)=0; v1(2)=1;
else
M1(2) = M1(1)+BM0; v1(1)=1; v1(2)=0;
end
vhat0=v0; vhat1=v1; % Path update
% ======== Trellis stage 3 ========
BM0 = rh(3); BM1 = 1 - BM0; % Branch metrics (Hamming distance)
% State 0
if M0(2)+BM0 < M1(2)+BM1
M0(3) = M0(2)+BM0; v0(3)=0;
else
M0(3) = M1(2)+BM1; v0(1:2) = vhat1(1:2); v0(3)=1;
end
% State 1
if M0(2)+BM1 < M1(2)+BM0
M1(3) = M0(2)+BM1; v1(1:2) = vhat0(1:2); v1(3)=1;
else
M1(3) = M1(2)+BM0; v1(3)=0;
end
vhat0=v0; vhat1=v1; % Path update
% ======== Trellis stage 4 ========
BM0 = rh(4); BM1 = 1 - BM0; % Branch metrics (Hamming distance)
% State 0
if M0(3)+BM0 < M1(3)+BM1
M0(4) = M0(3)+BM0; v0(4)=0;
else
M0(4) = M1(3)+BM1; v0(1:3) = vhat1(1:3); v0(4)=1;
end
% State 1
if M0(3)+BM1 < M1(3)+BM0
M1(4) = M0(3)+BM1; v1(1:3) = vhat0(1:3); v1(4)=1;
else
M1(4) = M1(3)+BM0; v1(4)=0;
end
vhat0=v0; vhat1=v1; % Path update
% ======== Trellis stage 5 (last) ========
BM0 = rh(5); BM1 = 1 - BM0; % Branch metrics (Hamming distance)
% State 0
if M0(4)+BM0 < M1(4)+BM1
M0(5) = M0(4)+BM0; v0(5)=0;
codehard = v0;
else
M0(5) = M1(4)+BM1; v0(1:4) = vhat1(1:4); v0(5)=1;
codehard = v0;
end
% ------------------------------------------------------------
% Viterbi decoding with SOFT DECISIONS
% ------------------------------------------------------------
% First trellis stage
M0(1) = rec(1); % State 0
M1(1) = -rec(1); % State 1
% ======== Trellis stage 2 ========
BM0 = rec(2); BM1 = -BM0; % Branch metrics (correlation)
% State 0
if M0(1)+BM0 > M1(1)+BM1
M0(2) = M0(1)+BM0; v0(1)=0; v0(2)=0;
else
M0(2) = M1(1)+BM1; v0(1)=1; v0(2)=1;
end
% State 1
if M0(1)+BM1 > M1(1)+BM0
M1(2) = M0(1)+BM1; v1(1)=0; v1(2)=1;
else
M1(2) = M1(1)+BM0; v1(1)=1; v1(2)=0;
end
vhat0=v0; vhat1=v1; % Path update
% ======== Trellis stage 3 ========
BM0 = rec(3); BM1 = -BM0; % Branch metrics (correlation)
% State 0
if M0(2)+BM0 > M1(2)+BM1
M0(3) = M0(2)+BM0; v0(3)=0;
else
M0(3) = M1(2)+BM1; v0(1:2) = vhat1(1:2); v0(3)=1;
end
% State 1
if M0(2)+BM1 > M1(2)+BM0
M1(3) = M0(2)+BM1; v1(1:2) = vhat0(1:2); v1(3)=1;
else
M1(3) = M1(2)+BM0; v1(3)=0;
end
vhat0=v0; vhat1=v1; % Path update
% ======== Trellis stage 4 ========
BM0 = rec(4); BM1 = -BM0; % Branch metrics (correlation)
% State 0
if M0(3)+BM0 > M1(3)+BM1
M0(4) = M0(3)+BM0; v0(4)=0;
else
M0(4) = M1(3)+BM1; v0(1:3) = vhat1(1:3); v0(4)=1;
end
% State 1
if M0(3)+BM1 > M1(3)+BM0
M1(4) = M0(3)+BM1; v1(1:3) = vhat0(1:3); v1(4)=1;
else
M1(4) = M1(3)+BM0; v1(4)=0;
end
vhat0=v0; vhat1=v1; % Path update
% ======== Trellis stage 5 (last) ========
BM0 = rec(5); BM1 = -BM0; % Branch metrics (correlation)
% State 0
if M0(4)+BM0 > M1(4)+BM1
M0(5) = M0(4)+BM0; v0(5)=0;
v = v0;
else
M0(5) = M1(4)+BM1; v0(1:4) = vhat1(1:4); v0(5)=1;
v = v0;
end
% Enumerate decoding __information__ errors and add them
error_HD = error_HD + sum(xor(codeword(1:k),codehard(1:k)));
error_ML = error_ML + sum(xor(codeword(1:k),v(1:k)));
end % for Ns
num = num + 1;
ber_HD(num) = error_HD/(nsim*k);
ber_ML(num) = error_ML/(nsim*k);
snr(num) = SNRdB;
fprintf('%5.2f\t%e %6.0f\t%e %6.0f\n', snr(num), ber_HD(num), ...
error_HD, ber_ML(num), error_ML);
end % for SNRdB
% Compute the union bound for ML decoding
EbNo = 0:0.5:6; EbNoratio = 10.^(EbNo/10);
bound = 0;
for i=1:n
bound = bound + i*W(i+1)/n * Q(sqrt(2*i*(k/n)*EbNoratio));
end
% and plot results
semilogy(snr,ber_HD,'-r^'), hold on, semilogy(snr,ber_ML,'-bs')
semilogy(EbNo,bound,'-bo'), axis tight, grid on
p=Q(sqrt(2*EbNoratio)); plot(EbNo,1-(1-p).^5)
xlabel('E_b/N_0 (dB)'), ylabel('BER'), legend('HD sim','ML sim',...
'Union bound','HD bound')
评论0