% Joshua R. Smith 11/20/11
% Sum-Product / Belief propagation iterative decoding of the
% 7,4 hamming code
% Syndrome decoding point of view
%   ==> Assume we know the codeword that was sent was all zeros
f = .05; % noise level...probability of 
pn1 = f;
pn0 = 1-f;


s = [1 0 1 0]'; % src bits
P = [[1 1 1 0]; [0 1 1 1]; [1 0 1 1];];
gt = [eye(4); P]; % Transpose of generator matrix
h = [P eye(3)]; % Parity check matrix
[m n] = size(h); % m: number of checks  n: number of variable bits
q0 = pn0*h; % prob that bit n has value 0, given check bits other than m
q1 = pn1*h;


t = mod(gt * s,2); % transmitted vector

nz = [0 1 0 0 0 0 0]'; % noise vector
                       % In the syndrome decoding perspective used here,
                       % the whole goal is to to estimate this noise
                       % vector.  The noise vector is everything the 
                       % receiver needs in order to successfuly decode, 
                       % since the RCV is given r, and knowledge of 
                       % r and n gives t and thus [trivially] s
r = mod(t+nz,2); % received vector (transmitted vector corrupted by noise)

z = mod(h*r,2) % syndrome: observed list of parity check results
bi2de(z')      % for the 7,4 Hamming code, there are 3 parity checks
               % each 1 in the syndrome is a parity check failure
q1

q0+q1
ITMAX = 100; % number of iterations to do
pp1vect = zeros(ITMAX,n); % allocate a big array to hold the time series of bit probabilities (for plotting)
pp0vect = zeros(ITMAX,n); % for any value of j, it should be the case that pp0vect(j,:) + pp1vect(j,:) = [1.0 1.0 ... 1.0] i.e. n 1s 
for t=1:ITMAX,

% Horizontal step
r0 = zeros(m,n); % prob that check bit m is satisfied (syndrome bit z(m)=0) if bit n is fixed at 0
r1 = zeros(m,n); % prob that check m is satisfied (syndrome bit m=0) if bit n is fixed at 1
for i=1:m, % for each check, i.e. for each row in parity check matrix
    ind = find(h(i,:)); % ind is a vector of indices (col #s) of non-zero entries in parity check matrix
    rweight = length(ind); % row weight: number of 1s in the row

    
    for k = 0:(2^rweight)-1; % Consider every setting of bits that participate in check
        hyp = de2bi(k,rweight); % de2bi=="decimal to binary"...hyp is a hypothetical setting of all bit values that participate in check
        xh=de2bi(sum(hyp.*(2.^(ind-1))),n);  % "unpacked" hypothesis...this maps hypothetical participating bit values into their proper places in a bit vector of length n

        % Precompute the "complete product term" qp, not excluding anything [i.e. no "\ n"]...
        % from this we will generate the actual product terms we want by
        % dividing out the one term we don't want each time through the 2nd b loop
        % below.  Note that this 1st b loop has to be inside the k loop
        % because it uses hyp, which is set by the k loop
        % ...which makes sense, because the "product term" is an estimate
        % of the probability of the hypothesis under consideration
        % (without considering parity checks)
        qp = 1.0;
        for b = 0:(rweight-1);
            jp = ind(b+1); % jp is a pointer into the row of the parity check matrix... Note the +1 is cuz Matlab indices start at 1
            qp = qp*( (1-hyp(b+1))*q0(i,jp) + hyp(b+1)*q1(i,jp)); % This is using hyp as an if statement to decide whether to multiply in q0 or q1
        end
        
        for b = 0:rweight-1; % pick a bit to knock out?
            jp = ind(b+1); % jp is a pointer into the row of the parity check matrix... Note the +1 is cuz Matlab indices start at 1
                            % jp is a pointer to the variable node to
                            % exclude from the present sum
            % NB: *hyp(b+1) is equivalent to "if bit b is 0, then 0"
            qpp = qp/( (1-hyp(b+1))*q0(i,jp) + hyp(b+1)*q1(i,jp)); % Divide out to exclude prob for currently "knocked out" bit
            %keyboard
            % This is coded a bit differently than MacKay's notation, but
            % (we believe!) does the same thing.  We want
            % P(z|x) = 1 when z = parity of x.  So, P(z|x)=1+sum(x)+z.
            % We don't explicitly separate out the xh(jp) bit [the one
            % "left out" in the parity computation...rather, we do that
            % with the prefactor [hyp(b+1) or (1-hyp(b+1))] which decides
            % whether we accumulate into r0 or r1
            % Note that hyp(b+1) = 1 when xh(jp) = 1
            r0(i,jp) = r0(i,jp) + (1-hyp(b+1))*(mod(1+sum(xh)+z(i),2))*qpp; % Coef before qpp verifies that parity check is satisfied
            r1(i,jp) = r1(i,jp) + (hyp(b+1))*(mod(1+sum(xh)+z(i),2))*qpp; % Coef before that verifies that 0 or 1 matches current hyp bit
        end % for b
    end % for k
end % for i

r1

% Vertical
q0 = pn0*h;
q1 = pn1*h;
for j=1:n,
  for i=1:m,
    for k=1:m;
        if ((k~=i)&&(h(k,j)==1)&&h(i,j)==1),
          q0(i,j)=q0(i,j)*r0(k,j);
          q1(i,j)=q1(i,j)*r1(k,j);
      end
    end
  end
end

% Normalize q0 and q1
for j=1:n,
  for i=1:m,
      sm = q0(i,j)+q1(i,j);
      if ((sm > 0)&&(h(i,j)==1)),
          q0(i,j)=q0(i,j)/sm;
          q1(i,j)=q1(i,j)/sm;
      end
  end
end


% Compute pseduo-posterior, i.e. our best estimates
% of bit probabilities
for j=1:n,
  pp0(j) = pn0;
  pp1(j) = pn1;
  for i=1:m,
      if (h(i,j)==1),
        pp0(j) = pp0(j) * r0(i,j);
        pp1(j) = pp1(j) * r1(i,j);
      end
  end
end

for j=1:n,
  for i=1:m,
      sm = pp0(j)+pp1(j);
      if ((sm > 0)&&(h(i,j)==1)),
          pp0(j)=pp0(j)/sm;
          pp1(j)=pp1(j)/sm;
      end
  end
end

q1
pp1
pp0vect(t,:) = pp0;
pp1vect(t,:) = pp1;
end % for t

plot(pp1vect(:,1),'k'); hold on
plot(pp1vect(:,2),'r'); hold on % this one is harcoded to be red
plot(pp1vect(:,3),'k'); hold on
plot(pp1vect(:,4),'k'); hold on
plot(pp1vect(:,5),'k'); hold on
plot(pp1vect(:,6),'k'); hold on
plot(pp1vect(:,7),'k');
hold off
