% The program can count the number of the girth 4 of LDPC codes.
% Copyright (C) Yang XIAO, Jun FAN, Beijing Jiaotong University, Aug.8,2007
% E-Mail: yxiao@bjtu.edu.cn.
% 
% Ref.[1] shows that it is necessary for the LDPC codes without girth 4.
% The theorem about girth 4 in Ref.[1] can test the existence of girth 4, 
% but it can not count the number of girth 4. 
% To solve the problem, Ref.[2] provided an algorithm of counting the number of girth 4.
% This program is following the algorithm of Ref [2], and counts the number
% of girth 4 for Gallager LDPC codes.
% 
% Ref:
% [1] Y. Xiao, M -H Lee,Low complexity MIMO-LDPC CDMA systems over multipath channels, 
% IEICE Transactions on Communications, v E89-B, n 5, May, 2006, p 1713-1717
% [2] J. Fan, Y. Xiao, A method of counting the number of cycles in LDPC codes,
% 8th International Conference on Signal Processing, ICSP 2006,Volume: 3,  
% ISBN: 0-7803-9737-1, Digital Object Identifier: 10.1109/ICOSP.2006.345906
% [3] MOHAMMED ALKHANBASH, Gallager LDPC, Matlab center, File Exchange, 2007
% 
% The papers [1] and [2] can be downloaded from Web site of IEICE and IEEE Explore,respectively.

clear;

% Example: To obtain a H, I revised the program of Ref [3].  
% The code generated by the program of Ref [3] not good, we only use it to demo 
% the counting procedure for girth 4.
% If you have a parity check matrix H, you need not use the sub-program.
%
n=1200                         %Code length
c=n/2                          %Length of parity check bits 
wc=3                           %wc: the number of ones in each columns
wr=(n.*wc)./c                  %wr:the number of ones in each row
h=zeros(c,n);                   %generate the empty matrix and start 
                                %assigning the 1s 
%Step 1:LDPC Code Generatio
    onesvector=ones(1,wr);
    start=1;
    finish=wr;   

    for ii = 1:c./wc

       h(ii,start:finish)=1;        
       start=start+wr;              
       finish=(ii.*wr)+wr;          
    end

  for i=1:wc-1
    r=1;

       for ii=1:n
       col_index = (round(rand(1) * (n-1) ))+1;
       randomCol=h(1:fix(c/wc),col_index);
            
       h((i*fix(c/wc))+1   :  ((i*(c/wc))+1  +  fix(c/wc)-1),ii )=randomCol;
               
       end
       r=r+1;
  end
  
%%%%%%% Girth4 Test %%%%%
% Here, we only concern whether a parity matrix H has girth4, instead of
% the number of girth4. If you need know the nember of girth 4, please
% refer Ref. [2].
% get the number of H
H=h;
rows=size(h,1);
cols=size(h,2);

% We can test girth 4 by using the matrix O, see [1].

girth_four_num=0;
for i=1:rows-1
    for j=1:cols-1
        if H(i,j)==1
            x1=i;y1=j;
            for j1=j+1:cols  
                if H(i,j1)==1
                   x2=i;y2=j1;
                   for i1=i+1:rows
                        if H(i1,j1)==1
                            x3=i1;y3=j1;
                            if H(i1,j)==1 %the forth point;
                               x4=i1;y4=j;
                               girth_four_num=girth_four_num+1;
                            end
                        end
                   end
                end
            end
        end
    end
end
girth_four_num
fprintf('girth_four_num=%d\n',girth_four_num)