format long
clear all;

r=2/3
H = dvbs2ldpc®;
%e=0
% spy(H); % Visualize the location of nonzero elements in H.

% Construct a default LDPC encoder object
henc = fec.ldpcenc(H);

for SNRdB = 0

for n=1:100;
% Generate a random binary message
msg = randint(1,henc.NumInfoBits,2);

% Encode the message
codeword = encode(henc, msg);

% Verify the parity checks (which should be a zero vector)
paritychecks = mod(henc.ParityCheckMatrix * codeword', 2);

hdec = fec.ldpcdec(H);

% Construct a BPSK modulator object
modObj = modem.pskmod('M',2,'InputType','Bit');

% Modulate the signal (map bit 0 to 1 + 0i, bit 1 to -1 + 0i)
modulatedsig = modulate(modObj, codeword);

% Noise parameters

sigma = sqrt(10^(-SNRdB/10));

% Transmit signal through AWGN channel
receivedsig = awgn(modulatedsig, SNRdB, 0); % Signal power = 0 dBW

% % Visualize received signal
% scatterplot(receivedsig)

% Construct a BPSK demodulator object to compute
% log-likelihood ratios
demodObj = modem.pskdemod(modObj,'DecisionType','LLR', ...
'NoiseVariance',sigma^2);

% Compute log-likelihood ratios (AWGN channel)
llr = demodulate(demodObj, receivedsig);

% Decode received signal
decodedmsg = decode(hdec, llr);

N_err=nnz(decodedmsg-msg);

% % Actual number of iterations executed
% disp(['Number of iterations executed = ' ...
% num2str(hdec.ActualNumIterations)]);
% % Number of parity-checks violated
% disp(['Number of parity-checks violated = ' ...
% num2str(sum(hdec.FinalParityChecks))]);
% % Compare with original message
% disp(['Number of bits incorrectly decoded = ' ...
% num2str(N_err)])


a(n)=N_err; % Общее количество ошибок

end

BER=sum(a)/(n*henc.NumInfoBits)


end
BER
FER=(nnz(a))/n

function [vHat n] = decodeLogMinSumOffsetSinglePassFpga(rx, H, iteration, beta)

llr_bit_width = 8;

max_val = (2^(llr_bit_width-1))-1;
min_val = -(2^(llr_bit_width-1));

[M N] = size(H);

Lci = rx';

% Initialization
Lrji = zeros(M, N);
Pibetaij = zeros(M, N);

Zn   = Lci;

for n = 1:iteration
  
   Lrji_prev =  Lrji;
  
   betaij  = ( repmat( Zn, M, 1 ).*H - Lrji );
  
   %Limit data by bit width value. Like in fpga.
    for i=1:size(betaij,1)*size(betaij,2)
              
        if (betaij(i) > max_val)
          betaij(i) = max_val;
        end
        
        if (betaij(i) < min_val)
          betaij(i) = min_val;
        end
    end    
  
   %In convertion two's compliment to absolute value
   %didn't add 1 to save some resources
    %for i=1:size(betaij,1)*size(betaij,2)
    %  if ( betaij(i) < 0 )
    %    betaij(i) = betaij(i) + 1;  
    %  end
    %end
  
   betaij  = abs( betaij );
   alphaij = nonzerosign( repmat( Zn, M, 1 ).*H - Lrji );

   % ----- Horizontal and Vertical step -----
   for i = 1:M
      
      % Find non-zeros in the column
      c1 = find(H(i, :));
      
      for k = 1:length(c1)

         % Maximum betaij\c1(k) - beta
         beta_row = betaij(i, c1);
         beta_row(k) = realmax;
         minOfbetaij = min( beta_row ) - beta;
         minOfbetaij = max( [minOfbetaij 0] );
        
         % Multiplication alphaij\c1(k) (use '*' since alphaij are -1/1s)
         prodOfalphaij = prod(alphaij(i, c1))*alphaij(i, c1(k));
        
         % Update L(rji)
         Lrji(i, c1(k)) = prodOfalphaij*( minOfbetaij );
                  
         Zn_inc = Lrji(i, c1(k)) - Lrji_prev(i, c1(k));
        
         if ( Zn_inc  > max_val )
             Zn_inc = max_val;
         end
        
         if ( Zn_inc  < min_val )
             Zn_inc = min_val;
         end
                
         Zn( c1(k) ) = Zn( c1(k) ) + Zn_inc;
        
         if ( Zn( c1(k) ) > max_val )
             Zn( c1(k) ) = max_val;
         end
        
         if ( Zn( c1(k) ) < min_val )
             Zn( c1(k) ) = min_val;
         end
        
      end % for k

   end % for i
  
   %check parity
   if ( n ~= 1 )
     bits = ( nonzerosign( Zn_prev ) - 1 ) / -2;
     if ( sum(mod(H*bits',2)) == 0 )
       break;
     end
   end
  
   Zn_prev = Zn;
end % for n

for i=1:length(Zn)
  if( (Zn(i)) < 0 )
     vHat(i) = 1;
  else
     vHat(i) = 0;
  end
end
%nn = n_last;