% ldpcut.m
%
% David Rowe 18 Dec 2013
%
% Octave LDPC unit test script using CML library, based on simulation
% by Bill Cowley VK5DSP

% Libraries we need

ldpc;
qpsk;

function sim_out = run_simulation(sim_in)

  % Note this is effective Eb/No of payload data bits, sorta thing we
  % plot on BER versus Eb/No graphs of decoded data.  So if we have a
  % rate 1/2 code, each codeword bit will have Eb/No - 3dB.

  EbNodBvec = sim_in.EbNodBvec;

  Ntrials = sim_in.Ntrials;
  verbose = sim_in.verbose;

  % Init LDPC code ------------------------------------

  mod_order = 4; bps = 2;
  modulation = 'QPSK';
  mapping = 'gray';

  demod_type = 0;
  decoder_type = 0;
  max_iterations = 100;

  if strcmp(sim_in.code,'wimax')
    rate = 0.5; framesize = 576*4;
    code_param = ldpc_init_builtin(sim_in.code, rate, framesize, modulation, mod_order, mapping);
  elseif strcmp(sim_in.code,'dvbs2')
    framesize = 16200; rate = 0.8;
    code_param = ldpc_init_builtin(sim_in.code, rate, framesize, modulation, mod_order, mapping);
    rate = code_param.ldpc_data_bits_per_frame/code_param.ldpc_coded_bits_per_frame;
  else
    % deal with H stored in different file formats
    tempStruct = load(sim_in.code);
    b = fieldnames(tempStruct);
    ldpcArrayName = b{1,1};
    % extract the array from the struct
    HRA = tempStruct.(ldpcArrayName);
    [code_param framesize rate] = ldpc_init_user(HRA, modulation, mod_order, mapping);
  end

  % optional 1's stuffing 
  if isfield(sim_in, "data_bits_per_frame")
    code_param.data_bits_per_frame = sim_in.data_bits_per_frame;
    code_param.coded_bits_per_frame = code_param.data_bits_per_frame + code_param.ldpc_parity_bits_per_frame;
    code_param.coded_syms_per_frame = code_param.coded_bits_per_frame/code_param.bits_per_symbol;
    rate = code_param.data_bits_per_frame/code_param.coded_bits_per_frame;
    printf("data_bits_per_frame = %d\n", code_param.data_bits_per_frame);
    printf("coded_bits_per_frame = %d\n", code_param.coded_bits_per_frame);
    printf("coded_syms_per_frame = %d\n", code_param.coded_syms_per_frame);
    printf("rate: %f\n",rate);
  end

  % ----------------------------------
  % run simulation at each Eb/No point
  % ----------------------------------

  for ne = 1:length(EbNodBvec)
    randn('seed',1);
    rand('seed',1);

    % Given Eb/No of payload data bits, work out Es/No we need to
    % apply to each channel symbol:
    %
    % i) Each codeword bit gets noise: Eb/No - 3 (for a rate 1/2 code)
    % ii) QPSK means two bits/symbol.: Es/No = Eb/No + 3
    %
    % -> which neatly cancel out ...... (at least for rate 1/2)

    EsNodB = EbNodBvec(ne) + 10*log10(rate) + 10*log10(bps);
    EsNo = 10^(EsNodB/10);
    variance = 1/EsNo;

    Tbits = Terrs = Ferrs = Terrs_raw = Tbits_raw = 0;

    tx_bits = [];
    tx_symbols = [];

    % Encode a bunch of frames

    for nn=1:Ntrials
      atx_bits = round(rand( 1, code_param.data_bits_per_frame));
      tx_bits = [tx_bits atx_bits];
      [tx_codeword atx_symbols] = ldpc_enc(atx_bits, code_param);
      tx_symbols = [tx_symbols atx_symbols];
    end

    rx_symbols = tx_symbols;

    % Add AWGN noise, 0.5 factor splits power evenly between Re & Im

    noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols)));
    rx_symbols += noise;

    % Decode a bunch of frames

    rx_bits_log = [];

    for nn = 1: Ntrials
      st = (nn-1)*code_param.coded_syms_per_frame + 1;
      en = (nn)*code_param.coded_syms_per_frame;

      % coded

      arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, ones(1,code_param.coded_syms_per_frame));
      st = (nn-1)*code_param.data_bits_per_frame + 1;
      en = (nn)*code_param.data_bits_per_frame;
      error_positions = xor(arx_codeword(1:code_param.data_bits_per_frame), tx_bits(st:en));
      Nerrs = sum(error_positions);
      rx_bits_log = [rx_bits_log arx_codeword(1:code_param.data_bits_per_frame)];

      % uncoded - to est raw BER compare data symbols as code is systematic

      raw_rx_bits = [];
      for s=1:code_param.coded_syms_per_frame*rate
        sym_st = (nn-1)*code_param.coded_syms_per_frame + 1;
        raw_rx_bits = [raw_rx_bits qpsk_demod(rx_symbols(sym_st+s-1))];
      end
      Nerrs_raw = sum(xor(raw_rx_bits, tx_bits(st:en)));
      Nbits_raw = code_param.data_bits_per_frame;

      if verbose == 2
        % print "." if frame decoded without errors, 'x' if we can't decode

        if Nerrs > 0, printf('x'),  else printf('.'),  end
      end

      if Nerrs > 0,  Ferrs = Ferrs + 1;  end
      Terrs     += Nerrs;
      Tbits     += code_param.ldpc_data_bits_per_frame;
      Terrs_raw += Nerrs_raw;
      Tbits_raw += Nbits_raw;
    end

    if verbose
      printf("\nCoded EbNodB: % 5.2f BER: %4.3f Tbits: %6d Terrs: %6d FER: %4.3f Tframes: %d Ferrs: %d\n",
             EbNodBvec(ne), Terrs/Tbits, Tbits, Terrs,  Ferrs/Ntrials, Ntrials, Ferrs);
      EbNodB_raw = EbNodBvec(ne) + 10*log10(rate);
      printf("Raw EbNodB..: % 5.2f BER: %4.3f Tbits: %6d Terrs: %6d\n",
             EbNodB_raw, Terrs_raw/Tbits_raw, Tbits_raw, Terrs_raw);
    end

    sim_out.rate = rate;
    sim_out.BER(ne)   = Terrs/Tbits;
    sim_out.PER(ne)   = Ferrs/Ntrials;
  end

