LDPC码的编译码程序

关于LDPC码编码和BP译码算法的matlab仿真！

%function vHat = decodeProbDomain（rx H N0 iteration）
% Probability-domain sum product algorithm LDPC decoder
%
%  rx        : Received signal vector （column vector）
%  H         : LDPC matrix
%  N0        : Noise variance
%  iteration : Number of iteration
%
%  vHat      : Decoded vector （0/1） 
%
%
% Copyright Bagawan S. Nugroho 2007 
% http://bsnugroho.googlepages.com
H=[1 1 1 0 0 0 0 0;0 0 0 1 1 1 0 0;1 0 0 1 0 0 1 0;0 1 0 0 1 0 0 1];
[M N]=size（H）;rate=M/N;N0=1;
en=10^（N0/10）; 
sigma=1/sqrt（2*rate*en）;

c=[1 0 1 0 1 1 1 1];
rx=c+sigma*randn（1N）;


% Prior probabilities
P1 = ones（size（rx））./（1 + exp（-2*rx./（N0/2）））;
P0 = 1 - P1;

% Initialization
K0 = zeros（M N）;
K1 = zeros（M N）;
rji0 = zeros（M N）;
rji1 = zeros（M N）;
qij0 = H.*repmat（P0 M 1）;
qij1 = H.*repmat（P1 M 1）;

iteration=100;
for n = 1:iteration
   fprintf（‘\n‘）;
   fprintf（‘\n‘）;
   fprintf（‘Iteration : %d\n‘ n）;
   
   % ----- Horizontal step -----
   for i = 1:M
      
      % Find non-zeros in the column
      c1 = find（H（i :））;
      
      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）;
      if（j==1） 
          fprintf（‘Q1=%f   ‘Qi1）;
      else
          fprintf（‘%f   ‘Qi1）;
      end
      % Decode Qj        
      if Qi1 > Qi0
         vHat（j） = 1;
      else
         vHat（j） = 0;
      end
         
   end % for j
  fprintf（‘\nc=‘）;
  
       fprintf（‘%d           ‘vHat）;
       
       
   if mod（vHat*H‘2）==0
       break;end
   
end % for n

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件       84992  2010-01-04 16:08  Ldpc 码.doc
     文件        2901  2010-01-04 16:08  decodeProbDomain.m