LDPC码编译码仿真

function [newH,g,invPhi] = preprocess(H) % preprocessing of LDPC encode performed only once for every Parity-Check matrix % Steps: % 1. Greedy Algorithm: approximate lower triangulation % 2. Gaussian Elimination: get submatrices A B T C D E % 3. Exchange Columns: make submatrix D~(notated by D) invertible % ChenHz 2013-12-3 [rows,cols] = size(H); newH = double(H); %% Step I: Greedy Algorithm % Initialization % search for the col which has the minimum column weight in H colweight = sum(newH); [mincolweight,c] = min(colweight); % permute cols to move col"c" to the right tmp = newH(:,cols); newH(:,cols) = newH(:,c); newH(:,c) = tmp; % permute rows to move 1s down in the current col rr = find(newH(:,cols)==1);% row indices of 1-elements in col"c" i=0; for ii=1:mincolweight if (rows-ii+1)==rr(mincolweight-i)% check if there is a 1 in the current position i = i+1; else tmp = newH(rows-ii+1,:); newH(rows-ii+1,:) = newH(rr(ii-i),:); newH(rr(ii-i),:) = tmp; end end rerows = mincolweight; recols = 1; g = rerows-1; % Loop for diagonal extension for i=1:cols % search for the col which has the minimum column weight in the cutted H clear colweight colweight = sum(newH(1:rows-rerows,1:cols-recols),1); colweight(colweight==0) = inf;% avert min(sum) getting all-zero cols [mincolweight,c] = min(colweight); % permute cols to move col"c" to the right tmp = newH(:,cols-i); newH(:,cols-i) = newH(:,c); newH(:,c) = tmp; % move 1s to diagonal position and the bottom of newH rr = find(newH(1:rows-rerows,cols-i)==1);% row indices of 1s if ~newH(rows-rerows,cols-i)% if the diagonal position is 0, move a 1 to diagonal position tmp = newH(rr(1),:); newH(rr(1),:) = newH(rows-rerows,:); newH(rows-rerows,:) = tmp; end for ii=2:mincolweight% move other 1s to the bottom of newH tmp = newH(rr(ii),:); newH(rr(ii):(rows-1),:) = newH((rr(ii)+1):rows,:);% all rows below rr(ii) move up 1 row newH(rows,:) = tmp;% move the current 1 to the last row rr = rr-1;% since all rows below rr(ii) move up, row indices of 1s substract 1 end rerows = rerows + mincolweight; recols = recols + 1; g = g+(mincolweight-1); % If ALT is formed, check and terminate the loop if rerows==rows for ii=1:(rows-g) err = sum(newH(ii,cols-g+ii:end),2); if err~=0 error('Matrix H cannot be approximate lower triangularized.'); end end fprintf('Approximate Lower Triangulation is finished.\n'); break end end %% Step II: Gaussian Elimination % top to down, right to left Htmp = newH;% newH doesn't change, this step only gets D for i=(rows-g+1):rows % for every row, ¡ý for j=cols:-1:(cols-(rows-g)+1)% for every elemnt in current row,¡û if Htmp(i,j)==1 Htmp(i,:) = mod(Htmp(i,:)+Htmp((rows-g)-(cols-j),:),2); end end end fprintf('Gaussian Elimination is finished.\n'); %% Step III: Exchange Columns & Compute invPhi % When D is singular, remove this singularity by exchanging columns in D with cols in A&C D = Htmp((rows-g+1):end,(cols-(rows-g)-g+1):(cols-(rows-g))); if rank(gf(D))==g Phi = D; fprintf('Submatrix D~ is invertible.\n'); else % when D is singular % find dependent cols echelon = gf2rref(D);% Reduced row echelon form singularset = ones(1,g);% dependent cols(1-dependent,0-independent) for i=1:g singularset(find(echelon(i,:),1,'first')) = 0;% the first nonzero entry of each row¡úindependent end depcols = find(singularset==1);% dependent cols in D indcols = find(singularset==0);% independent cols in D % find cols in submatix C for exchanging indD = D(:,indcols); % independent submatrix of D colsforex = zeros(1,sum(singularset,2)); for i=1:length(depcols)% for each dependent col for ii=1:cols-(rows-g)-g% replace the dependent col with cols in submatrix C Dtmp = [indD,Htmp((rows-g+1):end,ii)]; if rank(Dtmp)==rank(indD)+1 indD = Dtmp;% if Dtmp is invertible, flag col-ii colsforex(i) = ii; break end end end if rank(gf(indD))~=g error('Submatrix D~ cannot be invertible.'); end % Exchange cols in Htmp to get Phi for i=1:length(depcols) tmp = Htmp(:,colsforex(i)); Htmp(:,colsforex(i)) = Htmp(:,cols-rows+depcols(i)); Htmp(:,cols-rows+depcols(i)) = tmp; end Phi = Htmp((rows-g+1):end,(cols-(rows-g)-g+1):(cols-(rows-g))); % Exchange cols in newH for i=1:length(depcols) tmp = newH(:,colsforex(i)); newH(:,colsforex(i)) = newH(:,cols-rows+depcols(i)); newH(:,cols-rows+depcols(i)) = tmp; end fprintf('After exchanging columns, submatrix D~ is invertible.\n'); end % compute invPhi by Gaussian Elimination invPhi = gf2inv(Phi);