endfunction


% ---------------------------------------------------------------------------------
% 1/ Simplest possible one frame simulation
% ---------------------------------------------------------------------------------

function test1_single(code="wimax", data_bits_per_frame)
  printf("\nTest 1:Single -----------------------------------\n");

  mod_order = 4;
  modulation = 'QPSK';
  mapping = 'gray';
  demod_type = 0;
  decoder_type = 0;
  max_iterations = 100;

  % CML library has a bunch of different framesizes available
  if strcmp(code,'wimax') framesize = 576*2; rate = 0.5; end
  if strcmp(code,'dvbs2') framesize = 16200; rate = 0.6; end
  code_param = ldpc_init_builtin(code, rate, framesize, modulation, mod_order, mapping);
    
  % optional 1's stuffing 
  if nargin == 2
    code_param.data_bits_per_frame = data_bits_per_frame;
    code_param.coded_bits_per_frame = code_param.data_bits_per_frame + code_param.ldpc_parity_bits_per_frame;
    code_param.coded_syms_per_frame = code_param.coded_bits_per_frame/code_param.bits_per_symbol;
    framesize = code_param.coded_bits_per_frame;
  end
  
  % find out what rate we actually obtained ...
  rate = code_param.data_bits_per_frame/code_param.coded_bits_per_frame;
  printf("Ndata_bits: %d Nparity_bits: %d Ncodeword_bits: %d rate: %3.2f\n",
         code_param.data_bits_per_frame, code_param.ldpc_parity_bits_per_frame, 
         code_param.coded_bits_per_frame, rate);

  % decoder needs an estimated channel EsNo (linear ratio, not dB)
  EsNo = 10;

  tx_bits = round(rand(1, code_param.data_bits_per_frame));
  [tx_codeword, qpsk_symbols] = ldpc_enc(tx_bits, code_param);
  rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, qpsk_symbols, EsNo, ones(1,length(qpsk_symbols)));

  errors_positions = xor(tx_bits, rx_codeword(1:framesize*rate));
  Nerr = sum(errors_positions);
  printf("Nerr: %d\n", Nerr);
endfunction


% ---------------------------------------------------------------------------------
% 2/ Run a bunch of trials at just one EsNo point
% ---------------------------------------------------------------------------------

function test2_multiple(code, Ntrials=100, data_bits_per_frame)
  printf("\nTest 2: Multiple: %s ----------------------------\n", code);

  % these are inputs for Wimax mode, e.g. framesize defines code used

  sim_in.code = code;
  sim_in.verbose = 2;
  sim_in.Ntrials = Ntrials;
  sim_in.EbNodBvec = 3;
  if nargin == 3
    sim_in.data_bits_per_frame = data_bits_per_frame;
  end
  run_simulation(sim_in);
end


% ---------------------------------------------------------------------------------
% 3/ Lets draw some Eb/No versus BER curves
% ---------------------------------------------------------------------------------

function test3_curves(code,fg=1,Ntrials=100)
  printf("\nTest 3: Curves: %s -------------------------------------\n", code);

  sim_in.code = code;
  sim_in.verbose = 2;
  sim_in.Ntrials = Ntrials;
  sim_in.EbNodBvec = -2:10;
  sim_out = run_simulation(sim_in);

  EbNodB = sim_in.EbNodBvec;
  uncoded_awgn_ber_theory = 0.5*erfc(sqrt(10.^(EbNodB/10)));

  figure(fg); clf; title(code);
  semilogy(EbNodB, uncoded_awgn_ber_theory,'r-+;AWGN;')
  hold on;
  semilogy(EbNodB, sim_out.BER+1E-10,'g-+;AWGN LDPC;');
  hold off;
  grid('minor')
  xlabel('Eb/No (dB)')
  ylabel('BER')
  axis([min(EbNodB) max(EbNodB) 1E-3 1])
  legend('boxoff');
end


% --------------------------------------------------------------------------------
% START SIMULATIONS
% --------------------------------------------------------------------------------

more off;
format;

% ---------------------------------------------------------------------------------
% Start CML library (see CML set up instructions in ldpc.m)
% ---------------------------------------------------------------------------------

init_cml();

% Ctest kicks off these tests using env variables
if getenv("CTEST_SINGLE")
  test1_single
  return;
end
if getenv("CTEST_ONE_STUFFING")
  test2_multiple("wimax",10,576);
  return;
end

% Uncomment and try some of these tests if you like ....
%test3_curves("H_1024_2048_4f.mat",1)
%test1_single("dvbs2")
%test3_curves("dvbs2",1,10)
%test2_multiple("wimax")
%test2_multiple("H2064_516_sparse.mat")
%test3_curves("wimax",1)
%test3_curves("H2064_516_sparse.mat",2)
%test3_curves("H_256_768_22.txt",2)
%test3_curves("H_4096_8192_3d.mat")
%test3_curves("H_212_158.mat")