% WEBDEMO demo_LLR_4PAM.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: Klaus Oestreich, oestreich@inue.uni-stuttgart.de
% date: 30.10.2013
%
% Written for GNU Octave (works for MATLAB as well)
%
% EXAMPLE:  demo_LLR_4PAM(10, 0.5, 0.5, 1, 15, 0, 1, 1,'demo_LLR_4PAM')
% PARAMETERS:
%  EsN0,            % signal to noise ratio
%  Pa1,             % a priori propability of the first bit
%  Pa2,             % a priori propability of the second bit
%  mapping,         % 1: Gray, 2: anti-Gray, 3: natural mapping
%  LLRampl,			  % maximum magnitude of the LLR-axis 
%  plotMaxLog,      % plot max-log approximation of the LLRs 
%  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, main function (control)

function [] = demo_LLR_4PAM( ...
    EsN0, ...
    Pa1, ...
    Pa2, ...
    mapping, ...
    LLRampl, ...
    plotMaxLog, ...
    svg_en, ...
    webdemo_en, ...
    filename )

% default parameters 
switch nargin
    case 0
        EsN0 = 5;
        Pa1 = 0.5;
        Pa2 = 0.5;
        mapping = 1;
        LLRampl = 15;
        plotMaxLog = 1;
        svg_en = 0;
        webdemo_en = 0;
        filename = 'demo_LLR_4PAM';
    otherwise
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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');
end

landscape = '-S930,350';
portrait  = '-S640,480';
output_format = portrait; %portrait;landscape; 

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;
end

% RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm
gray 	   = [[0,0,0]; [32,32,32]; [64,64,64]; [96,96,96]; [128,128,128]; [160,160,160]; [192,192,192]; [224,224,224]; [255,255,255]]/255;
violett    = [[51,0,25]; [102,0,51]; [153,0,76]; [204,0,102]; [255,0,127]; [255,51,153]; [255,102,178]; [255,153,204]; [255,204,229]]/255;
magenta    = [[51,0,51]; [102,0,102]; [153,0,153]; [204,0,204]; [255,0,255]; [255,51,255]; [255,102,255]; [255,153,255]; [255,204,255]]/255;
lilac      = [[25,0,51]; [51,0,102]; [76,0,153]; [102,0,204]; [127,0,255]; [153,51,255]; [178,102,255]; [204,153,255]; [229,204,255]]/255;
blue       = [[0,0,51]; [0,0,102]; [0,0,153]; [0,0,204]; [0,0,255]; [51,51,255]; [102,102,255]; [153,153,255]; [204,204,255]]/255;
seablue    = [[0,25,51]; [0,51,102]; [0,76,153]; [0,102,204]; [0,128,255]; [51,153,255]; [102,178,255]; [153,204,255]; [204,229,255]]/255;
cyan       = [[0,51,51]; [0,102,102]; [0,153,153]; [0,204,204]; [0,255,255]; [51,255,255]; [102,255,255]; [153,255,255]; [204,255,255]]/255;
seagreen   = [[0,51,25]; [0,102,51]; [0,153,76]; [0,204,102]; [0,255,128]; [51,255,153]; [102,255,178]; [153,255,204]; [204,255,229]]/255;
green      = [[0,51,0]; [0,102,0]; [0,153,0]; [0,204,0]; [0,255,0]; [51,255,51]; [102,255,102]; [153,255,153]; [204,255,204]]/255;
lightgreen = [[25,51,0]; [51,102,0]; [76,153,0]; [102,204,0]; [128,255,0]; [153,255,51]; [178,255,102]; [204,255,153]; [229,255,204]]/255;
yellow     = [[51,51,0]; [102,102,0]; [153,153,0]; [204,204,0]; [255,255,0]; [255,255,51]; [255,255,102]; [255,255,153]; [255,255,204]]/255;
amber      = [[51,25,0]; [102,51,0]; [153,76,0]; [204,102,0]; [255,128,0]; [255,153,51]; [255,178,102]; [255,204,153]; [255,229,204]]/255;
red        = [[51,0,0]; [102,0,0]; [153,0,0]; [204,0,0]; [255,0,0]; [255,51,51]; [255,102,102]; [255,153,153]; [255,204,204]]/255;


% END, set up graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN check input parameters for consistency


% END check input parameters for consistency
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Constants:
xmin = -5;
xmax = 5;
x_step = .1;
  
