genh3.m

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

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

%rows需大于10

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
mean_bits_per_col=3;

%预置行重
bits_per_row=ceil(mean_bits_per_col*cols/rows);

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

for i=1:triangle_len
    if round(rand)==1
        set1=[i i+2 i+5];
    else
        set1=[i i+2];
    end
    parity_check(set1,leftcols+i)=1;
end

%计算行重
row_flag=zeros(1,rows);
for i=1:rows
   ind=find(parity_check(i,:)==1);
   row_flag(i)=length(ind);
end

%在左边的(cols-rows+5)列每行随机产生些1使行重为预置行重
for i=1:rows
    a=randperm(leftcols);
    parity_check(i,a([1:bits_per_row-row_flag(i)]))=1;
end

row_flag=repmat(bits_per_row,1,rows);

%计算列重
col_flag=zeros(1,leftcols);
for j=1:leftcols
   ind=find(parity_check(:,j)==1);
   col_flag(j)=length(ind);
end

%尝试在行上分散1的位置，使得每列列重不少于2
%消除列重0，先将其列重变成1
ind0=find(col_flag==0);
for j=ind0
   ind=find(col_flag>2);
   colpos=ind(ceil(length(ind)*rand));
   ind=find(parity_check(:,colpos)==1);
   rowpos=ind(ceil(length(ind)*rand));
   parity_check(rowpos,j)=1;
   col_flag(j)=1;
   parity_check(rowpos,colpos)=0;
   col_flag(colpos)=col_flag(colpos)-1;
end

%消除列重1，将其列重变成2
ind0=find(col_flag==1);
for j=ind0
   ind=find(col_flag>2);
   colpos=ind(ceil(length(ind)*rand));
   ind=find((parity_check(:,colpos)-parity_check(:,j))==1);
   rowpos=ind(ceil(length(ind)*rand));
   parity_check(rowpos,j)=1;
   col_flag(j)=2;
   parity_check(rowpos,colpos)=0;
   col_flag(colpos)=col_flag(colpos)-1;
end

%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)>bits_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)>bits_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);