ldpc_lintimeenc_parameter1.m

用Matlab实现的LDPC编码程序

%H is a uncomleted downtriangle matrix,but very alike;
%gap is a matrix gap between a truly downtriangle matrix and H;
%the function of this procedure can get parameter such as H and inv_phi1 for linetime encode
%But H is slightly modified to get parameter,but the H still has original (good) girth distribution  
%When H is full rank(M),this procedure not only can rearrange the columns of H so as to make phi invertible.
%When H is not full rank(M-1),this procedure can make H full rank(H),also can rearrange the columns of H so as to make phi invertible.
%function [H,inv_phi1,g]=LDPC_LinTimEnc_Parameter(H)
load regularHmatrix_N20_K15_j3_k4.mat%%HforOFDM222.mat%
H=Hsparse;

[m,n]=size(H);flag=0;
gap=max(find(H(:,1)));g=gap-1;
Hinverse=H(m:-1:1,:);H=Hinverse(:,n:-1:1);%uptriangle to downtriangle;
kk=0;T=H(1:m-g,n-m+g+1:n);E=H(m-g+1:m,n-m+g+1:n);
EmulInvT=sparse(g,m-g);EmulInvT=mod(E/T,2);rand('state',34);
%mod(EmulInvT*T-E,2)%check (E*inv(T))*T=I
A=H(1:m-g,1:n-m);B=H(1:m-g,n-m+1:n-m+g);
C=H(m-g+1:m,1:n-m);D=H(m-g+1:m,n-m+1:n-m+g);
CC=mod(-EmulInvT*A+C,2);DD=mod(-EmulInvT*B+D,2);
%*************compute the rank of phi*************
DDcheck=DD;[M,N]=size(DDcheck);Etr=spdiags(ones(M),0,M,M);
swapt=sparse(zeros(1,N));j=1;Recjfori=[];% j 指向台阶
for i=1:N,% i 指向列 
  a=find(DDcheck(j:M,i)~=0)+j-1;
  if ~isempty(a),
     L=length(a);
     if (DDcheck(j,i)==0),
        swapt=DDcheck(j,:);swapE=Etr(j,:);
        DDcheck(j,:)=DDcheck(a(1),:);Etr(j,:)=Etr(a(1),:);
        DDcheck(a(1),:)=swapt;Etr(a(1),:)=swapE;
     end;
     if L~=1,
        DDmat=ones(L-1,N)*spdiags(DDcheck(j,:)',0,N,N);
        DDcheck(a(2:L),:)=mod(DDcheck(a(2:L),:)+DDmat,2);
        Etrmat=ones(L-1,M)*spdiags(Etr(j,:)',0,M,M);
        Etr(a(2:L),:)=mod(Etr(a(2:L),:)+Etrmat,2);
     end;
     j=j+1;
  else
     l=M;flag1=0;
     while (l>=1)&(flag1==0),
        if DDcheck(l,i)==1,
            flag1=1;
        end;
        l=l-1;
     end;
     l=l+1;Recjfori=[i,l];%find linear correlate columns
  end;
end;
mid=any(DDcheck,2);
rank2=sum(mid);row0indx=find(mid==0);CCDD=[CC,DD];
[rankCCDD,maxnocorrcolgroup,diagCCDD,EtrCCDD]=gf2rankimpr(CCDD);
%**********make phi full rank  by add 1 in H****************
if rankCCDD==g,
   disp('A slight process,[CC,DD] is full rank,H succeess'); 
   Hs=H(:,n-m+1:n-m+g);H(:,n-m+1:n-m+g)=H(:,maxnocorrcolgroup);
   H(:,maxnocorrcolgroup)=Hs;
   inv_phi1=mod(diagCCDD(:,maxnocorrcolgroup)\EtrCCDD,2);%
elseif rankCCDD==g-1,
   disp('[CC,DD] is not full rank,H not success,but can be processed');
   DeltMat=sparse(g,g);
   DeltMat(row0indx,Recjfori(1))=1;%或者DeltMat(M-1,M-4)=1
   invEtr=mod(inv(Etr),2);
   if any(mod(Etr*invEtr,2)-spdiags(ones(g,1),0,g,g)~=0), %check function
      error('invE is wrong!');
   end;
   DeltD=invEtr*DeltMat;D1=mod(D+DeltD,2);
   H(m-g+1:m,n-m+1:n-m+g)=D1;
   [indx,indy]=find(DeltD);indx=indx+m-g;indy=indy+n-m;
   phi1=mod(-EmulInvT*B+D1,2);%gf2rank(phi1)==g%check phi1 is full rank
   %%%%%%%%%%%%%%%%get inv_phi1%%%%%%%%%%%%%%%%%
   inv_phi1=sparse(size(phi1));phi1_downdiag=mod(DeltMat+DDcheck,2);
   %%invEtr*phi1_downdiag=phi1;phi1*inv_phi1=I(eye matrix);
   %phi1_downdiag*inv_phi1=Etr;downdiag equ can be solve by matlab;
   inv_phi1=mod(phi1_downdiag\Etr,2),%mod(phi1*inv_phi1,2)-spdiags(ones(g,1),0,M,M)
   %%%%%%%%%%%girth modification by delete one entry in H%%%%%%%%%%%%%%%
   xset=find(H(indx,:));L=length(xset);
   [M,N]=size(H);VnodeDeg=sum(H);gir=2000*ones(1,N);
   for i=1:L,
      x=xset(i);
      U1=find(H(:,x)~=0);
      for y=1:VnodeDeg(x);
         flag=0;g0=2;Vnode1=x;Cnode1=U1(y);
         while (flag==0)&(g0<=20),
           Vnode2=find(any(H(Cnode1,:),1));
           Vnode2(find(ismember(Vnode2,Vnode1)))=[];
           Cnode2=find(any(H(:,Vnode2),2));
           Cnode2(find(ismember(Cnode2,Cnode1)))=[];
           g0=g0+2;
           if any(ismember(U1,Cnode2)),
             flag=1;
           end;  
           Vnode1=Vnode2;Cnode1=Cnode2;
         end;
         if g0<gir(x),
            gir(x)=g0;
         end;
      end;
   end;
   miny=find(g==min(gir));
   miny(find(ismember(miny,indy)))=[];
   if ~isempty(miny),
     H(indx,miny(1))=0;
   end;
   
   save HmodifiedforOFDM H
   %save invmodifield inv_phi1
   %save gmodifield g
   
   
   %gg=girthdistribution(H);%find girth is still high
else
   error('Please regenerate check matrix,this H fail to lintimEnc');
   
   %save invmodifield inv_phi1
   
end;
%mod(E*DD-H,2)