decode.m

ldpc系统建模

function [x_hat, success, k] = decode(f,H,maxiter)
% 以下程序针对二进制情况进行处理
% 本程序方法可参照MacMay1999年发表于IEEE TRANSACTIONS.VOL.45.NO.2的
% Good Error-Correcting Codes Based on Very Sparse Matrices

% We assume G is systematic G=[A|I] and G*H'=0 over GFq
% Binary case
%         sigma = 1;                          % AWGN noise deviation
%         x = (sign(randn(1,size(G,1)))+1)/2; % random bits
%         y = mod(x*G,2);                     % encoding 
%         z = 2*y-1;                          % 经过BPSK调制，二进制中的"1"对应于+1，"0"对应于-1。
%         z=z + sigma*randn(1,size(G,2));     % AWGN transmission
% % use the AWGN function directly to control the SNR
%
%         f1=1./(1+exp(-2*z/sigma^2));        % likelihoods
%         f1 = (f1(:))';                      % make it a row vector
%         f0=1-f1;
%         [z_hat, success, k] = ldpc_decode([f0;f1],H,2);
%         x_hat = z_hat(size(G,2)+1-size(G,1):size(G,2));
%         x_hat = x_hat'; 
%
% binary case first, just use the old code
   
   [m,n] = size(H); if m>n, H=H'; [m,n] = size(H); end
   if ~issparse(H) % make H sparse if it is not sparse yet
      [ii,jj,sH] = find(H);%just use
      H = sparse(ii,jj,sH,m,n);%  H=sparse(H)  to substitute for these two programes
   end                         % 但这里这样做是为了下面程序运用ii,jj等变量的方便
   
   f0 = f(1,:); % prob of 0
   f1 = f(2,:);% f1、f0都是从函数获得的参数f中提取的
   
   %initialization
   [ii,jj,sH] = find(H);          % subscript index to nonzero elements of H 这一句可以省略
   indx = sub2ind(size(H),ii,jj); % linear index to nonzero elements of H,将H中的所有1由全下标变单下标
   q0 = H * spdiags(f0(:),0,n,n); % q0的每一列是H的每一列乘上f0里的对应值，所以q0同一列里面的值是相同的等于f0里的对应值
   sq0 = full(q0(indx));          % sq0是H中为1的位置的初始概率的排列,是列向量
   sff0 = sq0;                    % 为后面的计算作准备

   q1 = H * spdiags(f1(:),0,n,n); 
   sq1 = full(q1(indx));
   sff1 = sq1;

   %iterations
   k=0;
   success = 0;
   max_iter = maxiter;   %最大迭代次数
   while ((success == 0) & (k < max_iter)),
      k = k+1;
   
      %horizontal step
      sdq = sq0 - sq1;               % deltaq
      sdq(find(sdq==0)) = 1e-20;     % if   f0 = f1 = .5
      dq = sparse(ii,jj,sdq,m,n);                     % 如果没有f0＝f1的情况，dq可直接由q0－q1得到。
                                                      % 但考虑上述情况，这种方法更方便
      Pdq_v = full(real(exp(sum(spfun('log',dq),2))));% this is ugly but works :)
                                                      % 将每行上不为0的概率差值相乘，Pdq_v应是一列向量
      Pdq = spdiags(Pdq_v(:),0,m,m) * H;              % 将概率差值之积分配到H的每一个1的位置，这时同一行的数值是一样的，等于对应的Pdq_v值
      sPdq = full(Pdq(indx));
      sr0 = (1+sPdq./sdq)./2;                         % 除以sdq就可排除当前位置的概率差值
      sr0(find(abs(sr0) < 1e-20)) = 1e-20;
      sr1 = (1-sPdq./sdq)./2; 
      sr1(find(abs(sr1) < 1e-20)) = 1e-20;
      r0 = sparse(ii,jj,sr0,m,n);
      r1 = sparse(ii,jj,sr1,m,n);
   
      %vertical step
      Pr0_v = full(real(exp(sum(spfun('log',r0),1))));% 每列不为零的概率相乘
      Pr0 = H * spdiags(Pr0_v(:),0,n,n);              % 将概率分配到H的每一个1的位置，这时每列的数值都是一样的
      sPr0 = full(Pr0(indx));
      Q0 = full(sum(sparse(ii,jj,sPr0.*sff0,m,n),1))';% 伪后验概率,用于估计码字
      sq0 = sPr0.*sff0./sr0;                          % 作为下一轮循环变量节点为0的概率
   
      Pr1_v = full(real(exp(sum(spfun('log',r1),1))));
      Pr1 = H * spdiags(Pr1_v(:),0,n,n);
      sPr1 = full(Pr1(indx)); 
      Q1 = full(sum(sparse(ii,jj,sPr1.*sff1,m,n),1))';
      sq1 = sPr1.*sff1./sr1;
   
      sqq = sq0+sq1;
      sq0 = sq0./sqq;
      sq1 = sq1./sqq; % 为了保证sq0+sq1=1
   
      %tentative decoding
      QQ = Q0+Q1;
      Q0 = Q0./QQ;
      Q1 = Q1./QQ;
   
      x_hat = ((sign(Q1-Q0)+1)/2)';
      if rem(H*x_hat',2) == 0, success = 1; end
   end
   % end of binary case