var = 1/(2*10^(EsN0/10));

% for numerical stability
switch Pa1
    case 1
        Pa1 = 1-eps;
    case 0
        Pa1 = eps;
end
switch Pa2
    case 1
        Pa2 = 1-eps;
    case 0
        Pa2 = eps;
end
% compute points
if( Pa1 >= 0 && Pa1 <=1 && Pa2 >= 0 && Pa2 <=1 )
    La1 = log(Pa1/(1-Pa1))
    La2 = log(Pa2/(1-Pa2))
    Pa1Vec = [1-Pa1 Pa1];
    Pa2Vec = [1-Pa2 Pa2];
    
    
    x_rx = (xmin:x_step:xmax);
    
    [m l] = get_mapping_and_labeling(mapping);
    Ps = sum(abs(l).^2 )/length(l);
    var = var * Ps;
    
    nb=2; % number of bits per symbol
    ns=2^nb; % number of symbols
    
    bit=zeros(ns,nb);
    for k=1:ns
        bit(k,:)=bitget(m(k),1:1:nb);    
    end
    
    y1_o = l(bit(:,1)==1);
    y1_u = l(bit(:,1)==0);
    y2_o = l(bit(:,2)==1);
    y2_u = l(bit(:,2)==0);
    
    c=-1/(2*var);
    thePDF = inline('exp(c*(arg).^2)', 'arg', 'c' );

    % Lch
    y_Lch_1 = log(  ( thePDF(x_rx-y1_o(1), c ) .* Pa2Vec(bit(find(l==y1_o(1)),2)+1) + thePDF(x_rx-y1_o(2), c) .* Pa2Vec( bit(find(l==y1_o(2)),2)+1) ) ./ ...
                    ( thePDF(x_rx-y1_u(1), c ) .* Pa2Vec(bit(find(l==y1_u(1)),2)+1) + thePDF(x_rx-y1_u(2), c) .* Pa2Vec( bit(find(l==y1_u(2)),2)+1) ) );
%     figure(2);
%     subplot(4,1,1);
%     plot(x_rx, thePDF(x_rx-y1_o(1), c),'r');
%     hold on;
%     plot(x_rx,thePDF(x_rx-y1_o(2), c) .* exp(La2),'g');
%     posLabeling = 0.1;
%     plotMapping(const_lw, bit, posLabeling);
%     hold off;
%     
%     subplot(4,1,2);
%     plot(x_rx,thePDF(x_rx-y1_u(1), c ),'r');
%     hold on;
%     plot(x_rx,thePDF(x_rx-y1_u(2), c ) .* exp(La2),'g' );
%     plotMapping(const_lw, bit, posLabeling);
%     hold off;

    y_L_1 = y_Lch_1 + La1;
    
    y_Lch_2 = log(  ( thePDF(x_rx-y2_o(2), c ) .* Pa1Vec(bit(find(l==y2_o(2)),1)+1) + thePDF(x_rx-y2_o(1), c) .* Pa1Vec(bit(find(l==y2_o(1)),1)+1) ) ./ ...
                    ( thePDF(x_rx-y2_u(1), c ) .* Pa1Vec(bit(find(l==y2_u(1)),1)+1) + thePDF(x_rx-y2_u(2), c) .* Pa1Vec(bit(find(l==y2_u(2)),1)+1) ) );
%     subplot(4,1,3);
%     plot(x_rx,thePDF(x_rx-y2_o(2) , c),'r');
%     hold on;
%     plot(x_rx,thePDF(x_rx-y2_o(1), c) .*  thePDF(La1, c) ,'g');
%     plotMapping(const_lw, bit, posLabeling);
%     hold off;
%
%     subplot(4,1,4);
%     plot(x_rx,thePDF(x_rx-y2_u(1) , c) ,'r');
%     hold on;
%     plot(x_rx,thePDF(x_rx-y2_u(2), c )  .* thePDF(La1, c),'g');
%     plotMapping(const_lw, bit, posLabeling);
%     hold off;

    y_L_2 = y_Lch_2 + La2;


    if( plotMaxLog )
        y_Lch_1_maxLog = max( c*(x_rx-y1_o(1)).^2 + La2*bit(find(l==y1_o(1)),2) , c*(x_rx-y1_o(2)).^2 + La2*bit(find(l==y1_o(2)),2) )...
                       - max( c*(x_rx-y1_u(1)).^2 + La2*bit(find(l==y1_u(1)),2) , c*(x_rx-y1_u(2)).^2 + La2*bit(find(l==y1_u(2)),2) );

        y_L_1_maxLog = y_Lch_1_maxLog + La1;

        y_Lch_2_maxLog = max( c*(x_rx-y2_o(1)).^2 + La1*bit(find(l==y2_o(2)),1) , c*(x_rx-y2_o(2)).^2 + La1*bit(find(l==y2_o(1)),1) ) ...
                       - max( c*(x_rx-y2_u(1)).^2 + La1*bit(find(l==y2_u(1)),1) , c*(x_rx-y2_u(2)).^2 + La1*bit(find(l==y2_u(2)),1) );

        y_L_2_maxLog = y_Lch_2_maxLog + La2;
    end

