ldpcber.m

LDPC encoder and decoder design in MATLAB

% Bit error rate of BPSK modulated LDPC codes under AWGN channel

clc; 
clear all;

% LDPC matrix size, rate must be 1/2
% Warning: encoding - decoding can be very long for large LDPC matrix!
M = input('the no. of rows in the generator matrix:');
N = input('the no. of cols in the generator matrix:');

% Method for creating LDPC matrix (0 = Evencol; 1 = Evenboth)
method = 1;

% Eliminate length-4 cycle
noCycle = 1;

% Number of 1s per column for LDPC matrix
onePerCol = 3;

% LDPC matrix reorder strategy (0 = First; 1 = Mincol; 2 = Minprod)
strategy = 2;

% EbN0 in dB
EbN0 = input('enter the value of EbN0:');

% Number of iteration;
iter = 5;

% Number of frame (N bits per frame)
frame = 1;

% Make the LDPC matrix
H = makeLdpc(M, N, 1, 1, onePerCol);

for i = 1:length(EbN0)
   ber0(i) = 0;
   ber1(i) = 0;
   ber2(i) = 0;
   ber3(i) =0;
   
   % Make random data (0/1)
   dSource = round(rand(M, frame));
   
   for j = 1:frame
      fprintf('Frame : %d\n', j);
      
      % Encoding message
      [c, newH] = makeParityChk(dSource(:, j), H, strategy);
      u = [c; dSource(:, j)];

      % BPSK modulation
      bpskMod = 2*u - 1;
      % Additional white gaussian noise
      N0 = 1/(exp(EbN0(i)*log(10)/10));
      tx = bpskMod + sqrt(N0)*randn(size(bpskMod));

      % Decoding (select decoding method)
      vhat0 = decodeProbDomain(tx, newH, N0, iter);
      vhat1 = decodeLogDomain(tx, newH, N0, iter);
      vhat2 = decodeLogDomainSimple(tx, newH, iter);
      vhat3 = decodeBitFlip(tx, newH, iter);
   
      % Get bit error rate (for brevity, BER calculation includes parity bits)
      [num0, rat0] = biterr(vhat0', u);
      ber0(i) = (ber0(i) + rat0);
      
      [num1, rat1] = biterr(vhat1', u);
      ber1(i) = (ber1(i) + rat1);
      
      [num2, rat2] = biterr(vhat2', u);
      ber2(i) = (ber2(i) + rat2);
     
      [num3, rat3] = biterr(vhat3', u);
      ber3(i) = (ber3(i) + rat3);
      
      
   end % for j
   
   % Get average of BER
   ber0(i) = ber0(i)/frame;
   ber1(i) = ber1(i)/frame;
   ber2(i) = ber2(i)/frame;
   ber3(i) = ber3(i)/frame;
   
end % for i

% Plot the result
semilogy(EbN0, ber0, 'r');
hold;
semilogy(EbN0, ber1, 'r');
%semilogy(EbN0, ber3, 'r');
semilogy(EbN0, ber2, 'o--');
grid on;
hold off;