%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .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: 25.05.2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [] = union_bound(code1, code2, code3, show_uncoded, show_ube, show_fer, snr_def, snr_start, snr_end, y_lim, filename)
  snr_b_real_db_points = 0:0.5:15 % do not change this
  
  % define formats
  landscape = "-S930,350";
  portrait  = "-S640,480";
  
  linewidth = 3;
  
  % 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);
  
  if snr_start > snr_end
    tmp = snr_start;
    snr_start = snr_end;
    snr_end = tmp;
  endif
  
  
  num_points = 25;
  
  snr_vals_db = linspace(snr_start, snr_end, num_points);
  legendCounter = 1;
  legend_entries = {};
  
  clf;
  hold on;
  
  
  xlim([snr_start,snr_end]);
  
  
  if show_uncoded
    % uncoded
    switch (snr_def)
      case {"eb_real", "es_real"}
      y = Perr(snr_vals_db);
      case {"eb_complex", "es_complex"}
      y = Perr(snr_vals_db+dB(2));
    endswitch
    semilogy(snr_vals_db, y, "color", 'k','linewidth',linewidth);
    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
    
    switch (code)
      case "rep"
      G = [1 1];
      code_name = "Repetition Code"
      case "four_two"
      G = [1 0 1 1; 0 1 1 0];
      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];
      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];
      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];
      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];
      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];
      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];
      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];
      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];
      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];
      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];
      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];
         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];
      otherwise
      continue
    endswitch
    
    % generate codeword table
    [K,N] = size(G);
    Rc = K/N;
    
    x_vals = snr_vals_db; % always use the defined range
    switch (snr_def)
      case "eb_real"
      snr_lin_s = Rc*10.^(snr_vals_db./10.0);
      case "es_real"
      snr_lin_s = 10.^(snr_vals_db./10.0);
      case "eb_complex"
      snr_lin_s = 2*Rc*10.^(snr_vals_db./10.0);
      case "es_complex"
      snr_lin_s = 2*10.^(snr_vals_db./10.0);
    endswitch
    
    sn=sqrt(1./snr_lin_s);
    
    sfname = strcat('spectra/', code, '_spectrum.mat')
    
    if (exist(sfname))
      load(sfname); % A_j
    else
      continue;
    endif
    
    
    if is_cached(code,'sfer')
      frame_error_rate = load_cached(code,'sfer')
    else
      disp("Didn't find curves");
    endif
    
    ub = zeros([1,num_points]);
    ube = zeros([1,num_points]);
    for j = N:-1:0
      A = A_j(j+1);
      if (j != 0 && A > 0)
        term = A*Q(sqrt(j)./sn);
        ub = ub + term;
        ube = term;
      endif
    endfor
    
    
    semilogy(x_vals, ub, "color", color,'linewidth',linewidth);
    legend_entries{legendCounter} = [code_name ', UB'];
    legendCounter = legendCounter + 1;
    
    if show_ube
      semilogy(x_vals, ube, "color", color,'linewidth',linewidth, "linestyle", ':');
      legend_entries{legendCounter} = [code_name ', UBE'];
      legendCounter = legendCounter + 1;
    endif
    
    if (exist("frame_error_rate", "var") == 1) && show_fer
      % compute offset
      switch (snr_def)
        case "eb_real"
        x_vals2 = snr_b_real_db_points;
        case "es_real"
        x_vals2 = snr_b_real_db_points+dB(Rc);
        case "eb_complex"
        x_vals2 = snr_b_real_db_points-dB(2);
        case "es_complex"
        x_vals2 = snr_b_real_db_points+dB(Rc)-dB(2);
      endswitch
      
      semilogy(x_vals2(snr_start<=x_vals2 & x_vals2 <= snr_end), frame_error_rate(snr_start<=x_vals2 & x_vals2 <= snr_end) , "color", color, "linestyle", '--',"marker", 'o', 'linewidth',linewidth, "markersize",3);
      legend_entries{legendCounter} = [code_name ', ML FER'];
      legendCounter = legendCounter + 1;
    endif
    
  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
  
  
  ylabel('FER');
  
  if y_lim < 0
    yyy = ylim();
    y_max = min(yyy(2),1);
    ylim([10^y_lim, y_max]);
  endif
  
  grid on;
  legend(legend_entries);
  
  
  % Create an image file. This png image is displayed on the webpage. Its size
  % must be 640x480.
  print(filename, "-dpng", output_format);
  
  
end

function error_rate = load_cached(code, simulation)
  erfname = strcat('error_curves_cache/', code, '_', simulation,'.mat')
  if exist(erfname)
    load( erfname); % error_rate
  else
    error_rate = zeros([1,length(error_rate_param)]);
  endif
endfunction

function ic = is_cached(code, simulation)
  erfname = strcat('error_curves_cache/', code, '_', simulation,'.mat')
  ic = exist(erfname)
endfunction


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


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

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

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