% WEBDEMO constellation_from_snr_animation.m-file, Institute of Telecommunications, University of Stuttgart, Germany
%
% Creates a constellation diagram from a given SNR value.
%
% www.inue.uni-stuttgart.de
% author: Michael Bernhard, michibernhard@gmx.de
% date: 15.04.2014
%
% Written for GNU Octave 
%
% EXAMPLE:  
%
%
% PARAMETERS:

%  qam_input,		   % e.g. 4, 16, 64, ...
%  snrStart_input, % in dB
%  snrStop_input,  % in dB
%  plot_style,     % in dB
%  svg_en,         % if 1, then also output svg-file
%  webdemo_en,     % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
%  filename        % output file for webdemo


% Beschreibung:
% -------------------------------------------------------------------------
% Aus dem SNR wird zuerst die Standardabweichung berechnet und damit dann
% die Gaußverteilung. Der Erwartungswert "mu" wird variiert, also bei QPSK 
% m=-1 und mu=1. Alle Gaußkurven werden dann aufaddiert und das ist der
% erste Plot. Daraus wird dann dieses 2-D-Diagramm erstellt.
% -------------------------------------------------------------------------

function constellation_from_snr(
  qam_input,	   	% e.g. 4, 16, 64, ...
  snrStart_input,	% in dB
  snrStop_input,  % in dB
  plot_style,     % in dB
  svg_en,         % if 1, then also output svg-file
  webdemo_en,     % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
  filename        % output file for webdemo
)

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, set up graphical output
  

  % Create an invisible figure.
  if (webdemo_en)
    fig = figure(1); set(fig, "visible", "off");
  else
    fig = figure(1); set(fig, "visible", "on");
  endif


  landscape = "-S930,350";
  portrait  = "-S640,480";

  output_format = portrait;     % default


  set(0, "defaultlinelinewidth", 2);


  if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
  endif

  hold off;
  % END, set up graphical output
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  qam = qam_input;

  numberOfPictures = 16;
  snr_list = linspace(snrStart_input, snrStop_input, numberOfPictures);
  max_value = 0;
  for index=0:numberOfPictures-1
    snr = snr_list(index+1);
    newFilename = [filename '_' num2str(index)];
  
    range_min = -1.5;
    range_max = 1.5;
    x=[range_min:0.005:range_max]; %x=[range_min:0.04:range_max];
    if qam > 2
      pointsPerLine = sqrt(qam);
    else
      pointsPerLine = 2;
    end
    mu = [-1:2/(pointsPerLine-1):1];
    sigma = getSigma(snr);
    f = gaussPlot(x, mu, sigma);
    f2 = zeros(1,size(x,2));
    f2 = sum(f);
  
    f3 = zeros(size(x,2));
    if qam > 2
      for k=1:size(x,2)
        for n=1:size(x,2)        
          f3(k,n)=f2(1,k) * f2(1,n);
        end
      end
    else % BPSK
      mu_psk = [0 0];
      sigma_psk = getSigma(snr);
      f_psk = gaussPlot(x, mu_psk, sigma_psk);
      f2_psk = zeros(1,size(x,2));
      f2_psk = sum(f_psk);
      for k=1:size(x,2)
        for n=1:size(x,2)        
          f3(k,n)=0.5 * f2(1,n) * f2_psk(1,k);
        end
      end
    end
    %colormap('jet');
    %colormap(flipud(colormap));
    my_color_map = jet(256); 
    M = 24;
    my_color_map(1:M,1) = linspace(1,0,M); % red
    my_color_map(1:M,2) = linspace(1,0,M); % green
    my_color_map(1:M,3) = linspace(1,1,M); % blue
    colormap(my_color_map);
    if plot_style == 0
      imagesc (f3);
  
      caxis([0 1]);
      axis equal;
      colorbar
    else
      surf(f3);
      %if index==0
      %  max_value = max(max(f3)) * 0.8;
      %end
      caxis([0 1]);
      zlim([0 1.5]);
      shading flat
      
    end
    if qam > 2
      title(['probabilty density (SNR = ', num2str(snr), ' dB, ', num2str(qam), '-QAM)']);
    else
      title(['probabilty density (SNR = ', num2str(snr), ' dB, BPSK)']);
    end 
    set(gca,'XTick',[100 300 500]); %set(gca,'XTick',[12.5 37.5 62.5]); % 37.5 = 1.5 / 0.04
    set(gca,'YTick',[100 300 500]); %set(gca,'YTick',[12.5 37.5 62.5]); 
    set(gca,'XTickLabel', {'-1', '0', '1'});
    set(gca,'YTickLabel', {'-1', '0', '1'});
    xlabel('real part');
    ylabel('imaginary part');
  
    % Create an image file. This image is displayed on the webpage. Its size
    % must be 640x480.
    %print(newFilename, "-dpng",output_format);
    print('-dpng', 'TextAlphaBits=4', 'GraphicsAlphaBits=4', 'mono', output_format, newFilename);
   
    if (svg_en)
      print(["fig_",newFilename], "-dsvg",output_format);
    endif
  end   
end


function sigma = getSigma(snr)
  sigma = 10^(-snr/20);
end

function f=gaussPlot(x, mu, sigma)
for k = 1:size(mu, 2)
		%f(k,:)= 1/(sqrt(2*pi)*sigma)*exp(-0.5.*((x-mu(1,k))./sigma).^2);
	
	% normalization maximal value to 1
	f(k,:)= exp(-0.5.*((x-mu(1,k))./sigma).^2);
end
end

