function vHat = decodeBitFlipping(rx, H, iteration)
%  rx        : Received signal vector (column vector)
%  H         : LDPC matrix
%  iteration : Number of iteration
%  vHat      : Decoded vector (0/1) 
[M N] = size(H);
% Prior hard-decision
ci = 0.5*(sign(rx') + 1);
% Initialization
rji = zeros(M, N);
% Asscociate the ci matrix with non-zero elements of H
qij = H.*repmat(ci, M, 1);
for n = 1:iteration
   % ----- Horizontal step -----
   for i = 1:M
      c1 = find(H(i, :)); 
      for k = 1:length(c1)
         rji(i, c1(k)) = mod(sum(qij(i, c1)) + qij(i, c1(k)), 2);        
      end
   end   
   % ------ Vertical step ------
   for j = 1:N
      r1 = find(H(:, j));
      numOfOnes = length(find(rji(r1, j)));
      for k = 1:length(r1)        
         % Update qij, set '1' for majority of 1s else '0', excluding r1(k)
         if numOfOnes + ci(j) >= length(r1) - numOfOnes + rji(r1(k), j)
            qij(r1(k), j) = 1;
         else
            qij(r1(k), j) = 0;
         end      
      end
      if numOfOnes + ci(j) >= length(r1) - numOfOnes
         vHat(j) = 1;
      else
         vHat(j) = 0;
      end       
   end 
end