ldpccodeanddecode.m

matlab环境下 LDPC码的编译码算法仿真

%************ preparation ******************************
M=100;
N=200;
iter = 5;
onePerCol=3;
sum=0;

%*********** creat H matrix *****************************
%H=makeLdpc(M, N, 1, 1, onePerCol);
if N/M ~= 2
   fprintf('Code rate must be 1/2\n');
end
onePerRow = (N/M)*onePerCol;
fprintf('Creating LDPC matrix...\n');

%case {1}
% Distribute 1s uniformly at random within column
      for i = 1:N
          onesInCol=zeros(M,N);
          onesInCol(: , i)= randperm(M)';
      end
        
      % Create non zero elements (1s) index
      r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
      tmp = repmat([1:N], onePerCol, 1);
      c = reshape(tmp, N*onePerCol, 1);
     
      % Make the number of 1s between rows as uniform as possible     
      
      % Order row index
      [r, ix] = sort(r);%升续排列这60个值
      
      % Order column index based on row index
      for i = 1:N*onePerCol
         cSort(i, :) = c(ix(i));
      end
      % Create new row index with uniform weight
      tmp = repmat([1:M], onePerRow, 1);
      r = reshape(tmp, N*onePerCol, 1);
      
      % Create sparse matrix H
      % Remove any duplicate non zero elements index using logical AND
      S = and(sparse(r, cSort, 1, M, N), ones(M, N));
      H = full(S); 
      
%*********** code ***************************************
seldata=round(rand(1,M)); 
checkcode=seldata*H;
for i=1:M
    checkcode(1,i)=mod(checkcode(1,i),2);
end
code=zeros(1,N);
code(1,M+1:N)=checkcode(1,1:M);
code(1,1:M)=seldata;               %1*N
rx=code;

%*********** decode **************************************
%vHat=decodeProbDomain(rx,H,iter);
[M N] = size(H);
% Prior probabilities
P1 = ones(size(rx))./(1 + exp(-2*rx));
P0 = 1 - P1;
K0 = zeros(M, N);
K1 = zeros(M, N);
rji0 = zeros(M, N);
rji1 = zeros(M, N);
qij0 = H*repmat(P0', 1, N);%M应改成N
qij1 = H*repmat(P1', 1, N);

   % ----- Horizontal step -----
   for i = 1:M
      
      % Find non-zeros in the column
      c1 = find(H(i, :));%找到非零，即1的位置
      
      for k = 1:length(c1)
         
         % Get column products of drji\c1(l)
         drji = 1;
         for l = 1:length(c1)
            if l~= k
               drji = drji*(qij0(i, c1(l)) - qij1(i, c1(l)));
            end
         end % for l
         
         rji0(i, c1(k)) = (1 + drji)/2;
         rji1(i, c1(k)) = (1 - drji)/2;
         
      end % for k
      
   end % for i
   % ------ Vertical step ------
   for j = 1:N
      
      % Find non-zeros in the row
      r1 = find(H(:, j));
      
      for k = 1:length(r1)
        
         % Get row products of prodOfrij\ri(l)
         prodOfrij0 = 1;
         prodOfrij1 = 1;   
         for l = 1:length(r1)
            if l~= k
               prodOfrij0 = prodOfrij0*rji0(r1(l), j);
               prodOfrij1 = prodOfrij1*rji1(r1(l), j);
            end
         end % for l
         
         % Update constants
         K0(r1(k), j) = P0(j)*prodOfrij0;
         K1(r1(k), j) = P1(j)*prodOfrij1;
         
         % Update qij0 and qij1
         qij0(r1(k), j) = K0(r1(k), j)./(K0(r1(k), j) + K1(r1(k), j));
         qij1(r1(k), j) = K1(r1(k), j)./(K0(r1(k), j) + K1(r1(k), j));
               
      end % for k
      % Update constants
      Ki0 = P0(j)*prod(rji0(r1, j));
      Ki1 = P1(j)*prod(rji1(r1, j));
      
      % Get Qj
      Qi0 = Ki0/(Ki0 + Ki1);
      Qi1 = Ki1/(Ki0 + Ki1);
      
      % Decode Qj        
      if Qi1 > Qi0
         vHat(j) = 1;
      else
         vHat(j) = 0;
      end
         
   end % for j
for i=1:M
chazhi=abs(vHat(1,i)-seldata(1,i));
sum=sum+chazhi;
end