bp_answer.m

clear;

pchkArray = [1, 1, 1, 0, 1, 0, 0;
             0, 1, 1, 1, 0, 1, 0;
             1, 0, 1, 1, 0, 0, 1];      % parity check matrix for (7,4) hamming code

noiseCodeword =  [1, 0, 1, 1, 0, 0, 0];      % input code word with 1 bit error
decodeword = zeros(1, numel(noiseCodeword));                  % array for storing decoded word
edgeArray = [eye(7); pchkArray];        % build the edge array which stores the connection between variable nodes and factor nodes
[i, j, s] = find(edgeArray);            % find the index of edgeArray;
edgeMap = sparse(i, j, [1:length(s)]);  % map edge id to node id;

max_iter = 10;                          % max number of iteration for bp;
p = 0.01;                               % error probability for BSC

[m,n] = size(pchkArray);                % get number of factor nodes and variable nodes, where the num of factor nodes is equal to n + m, and num of factor nodes is n;
numOfEdge = sum(edgeArray(:));          % get number of edges between variable nodes and factor nodes
numOfLabel = 2;                         % set the num of bp labels, where 2 means we only use 0 and 1 as the label in our bp algorithm.
LabelValue = [0, 1];                    % set the value
messageF2V = ones(numOfEdge, numOfLabel);  % message from variable nodes to factor nodes
messageV2F = ones(numOfEdge, numOfLabel);  % message from factor nodes to variable
belief     = ones(n, numOfLabel);          % belief of each variable node.
epsVal = 10E-5;                            % a small value to make sure there is no zero denominator;


%% update factor node 1 to 7. since they only have one variable node neighbor, they only need to be update once.
for i = 1:n
    for j = i:i  
        edgeID = edgeMap(i,j); % find the edgeID of current neighbor j for factor node i.
        for k = 1: numOfLabel        
            messageF2V(edgeID,k) = fun1(noiseCodeword(i), LabelValue(k), p);
        end;        
    end;
end;

%% start bp iteration
for iter = 1:max_iter

    %% update variable nodes 1 to 7
    for j = 1:n
        prodVal = messageF2V(edgeMap(j,j), :);  % get the message from Left side factor node;
        edgeIDArray = full(edgeMap(find(edgeMap(j+1:n+m,j)~= 0) + j, j));  % find neighbor's edgeID;
        
        % mutliply all the incoming message from its neighbor
        for i = 1:length(edgeIDArray)
            edgeID = edgeIDArray(i); % find the index of current neighbor i.    
            prodVal = prodVal .* messageF2V(edgeID,:);
        end;
        
        %send out the message, divide prodVal by the incoming message you
        %get the sending out message;
        for i = 1:length(edgeIDArray)
            edgeID = edgeIDArray(i); % find the index of current neighbor i.                   
            messageV2F(edgeID,:) = prodVal ./ messageF2V(edgeID,:);
            % normalization
            messageV2F(edgeID,:) = messageV2F(edgeID,:) / sum(messageV2F(edgeID,:));
        end;
        
    end;    
    
    %% update factor nodes  8 to 10
    for i = n+1 : n+m
        edgeIDArray = full(edgeMap(i,find(edgeMap(i, :)~= 0)));  % find neighbor's edgeID;
        for j = 1: length(edgeIDArray)
            edgeID = edgeIDArray(j); % find the index of current neighbor j.    
            for k = 1:numOfLabel
                messageF2V(edgeID,k) = sumMsg(LabelValue(k), edgeIDArray, messageV2F, j);
            end;
            % normalization
            messageF2V(edgeID,:) = messageF2V(edgeID,:) / sum(messageF2V(edgeID,:));
        end;
    end;
    
    %% calculate belief
    for j = 1:n
        belief(j, :) = messageF2V(edgeMap(j,j), :);  % get the message from Left side factor node;
        edgeIDArray = full(edgeMap(find(edgeMap(j+1:n+m,j)~= 0) + j, j));  % find neighbor's edgeID;
        % mutliply all the incoming message
        for i = 1 : length(edgeIDArray)
            edgeID = edgeIDArray(i); % find the edgeID of current neighbor i for variable node j.    
            belief(j, :) = belief(j, :) .* messageF2V(edgeID,:);
        end;
    end;    
       
    %% get decoded codeword by belief
    decodeword = fix(belief(:,2) ./ belief(:,1))';
    decodeword(find(decodeword ~= 0)) = 1;
    %% validate the decoded codeword
    checkVal = mod(sum(pchkArray * decodeword'),2);    
    if(iter >= max_iter || checkVal == 0)
        break;
    end;   

end;

if(checkVal ==0)
    disp('decode successful');
else
    disp('decode failed;')
    
end;
disp(['  input codeword: ', int2str(noiseCodeword)]);
disp(['decoded codeword: ', int2str(decodeword)]);