% Funtion to generate an irregular Gallager code with the specified check
% node and variable node degree distributions.Cycles of Length 4 are
% eliminated if possible.
%
%   Prasanna Sethuraman, March 2004.
%
% Example for an irregular code: 
% H=makeLDPC(2000*3.96/8, 2000, [8;1], [2, 4, 18; 0.5100, 0.4200,0.0700]);
%
% Checks:
% dv = [2, 0.5100;  4, 0.4200; 18,0.0700];
% dc = 8; 
% The equation: (2*0.5100) + (4*0.4200) + (18*0.0700) = 3.96
% and           (8*1) = 8
% and so, rate = 3.96/8
%
% Note: The resulting parity check matrix may not necessarily have the
% exact distributions specified by dc and dv. However, the actual
% distributions will be made as close as possible. Check the actual
% distributions with the 'hist' function.


 
% M is row dimenstion and N is coloumn dimenstion of H matrix.
% dc and dv are respectively check node and variable node degree
% distributions.
function [H] = makeLDPC(M, N, dc, dv)
 
% some preliminary checks and computations
cndist_sum1 = 0;
cndist_sum2 = 0;
for i=1:length(dc(1,:))
    cndist(1,i) = dc(1,i);
    cndist(2,i) = floor(dc(2,i)*M);
    cndist_sum1 = cndist_sum1 + dc(2,i);
    cndist_sum2 = cndist_sum2 + dc(1,i)*dc(2,i)*M;
end;
 
vndist_sum1 = 0;
vndist_sum2 = 0;
for i=1:length(dv(1,:))
    vndist(1,i) = dv(1,i);
    vndist(2,i) = floor(dv(2,i)*N);
    vndist_sum1 = vndist_sum1 + dv(2,i);
    vndist_sum2 = vndist_sum2 + dv(1,i)*dv(2,i)*N;
end;
 
if cndist_sum1 ~= 1
    error('Fraction of check node degrees does not form a distribution');
elseif vndist_sum1 ~= 1
    error('Fraction of variable node degrees does not form a distribution');
elseif cndist_sum2 ~= vndist_sum2
    error('Not a valid distribution pair. Equation not satisfied');
end;
 
H = sparse(M,N);
 
rowCount = zeros(1,M);
rowDegCount = sparse(1,N);
 
cycle =  false;
 
n=1;                % a local index variable for vndist
m=1;                % a local index variable for cndist
lsum=vndist(2,n);   % a local summation variable
for j=1:N
    if j<=lsum
        ;   
    else
        n = n+1;
        if n>length(vndist(1,:)), break; end;
        lsum = lsum + vndist(2,n);
    end;
    if rowDegCount(cndist(1,m))<=cndist(2,m)
        ;
    else
        m=m+1;
        if m>length(cndist(1,:)), break; end;
    end;
   
    for k=1:vndist(1,n)
        r=randperm(M);
        for lp=1:2*M
            if lp==M+1
                cycle=true;
            else
                i=r(mod(lp,M+1));
            end;
            if H(i,j)==0
                if rowCount(i)<cndist(1,m)
                    if ( no4cycles(H,i,j)==false & cycle==false) %& tooLong<floor(M/2)
                        continue;
                    else
                        H(i,j)=1;
                        rowCount(i)=rowCount(i)+1;
                        rowDegCount(rowCount(i))=rowDegCount(rowCount(i))+1;
                        cycle=false;
                        break;
                    end;
                end;
            end;
        end;
    end;
    % ---- debug -----
    colCount(j)=length(find(H(:,j)));
    % ---- debug -----
end;
disp('done');