jiao2.m

此代码为用比特填充法编写的ldpc码的校验矩阵。它的girth为6

H=zeros(300,600);
H(100,1)=1;
H(120,1)=1;
H(150,1)=1;
for j=2:600
 a=rand(1,300);%利用随机数组把校验矩阵中的某一位随机的设为1
 i=find(a==max(a));%把最大的那一位找到
 H(i,j)=1;
 c=H(i,:);%把一行赋给c
 U=H(:,j);%把一列赋给U
 Ua=(find(U>=0));%Ua是U中1的序号，对应U在H中的位置
  p=find(c>0);
 q=find(U>0);
 p=length(p);
 q=length(q);
 m=8;
 n=3;
 while(p>m)
  a=rand(1,300);%利用随机数组把校验矩阵中的某一位随机的设为1
  i=find(a==max(a));%把最大的那一位找到
  H(i,j)=1;
  c=H(i,:);%把一行赋给c
  U=H(:,j);%把一列赋给U
  Ua=(find(U>=0));%Ua是U中1的序号，对应U在H中的位置
  p=find(c)>0;
 end;
 while(p<m&q<n)
     Ua=(find(U>=0)); 
     ss=find(U>0);
     for pp=1:length(ss)
         s=find(H(ss(pp),:)>0);
       for b=1:length(s)
           r=find(H(:,s(b))>0);
           Ua(r)=0;
       end;
     end;
     if(Ua==0)
         disp('error');
         break;
     end;
     kong=find(Ua==0);
     Ua(kong)=[];
     l=length(Ua);
     qq=zeros(1,l);
     for x=1:l
         ds=find(H(Ua(x),:)>0);
         qq(1,x)=length(ds);%qq与Ua大小相同
     end;
     f=find(qq==min(qq));
     F1=Ua(f);%每一循环会进一步去掉由于新增节点产生的girth4
     zds=zeros(1,length(F1));
     if(length(F1)==1)
         H(F1,j)=1;
     else g=length(F1);
         for b=1:g%对于每一个节点
             sz=find(H(F1(b),:)>0);%每个节点那一行中1的位置
             if(length(sz)==0)
                 zds(1,b)=0;
             else
             for t=1:length(sz)
             y=find(H(:,sz(t)));%y是每一列中1的位置
             bb=zeros(1,length(y));
             for u=1:length(y)
             yy=find(H(y(u),:));%每一行中1的位置
             bb(1,u)=length(yy);
             end;%
             tt=zeros(1,length(sz));
             tt(1,t)=sum(bb);%第一个节点的第一列完成
             end;
             zds(1,b)=sum(tt);%把一个节点的度数和记录下来
             end;
         end;
         da=find(zds==min(zds));
         if(length(da)>1)
             v=rand(1,length(da));
             gg=find(v==min(v));
             H(F1(gg),j)=1;
             U(F1(gg))=1;
         else
         H(F1(da),j)=1;
         U(F1(da))=1;
         end;
     end;
   q=q+1;
 
             
 end;
end;


bl=0;
ks=1;
mm=1;
for j=1:599
    for i=ks:600-j
        a=find(H(:,i)>0);
        b=find(H(:,i+1)>0);
         if(a(1,1)>b(1,1))
            zj=zeros(1,300);
            zj=H(:,i);
            H(:,i)=H(:,i+1);
            H(:,i+1)=zj;
            else if a(1,1)==b(1,1)
                zj=H(:,i);
                H(:,i)=H(:,mm);
                H(:,mm)=zj;
                bl=bl+1;
                mm=mm+1;
                end;
         end;
    end;
    ks=ks+bl;
    bl=0;
end;
                      

[v,vv]=size(H);
flag=0;
ori=0;
while flag==0
temp=H([1:ori+1],[vv-ori:vv]);
vect=diag(temp);
if vect==0
    ori=ori+1;
else flag=1;
temp=H([1:ori+1],[vv-ori:vv]);    
end 
end

vect=diag(temp);
dyl=find(vect>0);%对角线大于0
cz=length(dyl);%对角线大于零的个数  令cz=50
%H[vv-ori vv-ori+cz-1]
%H[1 cz]
temp1=H([1:cz],[vv-ori:vv-ori+cz-1]);
jw=1;
for i=vv-ori:vv-ori+cz-1
    tp=H(:,vv-cz+jw);
    H(:,vv-cz+jw)=H(:,i);
    H(:,i)=tp;
    jw=jw+1;
end;

temp2=H([1:cz],[vv-cz+1:vv]);





%m-g=len
[len,wid]=size(temp2);
T=H([1:len],[600-len+1:600]);
A=H([1:len],[1:300]);
B=H([1:len],[301:600-len]);
E=H([len+1:300],[600-len+1:600]);
D=H([len+1:300],[301:600-len]);
C=H([len+1:300],[1:300]);


W=inv(T);
%Q=zeros(265,65);
fai1=-E*W*B+D;
fai2=-E*W*A+C;
P=inv(fai1);


%EP=-inv(TM).*(-E.*inv(T).*A+C)