%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo Normal Distribution
% Written for GNU Octave
%
% INSTITUTE OF TELECOMMUNICATIONS
% University of Stuttgart
% www.inue-uni-stuttgart.de
% author: Marvin Geiselhart
% date: 18.05.2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example script just prints the values of the given parameters.
% Make sure that the function prototype matches the prototype defined
% in the XML config file. (Parameter names are arbitrary but should be the same
% for better readability)


function [] = approx_demo(paramtype, varvalue, show_pnxi, show_pyxim1, show_pyxip1, show_pyxi, show_hists, px1, num_points, filename)
  
  max_xi = 4;
  
  % 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);
  
  p = varvalue;
  switch (paramtype)
    case "std"
    sigma = abs(p);
    case "var"
    sigma = abs(sqrt(p));
    case "snr_r"
    sigma = sqrt(10^(-p/10));
    case "snr_c"
    sigma = sqrt(10^(-p/10)/2);
  endswitch
  
  if sigma < eps
    sigma = eps;
  endif
  
  if px1 < 0
    px1 = 0;
  endif
  if px1 > 1
    px1 = 1;
  endif
  
  num_points = round(num_points);
  if num_points < 0
    num_points = 1;
  endif
  
  x = (-1:0.001:1)*abs(max_xi);
  pn = pdf(x, sigma);
  pym1 = pdf(x+1, sigma);
  pyp1 = pdf(x-1, sigma);
  py = (1-px1)*pym1 + px1*pyp1;
  
  hold on;
  
  Xrv = 2*(rand(1,num_points)<px1)-1;
  Nrv = sigma *randn(1,num_points);
  Yrv = Xrv + Nrv;
  Y_xp1 = Yrv(Xrv>0);
  Y_xm1 = Yrv(Xrv<0);
  
  
  delta_x = 0.166667;
  
  x_bins = -max_xi:delta_x:max_xi;
  
  
  ber = mean(Xrv.*Yrv<0);
  ber_ex = Q(1/sigma)
  
  if show_pnxi && show_hists
    [counts,centers] = hist(Nrv, x_bins);
    freqs = counts/sum(counts)/delta_x;
    bb = bar(centers, freqs, 'facecolor',[0.3, 0.3, 0.3],'HandleVisibility','off');
  endif
  
  if show_pyxim1 && show_hists
    [counts,centers] = hist(Y_xm1,x_bins-delta_x/3);
    freqs = counts/sum(counts)/delta_x;
    bb = bar(centers, freqs, .8, 'facecolor',[0.75,0,0],'HandleVisibility','off');
  endif
  
  if show_pyxip1 && show_hists
    [counts,centers] = hist(Y_xp1,x_bins+delta_x/3);
    freqs = counts/sum(counts)/delta_x;
    bb = bar(centers, freqs, .8, 'facecolor',[0, 0.45, 0],'HandleVisibility','off');
  endif
  
  if show_pyxi && show_hists
    [counts,centers] = hist(Yrv,x_bins);
    freqs = counts/sum(counts)/delta_x;
    bb = bar(centers, freqs, .8, 'facecolor',[0,0,0.75],'HandleVisibility','off');
  endif
  
  legendCounter = 1;
  legend_entries = {};
  
  if show_pnxi
    plot(x,pn,'color',[0.5, 0.5, 0.5],'linewidth',5);
    legend_entries{legendCounter} = 'p_n(\xi)';
    legendCounter = legendCounter + 1;
  endif
  if show_pyxim1
    plot(x,pym1,'color',[1,0,0],'linewidth',5);
    legend_entries{legendCounter} = 'p_y(\xi|X=-1)';
    legendCounter = legendCounter + 1;
  endif
  
  if show_pyxip1
    plot(x,pyp1,'color',[0, 0.7, 0],'linewidth',5);
    legend_entries{legendCounter} = 'p_y(\xi|X=+1)';
    legendCounter = legendCounter + 1;
  endif
  
  if show_pyxi
    plot(x,py,'color',[0,0,1],'linewidth',5);
    legend_entries{legendCounter} = 'p_y(\xi)';
    legendCounter = legendCounter + 1;
  endif
  
  
  xlim([-max_xi, max_xi]);
  
  snr_r = 10*log10(1/(sigma^2));
  snr_c = 10*log10(1/(2*sigma^2));
  
  switch (paramtype)
    case "std"
    title1 = ['\sigma = ' num2str(sigma,'%4.2f')]
    case "var"
    title1 = ['\sigma^2 = ' num2str(sigma^2,'%4.2f')]
    case "snr_r"
    title1 = ['SNR (real-val.) = ' num2str(snr_r,'%4.2f') ' dB']
    case "snr_c"
    title1 = ['SNR (complex-val.) = ' num2str(snr_c,'%4.2f') ' dB']
  endswitch
  
  title([title1 ', BER = ' num2str(ber) ' (simulated), BER = ' num2str(ber_ex) ' (exact)']);
  xlabel('\xi');
  legend(legend_entries);
  
  
  grid on;
  
  % Create an image file. This png image is displayed on the webpage. Its size
  % must be 640x480.
  print(filename, "-dpng", output_format);
  
  
end

function y = pdf(x, sigma)
  y=  1/(sqrt(2*pi)*sigma)*exp(-1/(2*sigma^2)*x.^2)
endfunction

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


