%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Insert title here"
% Written for GNU Octave
%
% INSTITUTE OF TELECOMMUNICATIONS
% University of Stuttgart
% www.inue-uni-stuttgart.de
% author: Marvin Geiselhart
% date: 28.04.2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Caching documentation:
% We save all curves for a real-valued SNR_b (Eb/sigma²) from 0 dB to 15 dB in .5 dB steps, 31 values
% we do not save the snr points!
% For each of the four error rates (Hard/Soft FER/BER) we save the error rate (0=NaN) and the number of errors collected (to estimate accuracy)
% If someone simulates a code in a region with more errors than cached, we update the cache
% Directory: "error_curves_cache/"
% Filenames: %codename%_hfer.m, %codename%_sber.m, %codename%_sber_errors.m, etc


function [] = ml_decoding_demo(code1, code2, code3, show_uncoded, show_fer, show_hd, snr_def, snr_start, snr_end, y_lim, num_errs, maxtime, filename)

  snr_b_real_db_points = 0:0.5:15 % do not change this
  
  % define formats
  landscape = "-S930,350";
  portrait  = "-S640,480";
  
  % select format
  output_format = portrait;
  
  % Create an invisible figure.
  fh = figure(1);
  set(fh, "visible", "off");
  
  % define style
  set(0, "defaultlinelinewidth", 2);
  set(0, "defaultaxesfontsize", 8);
  set(0, "defaultaxesfontname", "Helvetica");
  set(0, "defaulttextfontsize", 8);
  set(0, "defaulttextfontname", "Helvetica");
  
  % This to test and demonstrate the debug engine:
  %fprintf("This is a debug message. Checking format: %s", output_format);
  
  % swap limits if out of order
  if snr_start > snr_end
    tmp = snr_start;
    snr_start = snr_end;
    snr_end = tmp;
  endif
  
  snr_start_p = snr_start;
  snr_end_p = snr_end;
  clf;
  hold on;
  xlim([snr_start,snr_end]);
  
  % setup legend
  legendCounter = 1;
  legend_entries = {};
  if show_uncoded
    % uncoded
    x_vals = linspace(snr_start, snr_end, 31);
    switch (snr_def)
    case {"eb_real", "es_real"}
      y = Perr(x_vals);
    case {"eb_complex", "es_complex"}
      y = Perr(x_vals+dB(2));
    endswitch
    semilogy(x_vals, y, "color", 'k');
    legend_entries{legendCounter} = "Uncoded BPSK";
    legendCounter = legendCounter + 1;
  endif
  
  
  for i = 1:3
    % select a code
    if (i == 1)
      code = code1
      color = 'r'
    elseif (i == 2)
      code = code2
      color = [0,0.7,0]
    elseif (i == 3)
      code = code3
      color = 'b'
    endif
    nonlinear = false;
    switch (code)
    case "rep"
      G = [1 1];
      sys_bits = [1];
      code_name = "Repetition Code"
    case "four_two"
      G = [1 0 1 1; 0 1 1 0];
      sys_bits = [1,2];
      code_name = "(4,2) Code"
    case "six_three"
      G = [1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1];
      sys_bits = [1,2,3];
      code_name = "(6,3) Code"
    case "ext_hamming"
        code_name = "(8,4) Ext. Hamming Code"
        G = [1 0 0 0 0 1 1 1; 0 1 0 0 1 0 1 1; 0 0 1 0 1 1 0 1; 0 0 0 1 1 1 1 0];
        short_name = 'eHam(8,4)';
        sys_bits = [1, 2, 3, 4];
    case "hamming_seven"
        code_name = "(7,4) Hamming Code"
        G = [1 0 0 0 0 1 1;
           0 1 0 0 1 0 1;
           0 0 1 0 1 1 0;
           0 0 0 1 1 1 1];
        sys_bits = [1, 2, 3, 4];
        short_name = 'Ham(7,4)';
    case "hamming_fifteen"
        code_name = "(15,11) Hamming Code"
        G = [1 0 0 0 0 0 0 0 0 0 0 0 0 1 1;
           0 1 0 0 0 0 0 0 0 0 0 0 1 0 1;
           0 0 1 0 0 0 0 0 0 0 0 0 1 1 0;
           0 0 0 1 0 0 0 0 0 0 1 0 0 0 1;
           0 0 0 0 1 0 0 0 0 0 1 0 0 1 0;
           0 0 0 0 0 1 0 0 0 0 1 0 1 0 0;
           0 0 0 0 0 0 1 0 0 0 1 0 1 1 1;
           0 0 0 0 0 0 0 1 0 0 1 0 1 1 0;
           0 0 0 0 0 0 0 0 1 0 1 0 1 0 1;
           0 0 0 0 0 0 0 0 0 1 1 0 0 1 1;
           0 0 0 0 0 0 0 0 0 0 0 1 1 1 1];
        sys_bits = [1,2,3,4,5,6,7,8,9,10,12];
        short_name = 'Ham(15,11)';
    case "rm_one_four"
      code_name = "(16,5) Reed-Muller Code"
      G = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;      0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1;      0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1;      0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1;      0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1];
      sys_bits = [1,2,3,5,9];
      
    case "polar"
      code_name = "(16,8) Polar Code"
      G = [1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0;
      1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0;
      1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0;
      1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0;
      1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;
      1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0;
      1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0;
      1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
      sys_bits = [16,15,14,13,12,11,10,8];
   case "ldpc"
      code_name = "(16,8) LDPC Code"
      G = [1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1;
        0 1 0 0 0 0 0 0 0 0 1 0 1 1 1 0;
        0 0 1 0 0 0 0 0 1 1 1 1 0 1 1 0;
        0 0 0 1 0 0 0 0 1 1 1 1 1 0 0 1;
        0 0 0 0 1 0 0 0 1 1 1 0 0 0 0 1;
        0 0 0 0 0 1 0 0 1 1 0 1 0 0 1 0;
        0 0 0 0 0 0 1 0 1 0 0 1 1 1 1 1;
        0 0 0 0 0 0 0 1 0 1 1 0 1 1 1 1];
      sys_bits = [1,2,3,4,5,6,7,8];   
    case "egolay"
      code_name = "(24,12) Ext. Golay Code"
      G = [1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 1 1 0 0 0 1 1;
      0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 1 0 0 1 0;
      0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 0 1 0 1 1;
      0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 1 1 1 0 1 1 0;
      0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 1 1 0 1 1 0 0 1;
      0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 1 1 0 1 1 0 1;
      0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 1 1 0 1 1 1;
      0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 1 0 1 1 1 1 0 0 0;
      0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 1 0 1 1 1 1 0 0;
      0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 1 0 1 1 1 1 0;
      0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 1 1 0 0 0 1 1 0 1;
      0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 1 1 0 0 0 1 1 1];
      sys_bits = [1,2,3,4,5,6,7,8,9,10,11,12];
    case "bch"
      code_name = "(31,16) BCH Code"
      G = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 1 0 1 0 1 1 1 0 1;
      0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1;
      0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 1 0 0 1 0 0;
      0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 1 0 0 1 0;
      0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 1 0 0 1;
      0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 1 1 0 0 1 0 0 1;
      0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 1 1 0 1 1 1 0 0 1;
      0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 0 0 0 0 0 0 1;
      0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 1 1 1 0 1 0 0 1 1 1 0 1;
      0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 1 0 0 1 1;
      0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 0 0;
      0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 0;
      0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1;
      0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 1 1 1;
      0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 1 0 1 0 1 1 1 0 1 1 0;
      0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 1 0 1 0 1 1 1 0 1 1];
      sys_bits = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16];
      
      %load("error_curves/bch_snr_b_real_cached.mat"); % snr_vals_db
      %snr_b_real_cached = snr_vals_db
      %load("error_curves/bch_sber.mat");
      %load("error_curves/bch_hber.mat");
      %load("error_curves/bch_sfer.mat");
      %load("error_curves/bch_hfer.mat");
      %use_cached = true;
    case "sixteen_eight_five"
      code_name = "(16,8) Distance-5 Code"
      G = [1 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0;
      0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 1;
      0 0 1 0 0 0 0 0 0 1 1 0 1 1 0 0;
      0 0 0 1 0 0 0 0 0 0 1 0 1 0 1 1;
      0 0 0 0 1 0 0 0 1 1 0 0 1 0 1 0;
      0 0 0 0 0 1 0 0 0 0 1 1 0 1 1 0;
      0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 1;
      0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1];
      sys_bits = [1,2,3,4,5,6,7,8];
    case "bluetooth"
      code_name = "(15,10) Bluetooth Code"
      G = [1 0 0 0 0 0 0 0 0 0 0 1 0 1 1;
      0 1 0 0 0 0 0 0 0 0 1 0 1 1 0;
      0 0 1 0 0 0 0 0 0 0 0 0 1 1 1;
      0 0 0 1 0 0 0 0 0 0 0 1 1 1 0;
      0 0 0 0 1 0 0 0 0 0 1 1 1 0 0;
      0 0 0 0 0 1 0 0 0 0 1 0 0 1 1;
      0 0 0 0 0 0 1 0 0 0 0 1 1 0 1;
      0 0 0 0 0 0 0 1 0 0 1 1 0 1 0;
      0 0 0 0 0 0 0 0 1 0 1 1 1 1 1;
      0 0 0 0 0 0 0 0 0 1 1 0 1 0 1];
      sys_bits = [1,2,3,4,5,6,7,8,9,10];
    case "NR"
      code_name = "(16,8,6) nonlinear NR-code"
      cwt = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
      0 0 0 0 0 0 0 1 1 1 1 0 0 1 0 1;
      0 0 0 0 0 0 1 1 0 0 1 1 1 1 1 1;
      0 0 0 0 0 0 1 0 1 1 0 1 1 0 1 0;
      0 0 0 0 0 1 0 0 1 0 1 0 1 0 1 1;
      0 0 0 0 0 1 0 1 0 1 1 1 0 0 1 0;
      0 0 0 0 0 1 1 1 1 0 0 1 0 1 0 0;
      0 0 0 0 0 1 1 0 0 1 0 0 1 1 0 1;
      0 0 0 0 1 1 0 0 1 1 1 1 1 1 0 0;
      0 0 0 0 1 1 0 1 0 0 0 1 1 0 0 1;
      0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1;
      0 0 0 0 1 1 1 0 0 0 1 0 0 1 1 0;
      0 0 0 0 1 0 0 0 0 1 0 1 0 1 1 1;
      0 0 0 0 1 0 0 1 1 0 0 0 1 1 1 0;
      0 0 0 0 1 0 1 1 0 1 1 0 1 0 0 0;
      0 0 0 0 1 0 1 0 1 0 1 1 0 0 0 1;
      0 0 0 1 0 0 0 0 0 1 1 1 1 0 0 1;
      0 0 0 1 0 0 0 1 1 0 0 1 0 0 1 1;
      0 0 0 1 0 0 1 1 0 1 0 0 0 1 1 0;
      0 0 0 1 0 0 1 0 1 0 1 0 1 1 0 0;
      0 0 0 1 0 1 0 0 0 0 0 1 1 1 1 0;
      0 0 0 1 0 1 0 1 1 1 0 0 1 0 0 0;
      0 0 0 1 0 1 1 1 0 0 1 0 0 0 0 1;
      0 0 0 1 0 1 1 0 1 1 1 1 0 1 1 1;
      0 0 0 1 1 1 0 0 1 0 0 0 0 1 0 1;
      0 0 0 1 1 1 0 1 0 1 1 0 1 1 1 1;
      0 0 0 1 1 1 1 1 1 0 1 1 1 0 1 0;
      0 0 0 1 1 1 1 0 0 1 0 1 0 0 0 0;
      0 0 0 1 1 0 0 0 1 1 1 0 0 0 1 0;
      0 0 0 1 1 0 0 1 0 0 1 1 0 1 0 0;
      0 0 0 1 1 0 1 1 1 1 0 1 1 1 0 1;
      0 0 0 1 1 0 1 0 0 0 0 0 1 0 1 1;
      0 0 1 1 0 0 0 0 1 1 0 0 1 1 1 1;
      0 0 1 1 0 0 0 1 0 0 1 0 1 0 1 0;
      0 0 1 1 0 0 1 1 1 1 1 1 0 0 0 0;
      0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 1;
      0 0 1 1 0 1 0 0 0 1 1 0 0 1 0 0;
      0 0 1 1 0 1 0 1 1 0 1 1 1 1 0 1;
      0 0 1 1 0 1 1 1 0 1 0 1 1 0 1 1;
      0 0 1 1 0 1 1 0 1 0 0 0 0 0 1 0;
      0 0 1 1 1 1 0 0 0 0 1 1 0 0 1 1;
      0 0 1 1 1 1 0 1 1 1 0 1 0 1 1 0;
      0 0 1 1 1 1 1 1 0 0 0 0 1 1 0 0;
      0 0 1 1 1 1 1 0 1 1 1 0 1 0 0 1;
      0 0 1 1 1 0 0 0 1 0 0 1 1 0 0 0;
      0 0 1 1 1 0 0 1 0 1 0 0 0 0 0 1;
      0 0 1 1 1 0 1 1 1 0 1 0 0 1 1 1;
      0 0 1 1 1 0 1 0 0 1 1 1 1 1 1 0;
      0 0 1 0 0 0 0 0 1 0 1 1 0 1 1 0;
      0 0 1 0 0 0 0 1 0 1 0 1 1 1 0 0;
      0 0 1 0 0 0 1 1 1 0 0 0 1 0 0 1;
      0 0 1 0 0 0 1 0 0 1 1 0 0 0 1 1;
      0 0 1 0 0 1 0 0 1 1 0 1 0 0 0 1;
      0 0 1 0 0 1 0 1 0 0 0 0 0 1 1 1;
      0 0 1 0 0 1 1 1 1 1 1 0 1 1 1 0;
      0 0 1 0 0 1 1 0 0 0 1 1 1 0 0 0;
      0 0 1 0 1 1 0 0 0 1 0 0 1 0 1 0;
      0 0 1 0 1 1 0 1 1 0 1 0 0 0 0 0;
      0 0 1 0 1 1 1 1 0 1 1 1 0 1 0 1;
      0 0 1 0 1 1 1 0 1 0 0 1 1 1 1 1;
      0 0 1 0 1 0 0 0 0 0 1 0 1 1 0 1;
      0 0 1 0 1 0 0 1 1 1 1 1 1 0 1 1;
      0 0 1 0 1 0 1 1 0 0 0 1 0 0 1 0;
      0 0 1 0 1 0 1 0 1 1 0 0 0 1 0 0;
      0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1;
      0 1 0 1 0 1 1 1 1 0 0 0 1 1 1 1;
      0 1 0 1 0 1 1 0 0 1 1 0 1 0 1 0;
      0 1 0 1 0 1 0 0 1 0 1 1 0 0 0 0;
      0 1 0 1 1 1 0 1 0 0 0 0 0 0 1 0;
      0 1 0 1 1 1 1 1 1 1 1 0 0 1 0 0;
      0 1 0 1 1 1 1 0 0 0 1 1 1 1 0 1;
      0 1 0 1 1 1 0 0 1 1 0 1 1 0 1 1;
      0 1 0 1 1 0 0 1 1 0 1 0 1 0 0 1;
      0 1 0 1 1 0 1 1 0 1 1 1 0 0 1 1;
      0 1 0 1 1 0 1 0 1 0 0 1 0 1 1 0;
      0 1 0 1 1 0 0 0 0 1 0 0 1 1 0 0;
      0 1 0 1 0 0 0 1 1 1 1 1 1 1 1 0;
      0 1 0 1 0 0 1 1 0 0 0 1 1 0 0 0;
      0 1 0 1 0 0 1 0 1 1 0 0 0 0 0 1;
      0 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1;
      0 1 1 1 0 1 0 1 1 1 1 0 0 0 1 1;
      0 1 1 1 0 1 1 1 0 0 1 1 0 1 1 0;
      0 1 1 1 0 1 1 0 1 1 0 1 1 1 0 0;
      0 1 1 1 0 1 0 0 0 0 0 0 1 0 0 1;
      0 1 1 1 1 1 0 1 0 1 1 1 1 0 0 0;
      0 1 1 1 1 1 1 1 1 0 0 1 0 0 0 1;
      0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1;
      0 1 1 1 1 1 0 0 1 0 1 0 1 1 1 0;
      0 1 1 1 1 0 0 1 0 0 0 1 1 1 1 1;
      0 1 1 1 1 0 1 1 1 1 0 0 1 0 1 0;
      0 1 1 1 1 0 1 0 0 0 1 0 0 0 0 0;
      0 1 1 1 1 0 0 0 1 1 1 1 0 1 0 1;
      0 1 1 1 0 0 0 1 1 0 0 0 0 1 0 0;
      0 1 1 1 0 0 1 1 0 1 1 0 1 1 0 1;
      0 1 1 1 0 0 1 0 1 0 1 1 1 0 1 1;
      0 1 1 1 0 0 0 0 0 1 0 1 0 0 1 0;
      0 1 1 0 0 1 0 1 1 0 0 1 1 0 1 0;
      0 1 1 0 0 1 1 1 0 1 0 0 0 0 0 0;
      0 1 1 0 0 1 1 0 1 0 1 0 0 1 0 1;
      0 1 1 0 0 1 0 0 0 1 1 1 1 1 1 1;
      0 1 1 0 1 1 0 1 1 1 0 0 1 1 0 1;
      0 1 1 0 1 1 1 1 0 0 1 0 1 0 1 1;
      0 1 1 0 1 1 1 0 1 1 1 1 0 0 1 0;
      0 1 1 0 1 1 0 0 0 0 0 1 0 1 0 0;
      0 1 1 0 1 0 0 1 0 1 1 0 0 1 1 0;
      0 1 1 0 1 0 1 1 1 0 1 1 1 1 0 0;
      0 1 1 0 1 0 1 0 0 1 0 1 1 0 0 1;
      0 1 1 0 1 0 0 0 1 0 0 0 0 0 1 1;
      0 1 1 0 0 0 0 1 0 0 1 1 0 0 0 1;
      0 1 1 0 0 0 1 1 1 1 0 1 0 1 1 1;
      0 1 1 0 0 0 1 0 0 0 0 0 1 1 1 0;
      0 1 1 0 0 0 0 0 1 1 1 0 1 0 0 0;
      0 1 0 0 0 1 0 1 0 0 1 0 1 1 0 0;
      0 1 0 0 0 1 1 1 1 1 1 1 1 0 0 1;
      0 1 0 0 0 1 1 0 0 0 0 1 0 0 1 1;
      0 1 0 0 0 1 0 0 1 1 0 0 0 1 1 0;
      0 1 0 0 1 1 0 1 1 0 1 1 0 1 1 1;
      0 1 0 0 1 1 1 1 0 1 0 1 1 1 1 0;
      0 1 0 0 1 1 1 0 1 0 0 0 1 0 0 0;
      0 1 0 0 1 1 0 0 0 1 1 0 0 0 0 1;
      0 1 0 0 1 0 0 1 1 1 0 1 0 0 0 0;
      0 1 0 0 1 0 1 1 0 0 0 0 0 1 0 1;
      0 1 0 0 1 0 1 0 1 1 1 0 1 1 1 1;
      0 1 0 0 1 0 0 0 0 0 1 1 1 0 1 0;
      0 1 0 0 0 0 0 1 0 1 0 0 1 0 1 1;
      0 1 0 0 0 0 1 1 1 0 1 0 0 0 1 0;
      0 1 0 0 0 0 1 0 0 1 1 1 0 1 0 0;
      0 1 0 0 0 0 0 0 1 0 0 1 1 1 0 1;
      1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
      1 1 1 1 1 1 1 0 0 0 0 1 1 0 1 0;
      1 1 1 1 1 1 0 0 1 1 0 0 0 0 0 0;
      1 1 1 1 1 1 0 1 0 0 1 0 0 1 0 1;
      1 1 1 1 1 0 1 1 0 1 0 1 0 1 0 0;
      1 1 1 1 1 0 1 0 1 0 0 0 1 1 0 1;
      1 1 1 1 1 0 0 0 0 1 1 0 1 0 1 1;
      1 1 1 1 1 0 0 1 1 0 1 1 0 0 1 0;
      1 1 1 1 0 0 1 1 0 0 0 0 0 0 1 1;
      1 1 1 1 0 0 1 0 1 1 1 0 0 1 1 0;
      1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0;
      1 1 1 1 0 0 0 1 1 1 0 1 1 0 0 1;
      1 1 1 1 0 1 1 1 1 0 1 0 1 0 0 0;
      1 1 1 1 0 1 1 0 0 1 1 1 0 0 0 1;
      1 1 1 1 0 1 0 0 1 0 0 1 0 1 1 1;
      1 1 1 1 0 1 0 1 0 1 0 0 1 1 1 0;
      1 1 1 0 1 1 1 1 1 0 0 0 0 1 1 0;
      1 1 1 0 1 1 1 0 0 1 1 0 1 1 0 0;
      1 1 1 0 1 1 0 0 1 0 1 1 1 0 0 1;
      1 1 1 0 1 1 0 1 0 1 0 1 0 0 1 1;
      1 1 1 0 1 0 1 1 1 1 1 0 0 0 0 1;
      1 1 1 0 1 0 1 0 0 0 1 1 0 1 1 1;
      1 1 1 0 1 0 0 0 1 1 0 1 1 1 1 0;
      1 1 1 0 1 0 0 1 0 0 0 0 1 0 0 0;
      1 1 1 0 0 0 1 1 0 1 1 1 1 0 1 0;
      1 1 1 0 0 0 1 0 1 0 0 1 0 0 0 0;
      1 1 1 0 0 0 0 0 0 1 0 0 0 1 0 1;
      1 1 1 0 0 0 0 1 1 0 1 0 1 1 1 1;
      1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 1;
      1 1 1 0 0 1 1 0 1 1 0 0 1 0 1 1;
      1 1 1 0 0 1 0 0 0 0 1 0 0 0 1 0;
      1 1 1 0 0 1 0 1 1 1 1 1 0 1 0 0;
      1 1 0 0 1 1 1 1 0 0 1 1 0 0 0 0;
      1 1 0 0 1 1 1 0 1 1 0 1 0 1 0 1;
      1 1 0 0 1 1 0 0 0 0 0 0 1 1 1 1;
      1 1 0 0 1 1 0 1 1 1 1 0 1 0 1 0;
      1 1 0 0 1 0 1 1 1 0 0 1 1 0 1 1;
      1 1 0 0 1 0 1 0 0 1 0 0 0 0 1 0;
      1 1 0 0 1 0 0 0 1 0 1 0 0 1 0 0;
      1 1 0 0 1 0 0 1 0 1 1 1 1 1 0 1;
      1 1 0 0 0 0 1 1 1 1 0 0 1 1 0 0;
      1 1 0 0 0 0 1 0 0 0 1 0 1 0 0 1;
      1 1 0 0 0 0 0 0 1 1 1 1 0 0 1 1;
      1 1 0 0 0 0 0 1 0 0 0 1 0 1 1 0;
      1 1 0 0 0 1 1 1 0 1 1 0 0 1 1 1;
      1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 0;
      1 1 0 0 0 1 0 0 0 1 0 1 1 0 0 0;
      1 1 0 0 0 1 0 1 1 0 0 0 0 0 0 1;
      1 1 0 1 1 1 1 1 0 1 0 0 1 0 0 1;
      1 1 0 1 1 1 1 0 1 0 1 0 0 0 1 1;
      1 1 0 1 1 1 0 0 0 1 1 1 0 1 1 0;
      1 1 0 1 1 1 0 1 1 0 0 1 1 1 0 0;
      1 1 0 1 1 0 1 1 0 0 1 0 1 1 1 0;
      1 1 0 1 1 0 1 0 1 1 1 1 1 0 0 0;
      1 1 0 1 1 0 0 0 0 0 0 1 0 0 0 1;
      1 1 0 1 1 0 0 1 1 1 0 0 0 1 1 1;
      1 1 0 1 0 0 1 1 1 0 1 1 0 1 0 1;
      1 1 0 1 0 0 1 0 0 1 0 1 1 1 1 1;
      1 1 0 1 0 0 0 0 1 0 0 0 1 0 1 0;
      1 1 0 1 0 0 0 1 0 1 1 0 0 0 0 0;
      1 1 0 1 0 1 1 1 1 1 0 1 0 0 1 0;
      1 1 0 1 0 1 1 0 0 0 0 0 0 1 0 0;
      1 1 0 1 0 1 0 0 1 1 1 0 1 1 0 1;
      1 1 0 1 0 1 0 1 0 0 1 1 1 0 1 1;
      1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0;
      1 0 1 0 1 0 0 0 0 1 1 1 0 0 0 0;
      1 0 1 0 1 0 0 1 1 0 0 1 0 1 0 1;
      1 0 1 0 1 0 1 1 0 1 0 0 1 1 1 1;
      1 0 1 0 0 0 1 0 1 1 1 1 1 1 0 1;
      1 0 1 0 0 0 0 0 0 0 0 1 1 0 1 1;
      1 0 1 0 0 0 0 1 1 1 0 0 0 0 1 0;
      1 0 1 0 0 0 1 1 0 0 1 0 0 1 0 0;
      1 0 1 0 0 1 1 0 0 1 0 1 0 1 1 0;
      1 0 1 0 0 1 0 0 1 0 0 0 1 1 0 0;
      1 0 1 0 0 1 0 1 0 1 1 0 1 0 0 1;
      1 0 1 0 0 1 1 1 1 0 1 1 0 0 1 1;
      1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 1;
      1 0 1 0 1 1 0 0 1 1 1 0 0 1 1 1;
      1 0 1 0 1 1 0 1 0 0 1 1 1 1 1 0;
      1 0 1 0 1 1 1 1 1 1 0 1 1 0 0 0;
      1 0 0 0 1 0 1 0 0 0 0 1 1 1 0 0;
      1 0 0 0 1 0 0 0 1 1 0 0 1 0 0 1;
      1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 1;
      1 0 0 0 1 0 1 1 1 1 1 1 0 1 1 0;
      1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 1;
      1 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0;
      1 0 0 0 0 0 0 1 1 0 1 1 1 0 0 0;
      1 0 0 0 0 0 1 1 0 1 0 1 0 0 0 1;
      1 0 0 0 0 1 1 0 1 1 1 0 0 0 0 0;
      1 0 0 0 0 1 0 0 0 0 1 1 0 1 0 1;
      1 0 0 0 0 1 0 1 1 1 0 1 1 1 1 1;
      1 0 0 0 0 1 1 1 0 0 0 0 1 0 1 0;
      1 0 0 0 1 1 1 0 0 1 1 1 1 0 1 1;
      1 0 0 0 1 1 0 0 1 0 0 1 0 0 1 0;
      1 0 0 0 1 1 0 1 0 1 0 0 0 1 0 0;
      1 0 0 0 1 1 1 1 1 0 1 0 1 1 0 1;
      1 0 0 1 1 0 1 0 0 1 1 0 0 1 0 1;
      1 0 0 1 1 0 0 0 1 0 1 1 1 1 1 1;
      1 0 0 1 1 0 0 1 0 1 0 1 1 0 1 0;
      1 0 0 1 1 0 1 1 1 0 0 0 0 0 0 0;
      1 0 0 1 0 0 1 0 0 0 1 1 0 0 1 0;
      1 0 0 1 0 0 0 0 1 1 0 1 0 1 0 0;
      1 0 0 1 0 0 0 1 0 0 0 0 1 1 0 1;
      1 0 0 1 0 0 1 1 1 1 1 0 1 0 1 1;
      1 0 0 1 0 1 1 0 1 0 0 1 1 0 0 1;
      1 0 0 1 0 1 0 0 0 1 0 0 0 0 1 1;
      1 0 0 1 0 1 0 1 1 0 1 0 0 1 1 0;
      1 0 0 1 0 1 1 1 0 1 1 1 1 1 0 0;
      1 0 0 1 1 1 1 0 1 1 0 0 1 1 1 0;
      1 0 0 1 1 1 0 0 0 0 1 0 1 0 0 0;
      1 0 0 1 1 1 0 1 1 1 1 1 0 0 0 1;
      1 0 0 1 1 1 1 1 0 0 0 1 0 1 1 1;
      1 0 1 1 1 0 1 0 1 1 0 1 0 0 1 1;
      1 0 1 1 1 0 0 0 0 0 0 0 0 1 1 0;
      1 0 1 1 1 0 0 1 1 1 1 0 1 1 0 0;
      1 0 1 1 1 0 1 1 0 0 1 1 1 0 0 1;
      1 0 1 1 0 0 1 0 0 1 0 0 1 0 0 0;
      1 0 1 1 0 0 0 0 1 0 1 0 0 0 0 1;
      1 0 1 1 0 0 0 1 0 1 1 1 0 1 1 1;
      1 0 1 1 0 0 1 1 1 0 0 1 1 1 1 0;
      1 0 1 1 0 1 1 0 0 0 1 0 1 1 1 1;
      1 0 1 1 0 1 0 0 1 1 1 1 1 0 1 0;
      1 0 1 1 0 1 0 1 0 0 0 1 0 0 0 0;
      1 0 1 1 0 1 1 1 1 1 0 0 0 1 0 1;
      1 0 1 1 1 1 1 0 1 0 1 1 0 1 0 0;
      1 0 1 1 1 1 0 0 0 1 0 1 1 1 0 1;
      1 0 1 1 1 1 0 1 1 0 0 0 1 0 1 1;
      1 0 1 1 1 1 1 1 0 1 1 0 0 0 1 0];
      sys_bits = [1,2,3,4,5,6,7,8]; % assume we can use these as linear independent columns
      systematic = false;
      nonlinear = true;
      K = 8;
      N = 16;
      short_name = 'NR';
    case "CCSDS"
    code_name = "(32,16) LDPC Code"
    G = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 1 0 0 1 1 1 1 1 0 0;
        0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1 0 1 1 0 1 1 0;
        0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 1 1 0 1 0 0 1 1;
        0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 1 1 1 0 1 0 0 1;
        0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 0 1 1 0 0 1 0 1 1;
        0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 1 0 1 0 1 1 0 1 1 0 1;
        0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 0 0 0 1 1 1 1 1 0;
        0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 1 1 1;
        0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 1 1 0 1 0 0 1 0 1 1 0 1 0;
        0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 1 0 1 0 1 1 0 1 0 0 1 0 1;
        0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 1 0 1 0 0 1 0 1 1 0 1 0;
        0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 1 0 1;
        0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 1 1 0 0 0 1 0 1 1 0 0 1 0 1;
        0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 0 1 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 0 0 1 0 0 1 0 1 0 1;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 1 0 1 1 0 0 1 0 1 0];
    systematic = true;
    sys_bits = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16];
    short_name = 'LDPC';
    otherwise
      continue
    endswitch
    
    % generate codeword table
    if !nonlinear
      [K,N] = size(G);
      A = dec2bin(0:2^length(sys_bits) - 1) - '0';
      cwt_bpsk = 1-2*mod(A*G, 2);
    else
      cwt_bpsk = 1-2*cwt;
    endif
    Rc = K/N;
    
    % Map axis to correct region
    switch (snr_def)
    case "eb_real"
      offset = 0;
    case "es_real"
      offset = dB(Rc);
    case "eb_complex"
      offset = -dB(2);
    case "es_complex"
      offset = dB(Rc)-dB(2);
    endswitch
    
    
    x_vals = snr_b_real_db_points+offset;
    num_points = length(snr_b_real_db_points);
    start_index = max(floor((snr_start_p-offset)*2), 0)+1
    end_index = min(ceil((snr_end_p-offset)*2), 2*15)+1

    %x_vals = snr_vals_db; % always use the defined range

    
    % snr loop
    
    bit_error_rate = zeros([1,num_points]);
    frame_error_rate = zeros([1,num_points]);
    bit_errs = zeros([1,num_points]);
    frame_errs = zeros([1,num_points]);
    
    
    
    if K <= 8
      batch_size = 512
    else
      batch_size = 64
    endif
    tic;
    t = 0;
    
    for i = start_index:end_index
      rand ("state", 0)
      snr_db = snr_b_real_db_points(i)
      snr_lin_s = Rc*10^(snr_db/10.0);
      sn=sqrt(1/snr_lin_s);   
      num_bits = 0;
      num_frames = 0;
      while bit_errs(i) < num_errs || num_bits < 10000
        %if nonlinear  % select batch_size many random codewords
        rndIDX = randi(2^K, batch_size, 1); 
        x = cwt_bpsk(rndIDX,:);
        %else
        %  x = ones(batch_size,N);    % all zero codeword suffices (linear code)
        %endif  

        n = sn*randn(batch_size,N); % generate noise
        y = x + n;
        
        corr=cwt_bpsk*y';
        [corr_val, cw_index] = max(corr);
        xhat=cwt_bpsk(cw_index,:);

        frame_errs(i) = frame_errs(i) + sum(any(x!=xhat,2));
        num_frames = num_frames + batch_size;

        bit_errs(i) = bit_errs(i) + sum(sum(x(:,sys_bits)!=xhat(:,sys_bits)));
        num_bits = num_bits + K*batch_size;
        t = toc;
        if (t > maxtime)
          break
        endif

      endwhile
      if bit_errs(i) > num_errs/8 % ignore results of much lower than desired accuracy
        bit_error_rate(i) = bit_errs(i) /num_bits;
        frame_error_rate(i) = frame_errs(i)/num_frames;
      endif
      
      if (t > maxtime)
        break
      endif
    endfor
    update_cache(code, 'sber', bit_error_rate, bit_errs);
    update_cache(code, 'sfer', frame_error_rate, frame_errs);
    
    
    if show_hd
      tic;
      t = 0;
      hard_bit_error_rate = zeros([1,num_points]);
      hard_frame_error_rate = zeros([1,num_points]);
      hard_bit_errs = zeros([1,num_points]);
      hard_frame_errs = zeros([1,num_points]);
      for i = start_index:end_index
        rand ("state", 0)
        snr_db = snr_b_real_db_points(i)
        snr_lin_s = Rc*10^(snr_db/10.0);
        sn=sqrt(1/snr_lin_s);
        
        
        num_bits = 0;
        num_frames = 0;
        while hard_bit_errs(i) < num_errs || num_bits < 10000
          %if nonlinear  % select batch_size many random codewords
          rndIDX = randi(2^K, batch_size, 1); 
          x = cwt_bpsk(rndIDX,:);
          %else
          %  x = ones(batch_size,N);    % all zero codeword suffices (linear code)
          %endif  


          n = sn*randn(batch_size,N);
          y = x + n;
          
          hard_y = sign(y);
          hard_corr=cwt_bpsk*hard_y';
          [corr_val, hard_cw_index] = max(hard_corr);
          xhat_hard=cwt_bpsk(hard_cw_index,:);
          
          hard_frame_errs(i) = hard_frame_errs(i) + sum(any(x!=xhat_hard,2));
          hard_bit_errs(i) = hard_bit_errs(i) + sum(sum(x(:,sys_bits)!=xhat_hard(:,sys_bits)));
          num_frames = num_frames + batch_size;
          num_bits = num_bits + K*batch_size;
          t = toc;
          if (t > maxtime)
            break
          endif
        endwhile
        if hard_bit_errs(i) > num_errs/8
          hard_bit_error_rate(i) = hard_bit_errs(i)/num_bits;
          hard_frame_error_rate(i) = hard_frame_errs(i)/num_frames;
        endif
        if (t > maxtime)
          break
        endif
      endfor
      update_cache(code, 'hber', hard_bit_error_rate, hard_bit_errs);
      update_cache(code, 'hfer', hard_frame_error_rate, hard_frame_errs);
    endif
    
    % PLOT
    if sum(snr_start_p>=x_vals)
      snr_start = max(x_vals(snr_start_p>=x_vals));
    else
      snr_start = snr_start_p;
    endif
    
    if sum(snr_end_p<=x_vals)
      snr_end = min(x_vals(snr_end_p<=x_vals));
    else
      snr_end = snr_end_p;
    endif
    
    if show_hd && show_fer
      semilogy(x_vals(snr_start<=x_vals & x_vals <= snr_end), hard_frame_error_rate(snr_start<=x_vals & x_vals <= snr_end) , "color", color, "linestyle", '-.',"marker", 'x', "markersize",3);
      legend_entries{legendCounter} = [code_name ', Hard Dec. FER'];
      legendCounter = legendCounter + 1;
    endif
    
    if show_hd
      semilogy(x_vals(snr_start<=x_vals & x_vals <= snr_end), hard_bit_error_rate(snr_start<=x_vals & x_vals <= snr_end) , "color", color, "linestyle", ':');
      legend_entries{legendCounter} = [code_name ', Hard Dec. BER'];
      legendCounter = legendCounter + 1;
    endif
    
    if show_fer
      semilogy(x_vals(snr_start<=x_vals & x_vals <= snr_end), frame_error_rate(snr_start<=x_vals & x_vals <= snr_end) , "color", color, "linestyle", '--',"marker", 'o', "markersize",3);
      legend_entries{legendCounter} = [code_name ', Soft Dec. FER'];
      legendCounter = legendCounter + 1;
    endif
    
    semilogy(x_vals(snr_start<=x_vals & x_vals <= snr_end), bit_error_rate(snr_start<=x_vals & x_vals <= snr_end), "color", color);
    
    if show_hd || show_fer
      legend_entries{legendCounter} = [code_name ', Soft Dec. BER'];
    else
      legend_entries{legendCounter} = code_name;
    endif
    legendCounter = legendCounter + 1;
    
  endfor
  
  switch (snr_def)
  case "eb_real"
    xlabel('E_b/N_0 in dB (real-valued channel)');
  case "es_real"
    xlabel('E_s/N_0 in dB (real-valued channel)');
  case "eb_complex"
    xlabel('E_b/N_0 in dB (complex-valued channel)');
  case "es_complex"
    xlabel('E_s/N_0 in dB (complex-valued channel)');
  endswitch
  
  if y_lim < 0
    yyy = ylim();
    ylim([10^y_lim, yyy(2)]);
  endif
  
  if show_fer
    ylabel('BER/FER');
  else
    ylabel('BER');
  endif
  grid on;
  hLg = legend(legend_entries);
  legend("boxoff");
  %set(hLg, "Position", [0,0,10,10]);
  print(filename, "-dpng", output_format);