figure(1);
    plot(x_rx, y_Lch_1, '--', 'color', blue(2,:));
    hold on;
    plot(x_rx, y_L_1, 'color', blue(2,:));

    hold on;
    plot(x_rx, y_Lch_2, '--', 'color', green(2,:));
    hold on;
    plot(x_rx, y_L_2, 'color', green(2,:));
    hold on;

    if( plotMaxLog )
        plot(x_rx, y_L_1_maxLog, '-.', 'color', cyan(2,:));
        hold on;
        plot(x_rx, y_L_2_maxLog, '-.', 'color', lightgreen(2,:));
    end

    grid on;

    str = ['L_A(x_1) = ' num2str(La1) ', P_{apriori}[x_1=1] = ' num2str(Pa1), '  L_A(x_2) = ' num2str(La2) ', P_{apriori}[x_2=1] = ' num2str(Pa2) ];
    title(str);

    xlim([-5 5]);
    if(LLRampl~=0)
        ylim([-abs(LLRampl) abs(LLRampl)]);
    end
    ylimits = ylim;
    posLabeling = (abs(ylimits(1))+abs(ylimits(2)))/20;

    plotMapping(const_lw, bit, posLabeling);

    xlabel('y');
    ylabel('L-value');

    if (plotMaxLog)
        legend('L_{ch}(x_1|y)','L(x_1|y)','L_{ch}(x_2|y)','L(x_2|y)','max-log approx. L(x_1|y)','max-log approx. L(x_2|y)','4-PAM-Symbols','location', 'southeast');
    else
        legend('L_{ch}(x_1|y)','L(x_1|y)','L_{ch}(x_2|y)','L(x_2|y)', '4-PAM-Symbols','location', 'southeast');
    end
    
    hold off

else
    
    % Create an invisible figure.
    fh = figure(1);
    set(fh, 'visible', 'off');
    
    str = ['\nP_{apriori}[y=+1] has to be in the range [0 1]!'];
    text(0.1, 0.8, str, 'fontsize', 12, 'color', 'red'); 

end  
  
% Create an image file. This image is displayed on the webpage.
if (webdemo_en)
    print(filename, '-dpng',output_format);
end
if (svg_en)
    print(['fig_',filename], '-dsvg',output_format);
end

end

% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create TX signal
function [m l] = get_mapping_and_labeling(map)

if ( map == 1 ) % Gray
    m = [0, 1, 3, 2];
elseif ( map == 2 ) % anti-Gray
    m = [0, 3, 2, 1]; 
else % natural
    m = [0, 1, 2, 3]; 
end

l = [-3, -1, 1, 3];
%l = l ./(sum(abs(l).^2));

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function plotMapping(const_lw, bit, posLabeling)

% plot mapping
    plot((-3:2:3), zeros(1,4), '.', 'linewidth', const_lw, 'color', 'm', 'markersize',10);
    hold on;
    
    str=[num2str(bit(1,2)),num2str(bit(1,1))];
    text(-3,-posLabeling,str,'horizontalalignment','center','fontsize',6,'color','magenta');
    str=[num2str(bit(2,2)),num2str(bit(2,1))];
    text(-1,-posLabeling,str,'horizontalalignment','center','fontsize',6,'color','magenta');
    str=[num2str(bit(3,2)),num2str(bit(3,1))];
    text(+1,-posLabeling,str,'horizontalalignment','center','fontsize',6,'color','magenta');
    str=[num2str(bit(4,2)),num2str(bit(4,1))];
    text(+3,-posLabeling,str,'horizontalalignment','center','fontsize',6,'color','magenta');
 
end

% END, function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
