awgn_vd_spc.m

详细讲述纠错码的书籍

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