end

function [] = update_cache(code, simulation, error_rate_param, num_errors_param)
  nfname  = strcat('error_curves_cache/', code, '_', simulation, '_errors.mat')
  erfname = strcat('error_curves_cache/', code, '_', simulation,'.mat')
  % load previous results
  if exist(nfname)
    load( nfname); % num_errors
  else
    num_errors = zeros([1,length(num_errors_param)]);
  endif
  if exist(erfname)
    load( erfname); % error_rate
  else
    error_rate = zeros([1,length(error_rate_param)]);
  endif
  % check where it is better than before:
  update_indices = num_errors_param > num_errors;
  num_errors(update_indices) = num_errors_param(update_indices);
  error_rate(update_indices) = error_rate_param(update_indices);
  % write cache
  try
    save(nfname, "num_errors")
  catch
    delete(nfname)
    save(nfname, "num_errors")
  end_try_catch
  
  try
    save(erfname, "error_rate")
  catch
    delete(erfname)
    save(erfname, "error_rate")
  end_try_catch
  
endfunction



function y=round_down_half(x)
  y=floor(x*2)/2;
endfunction

function y=round_up_half(x)
  y=ceil(x*2)/2;
endfunction

function y=Perr(EsN0dB)
  EsN0=10.^(EsN0dB/10.0);
  y=Q(sqrt(EsN0));       % Perr
endfunction

function y=dB(x)
  y=10*log10(x);
endfunction

function y=nodB(x)
  y=10.^(x./10);
endfunction

function y = Q(x)
  y = 0.5*(1-erf(x/sqrt(2)));
endfunction
