genh2.m

包含五种LDPC码编码生成矩阵的编码算法及其比较

function [H,newH,Gp1,rearranged_cols,InvT]=genH2(rows,cols);
%[H,newH,Gp1,rearranged_cols,InvT]=genH2(rows,cols);

%rows需大于10
%对应§2.5.7 一种消除4环且可近似线性编码的随机构造法

parity_check=zeros(rows,cols);

% Set Random number generators initial state
% reset random number generators based on current clock value
rand('state',sum(100*clock));
randn('state',sum(100*clock));

%固定列重为3
bits_per_col=3;

%计算每行1的最多个数
max_ones_per_row=ceil(cols*bits_per_col/rows);

triangle_len=rows-5;
leftcols=cols-triangle_len;%=cols-rows+5

%先在右边的(rows-5)列构造图2.13所示的子矩阵
for i=1:triangle_len
    parity_check([i i+2 i+5],leftcols+i)=1;
end

row_flag=[1 1 2 2 2 repmat(3,1,rows-10) 2 2 1 1 1];

%在左边的(cols-rows+5)列每列随机产生3个1即列重为3
for i=1:leftcols
    a=randperm(rows);
    parity_check(a([1:bits_per_col]),i)=1;
    row_flag(a([1:bits_per_col]))=row_flag(a([1:bits_per_col]))+1;
end

%try to distribute the ones so that the number of ones per row is as uniform as possible
%尝试在列上分散1的位置，使得每行1的个数均衡（相近或相一致）
for i=1:rows
   j=1;
   a=randperm(leftcols);
   while row_flag(i)>max_ones_per_row;  %如果该行行重大于允许的最大行重，则进行处理
      if parity_check(i,a(j))==1 %在左边的(cols-rows+5)列随机选择某一该行上为1的列来处理，将该列该行上的1分散到其他的行
         %随机查找该列上适合放置1（行重小于允许的最大行重，且该位置上为0）的行
         newrow=unidrnd(rows);
         k=0;
         while (row_flag(newrow)>=max_ones_per_row | parity_check(newrow,a(j))==1) & k<rows
            newrow=unidrnd(rows);
            k=k+1;
         end
         if parity_check(newrow,a(j))==0
            %将待处理行上的1转放到找到的行上
            parity_check(newrow,a(j))=1;
            row_flag(newrow)=row_flag(newrow)+1;
            parity_check(i,a(j))=0;
            row_flag(i)=row_flag(i)-1;
         end
      end%if test
      j=j+1;
   end%while loop
end%for loop

%try to eliminate cycles of length 4 in the factor graph
%尝试删除4环
for loop=1:20
   chkfinish=1;
   for r=1:rows
      ones_position=find(parity_check(r,:)==1);
      ones_count=length(ones_position);
      for i=[1:r-1 r+1:rows]
         common=0;
         for j=1:ones_count
            if parity_check(i,ones_position(j))==1
               common=common+1 ;
               if common==1
                  thecol=ones_position(j);
               end
            end
            if common==2
               chkfinish=0; %如果还存在4环，则不结束循环，还进入下一次循环
               common=common-1;
               %如果4环涉及到右边子矩阵的列，则保留右边子矩阵中的列，交换前面的列
               %否则随机决定是保留前面的列还是后面的列
               if ones_position(j)>leftcols | round(rand)==0
                  coltoberearranged=thecol;           %保留后面的列，交换前面的列
                  thecol=ones_position(j);
               else
                  coltoberearranged=ones_position(j); %保留前面的列，交换后面的列
               end
               parity_check(i,coltoberearranged)=3; %make this entry 3 so that we dont use
                                                    %of this entry again while getting rid
                                                    %of other cylces
               row_flag(i)=row_flag(i)-1;
               newrow=unidrnd(rows);
               iteration=0;     %尝试20次在待交换的列中随机查找0
               while (parity_check(newrow,coltoberearranged)~=0 | row_flag(newrow)>max_ones_per_row) & iteration<20
                  newrow=unidrnd(rows);
                  iteration=iteration+1;
               end
               if iteration>=20 %超过20次后则扩大范围随机查找非1的0或3，直到找到为止
                  while parity_check(newrow,coltoberearranged)==1 | row_flag(newrow)>max_ones_per_row | newrow==r | newrow==i
                     newrow=unidrnd(rows);
                  end
               end
               %把该列中找到的0或3置为1
               parity_check(newrow,coltoberearranged)=1;
               row_flag(newrow)=row_flag(newrow)+1;
            end%if common==2
         end%for j=1:ones_count
      end%for i=[1:r-1 r+1:rows]
   end%for r=1:rows

   %如果本次循环已不存在4环，则结束循环，不进入下一次循环
   if chkfinish
      break
   end
end%for loop=1:20

%replace the 3's with 0's
parity_check=parity_check==1;

%Get the Parity Checks
newH=double(parity_check);

%保证φ = -F·inv(T)·B+D可逆
T=newH(1:triangle_len, cols-triangle_len+1:cols);
InvT=qinv_GF2(T);
F=newH(triangle_len+1:rows, cols-triangle_len+1:cols);
FxInvT=mul_GF2(F, InvT);

Acols=cols-rows;

i=1;
B=newH(1:triangle_len, Acols+1:cols-triangle_len);
D=newH(triangle_len+1:rows, Acols+1:cols-triangle_len);
FxInvTxB_D=add_GF2(mul_GF2(FxInvT, B), D);
[inv_FxInvTxB_D stopcol]=qinv_GF2(FxInvTxB_D);

while stopcol~=0 & i<=Acols
   %交换newH的Acols+stopcol和i列
   temp=newH(:,Acols+stopcol);
   newH(:,Acols+stopcol)=newH(:,i);
   newH(:,i)=temp;

   B=newH(1:triangle_len, Acols+1:cols-triangle_len);
   D=newH(triangle_len+1:rows, Acols+1:cols-triangle_len);
   FxInvTxB_D=add_GF2(mul_GF2(FxInvT, B), D);
   [inv_FxInvTxB_D stopcol]=qinv_GF2(FxInvTxB_D);
   i=i+1;
end

if stopcol~=0 & i>Acols
   error('不能使φ = -F·inv(T)·B+D可逆');
end

A=newH(1:triangle_len, 1:Acols);
C=newH(triangle_len+1:rows, 1:Acols);
FxInvTxA_C=add_GF2(mul_GF2(FxInvT, A), C);

Gp1=mul_GF2(inv_FxInvTxB_D, FxInvTxA_C);

%随机重排各列，即随机分散信息比特和校验比特的位置
rearranged_cols=randperm(cols);

H(:, rearranged_cols)=newH;

InvT=InvT==1;
InvT=double(InvT);

Gp1=Gp1==1;
Gp1=double(Gp1);

newH=newH==1;
newH=double(newH);
newH=sparse(newH);

H=H==1;
H=double(H);
H=sparse(H);

%%%%%下面的求方差仅用作评估%%%%
%%计算列重
%col_flag(1:cols)=0;
%for j=1:cols
%   ind=find(H(:,j)==1);
%   col_flag(j)=length(ind);
%end
%%计算行重
%row_flag(1:rows)=0;
%for i=1:rows
%   ind=find(H(i,:)==1);
%   row_flag(i)=length(ind);
%end
%%每行1的个数的方差
%variance=var(row_flag);