%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function [H_sys, G_sys] = systematic(H);
%% This routine removes any dependent rows from an H parity matrix.
%% I then finds the systematic G and H matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Definition of input variables
%% H      : parity-check matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Definition of output variables
%% H_sys      : Systematic H that corresponds to the systematic G matrix
%% G_sys      : Systematic G matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% EE379, 06/2006, Chien-Hsin Lee ; subsequent edits by J. Cioffi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [H_sys, G_sys] = systematic(H);
[row col] = size(H);
H_in      = H;
offset    = 0;
r = 1;  c = 1;
num_dep_row   = 0;
%% Gaussion Elimination PART1 
%% This find the leading 1 in each column 
%% and eliminate the 1s in the row below
for i = 1:row
    row_order(i) = i;
end
while( r <= row && c <= col)
    row_lst = find(H(r:row,c));
    if( isempty(row_lst))
        c = c + 1;
    else
        if( length(row_lst) > 1)  
            for i = 2:length(row_lst)
                H(row_lst(i)+r-1,:) = mod(H(row_lst(1)+r-1,:) + H(row_lst(i)+r-1,:),2);
            end
        end
        H_temp                   = H;
        H_temp(r,:)              = H(row_lst(1)+r-1,:);
        H_temp(row_lst(1)+r-1,:) = H(r,:);
        H                        = H_temp;
        row_temp                 = row_order;
        row_temp(r)              = row_order(row_lst(1)+r-1);
        row_temp(row_lst(1)+r-1) = row_order(r);
        row_order                = row_temp;
        r = r + 1;
        c = c + 1;
    end
end
r = row;
while( H(r,:) == 0 & r ~= 0)
    H(r,:) = [];
    r = r-1;
end
num_dep_row = row-r;
%flg = sprintf('Report: %d dependent rows are removed;\n  Old matrix is %d by %d;\n  New matrix is %d by %d',row-r, row, col, r, col);
%disp(flg);
row = r;

%% Column Permutation To Put Diaginal In All 1's
for c = 1:col
    column(c) = c;
end 
r = 1; c = 1;
while( r <= row && c <= col)
    if( H(r,c) ~= 1)
        H(:,c:col) = [H(:,c+1:col),H(:,c)];
        column(c:col) = [column(c+1:col),column(c)];
        r = r-1; c= c-1;
    end    
    r = r+1; c=c+1;
end
 
%% Gaussian Elimination Part II To Get The Systematic H
for i = 0:row-1
    %% find backward leading cols and eliminate the upper trangle 1s
    col_lst = find(H(:,row-i));
    if( length(col_lst) > 1)
        for j = 1:length(col_lst)-1
            H(col_lst(j),:)=mod(H(col_lst(j),:)+H(row-i,:),2);
        end
    end
end

%% Get G And H Matrix
G_sys = [H(:,row+1:col)',eye(col-row)];
H_no_dep = [];
if(num_dep_row > 0)
    x = row_order(row+1:row+num_dep_row);
    x = sort(x); 
    j = 1;
	for i = 1:row+num_dep_row
        if(i~=x(j))
            H_no_dep = [H_no_dep;H_in(i,:)];
        else
            j=j+1;
        end
    end
else
    H_no_dep = H_in;
end

for i = 1:col
    H_sys(:,i)=H_no_dep(:,column(i));
end