% WEBDEMO demo_LLR_PaM.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_PaM(0, 0.5, 4, 1, 1, 15, 0, 1, 1,'demo_LLR_PaM')
% PaRAMETERS:
%  EsN0,            % signal to noise ratio
%  Pa,              % a priori propability
%  nb,              % number of bits per symbol
%  mapping,         % 1: Gray, 2: anti-Gray, 3: natural mapping
%  keepRandomLabeling, % if selected, the current random labeling is stored for further investigations
%  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_PaM( ...
    EsN0, ...
    Pa, ...
    nb, ...
    mapping, ...
    keepRandomLabeling, ...
    LLRampl, ...
    plotMaxLog, ...
    svg_en, ...
    webdemo_en, ...
    filename )

% default parameters
switch nargin
    case 0
        EsN0 = 0;
        Pa = [0.5 0];
        nb = 2;
        mapping = 1;
        keepRandomLabeling = 0;
        LLRampl = 15;
        plotMaxLog = 0;
        svg_en = 0;
        webdemo_en = 0;
        filename = 'demo_LLR_PaM';
    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 = landscape; %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;

myColors = {blue(3,:),green(3,:),amber(3,:),red(3,:),cyan(3,:),lilac(3,:),gray(1,:)};
myLinestyles = {'-','--','-.','.-',':'};


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

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


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

% Constants:
[xmin xmax] = get_xminmax(nb);
x_step = .1;

var = 1/(2*10^(EsN0/10));


if ischar(Pa) && length(find(Pa == '0' | Pa == '1' | Pa == '2' | Pa == '3' | Pa == '4' | Pa == '5' | Pa == '6' | Pa == '7' | Pa == '8' | Pa == '9'...
                             | Pa == '.' | Pa == ' ' | Pa == ']' | Pa == '[') ) == length(Pa)
    Pa = eval(Pa);
end

% compute points
if length(Pa) == 1
    Pa = ones(1, nb) .* Pa
end

% for numerical stability
Pa(Pa == 1) = 1-eps;
Pa(Pa == 0) = eps;

PaVec = [1-Pa; Pa]';

if( length(Pa(Pa >= 0)) == length(Pa) && length(Pa(Pa <= 1)) == length(Pa)  )
    La = log(Pa./(1-Pa));

    x_rx = (xmin:x_step:xmax);

    [m l] = get_mapping_and_labeling(nb,mapping,~keepRandomLabeling);
    Ps = sum(abs(l).^2)/length(l);
    var = var * Ps;

    ns=2^nb; % number of symbols

    bit=zeros(ns,nb);
    for k=1:ns
        bit(k,:)=bitget(m(k),1:1:nb);
    end

    % compute and plot LLR

    for k=1:nb
        y_o = l(bit(:,k)==1);
        y_u = l(bit(:,k)==0);

        % Lch:
        c=-1/(2*var);

        exp_o = 0;
        exp_u = 0;
        for n=1:(ns/2)

            pOtherBitsO = 1;
            pOtherBitsU = 1;
            for bitIndex = 1:nb
                if k ~= bitIndex
                    pOtherBitsO = pOtherBitsO * PaVec(bitIndex,bit(find(l==y_o(n)),bitIndex)+1);
                    pOtherBitsU = pOtherBitsU * PaVec(bitIndex,bit(find(l==y_u(n)),bitIndex)+1);
                end
            end


            exp_o = exp_o + exp(c*(x_rx-y_o(n)).^2) .* pOtherBitsO;
            exp_u = exp_u + exp(c*(x_rx-y_u(n)).^2) .* pOtherBitsU;
        end
        y_Lch = log( exp_o ./ exp_u );
        y_L = y_Lch + La(k);

        if( plotMaxLog )
            max_o = -1000;
            max_u = -1000;
            for n=1:(ns/2)
                max_o = max ( max_o, c*(x_rx-y_o(n)).^2 );
                max_u = max ( max_u, c*(x_rx-y_u(n)).^2 );
            end
            y_Lch_maxLog = max_o - max_u;
            y_L_maxLog = y_Lch_maxLog + La(k);
        end

        % plot(x_rx, y_Lch, '+-', 'linewidth', const_lw, 'color', 'g');
        % hold on;
        plot(x_rx, y_L, myLinestyles{mod(k-1,length(myLinestyles))+1}, 'color',myColors{mod(k-1,length(myColors))+1}, 'linewidth', const_lw, 'markersize', 3);
        hold on;


        if( plotMaxLog )
            plot(x_rx, y_L_maxLog, ':', 'color',myColors{mod(k-1,length(myColors))+1}, 'linewidth', const_lw, 'markersize', 3);
            hold on;
        end

    end

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

    % plot mapping
    plot(l, zeros(1,ns), '.', 'color', 'm', 'markersize',10);
    hold on;

    for k=1:ns
        str='';
        for n=1:nb
            str=[num2str(bit(k,n)),str];
        end
        text(l(k),-posLabeling,str,'horizontalalignment','center','fontsize',6,'color','magenta');
    end

    grid on;
    str = ['L_A(x_k) = [' num2str(La) '], P_{apriori}[x_k=-1] = [' num2str(1-Pa) ']'];
    title(str);

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

    plot_legend(nb,plotMaxLog)
    hold off

else

    % Create an invisible figure.
    fh = figure(1);
    set(fh, 'visible', 'off');

    str = ['P_{apriori}[x=+1] has to of length 1 or number of Bits/Symbol!' num2str(Pa)];
    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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select Mapping and labeling
function [m l] = get_mapping_and_labeling(nb,map,newRandomLabeling)

if( nb == 1 )
    if ( map ~= 4 )
        m = [0, 1]; % Gray, anti-Gray, natural
    else
        if(newRandomLabeling)
            m = randperm(2^nb)-ones(1,2^nb);
            writeRandomLabelingInFile(nb,m);
        else
            m = readRandomLabelingFromFile(nb);
        end
    end
    l = [-1, 1];
elseif( nb==2 )
    if ( map == 1 ) % Gray
        m = [0, 1, 3, 2];
    elseif ( map == 2 ) % anti-Gray
        m = [0, 3, 2, 1];
    elseif ( map == 3 ) % natural
        m = [0, 1, 2, 3];
    else % random
        if(newRandomLabeling)
            m = randperm(2^nb)-ones(1,2^nb);
            writeRandomLabelingInFile(nb,m);
        else
            m = readRandomLabelingFromFile(nb);
        end
    end
    l = [-3, -1, 1, 3];
elseif( nb == 3 )
    if ( map == 1 ) % Gray
        m = [0, 1, 3, 2, 6, 7, 5, 4];
    elseif ( map == 2 ) % anti-Gray
        m = [0, 7, 1, 6, 3, 4, 2, 5];
    elseif ( map == 3 ) % natural
        m = [0, 1, 2, 3, 4, 5, 6, 7];
    else % random
        if(newRandomLabeling)
            m = randperm(2^nb)-ones(1,2^nb);
            writeRandomLabelingInFile(nb,m);
        else
            m = readRandomLabelingFromFile(nb);
        end
    end
    l = [-7, -5, -3, -1, 1, 3, 5, 7];
else
    if ( map == 1 ) % Gray
        m = [0, 1, 3, 2, 6, 7, 5, 4, 12, 13, 15, 14, 10, 11, 9, 8];
    elseif ( map == 2 ) % anti-Gray
        m = [0, 15, 1, 14, 3, 12, 7, 8, 6, 9, 4, 11, 5, 10, 13, 2]; 
 	elseif ( map == 3 ) % natural
        m = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]; 
    else % random   
        if(newRandomLabeling)
            m = randperm(2^nb)-ones(1,2^nb);
            writeRandomLabelingInFile(nb,m);
        else
            m = readRandomLabelingFromFile(nb);
        end
    end
    l = [-15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15];

end 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select Mapping and labeling
function [xmin xmax] = get_xminmax(nb)

if( nb == 1 )
    xmin = -5;
    xmax = 5;
elseif( nb == 2 )
    xmin = -5;
    xmax = 5;
elseif( nb == 3 )
    xmin = -9;
    xmax = 9;
else
    xmin = -17;
    xmax = 17;
end 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot legend
function plot_legend(nb,plotMaxLog)

if(plotMaxLog == 0)
if( nb == 1 )
    legend('L(x|y)','2-PAM-Symbols','location', 'southeast');
elseif( nb == 2 )
    legend('L(x_1|y)','L(x_2|y)','4-PAM-Symbols','location', 'southeast');
elseif( nb == 3 )
    legend('L(x_1|y)','L(x_2|y)','L(x_3|y)','8-PAM-Symbols','location', 'southeast');
else
    legend('L(x_1|y)','L(x_2|y)','L(x_3|y)','L(x_4|y)','16-PAM-Symbols','location', 'southeast');
end

else
if( nb == 1 )
    legend('L(x|y)','max-log approx. L(x|y)','2-PAM-Symbols','location', 'southeast');
elseif( nb == 2 )
    legend('L(x_1|y)','max-log approx. L(x_1|y)','L(x_2|y)','max-log approx. L(x_2|y)','4-PAM-Symbols','location', 'southeast');
elseif( nb == 3 )
    legend('L(x_1|y)','max-log approx. L(x_1|y)','L(x_2|y)','max-log approx. L(x_2|y)','L(x_3|y)','max-log approx. L(x_3|y)','8-PAM-Symbols','location', 'southeast');
else
    legend('L(x_1|y)','max-log approx. L(x_1|y)','L(x_2|y)','max-log approx. L(x_2|y)','L(x_3|y)','max-log approx. L(x_3|y)','L(x_4|y)','max-log approx. L(x_4|y)','16-PAM-Symbols','location', 'southeast');
end    
end
 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ...
function [] = writeRandomLabelingInFile(nb,m)

str_file = ['storeRandomLabelingPaM',num2str(nb),'Bit.txt'];
fileID = fopen(str_file, 'w');
fwrite(fileID, m);
fclose(fileID);

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ...
function [m] = readRandomLabelingFromFile(nb)

str_file = ['storeRandomLabelingPaM',num2str(nb),'Bit.txt'];
[fileID,message] = fopen(str_file);
if(isempty(message))
    m = fread(fileID);
    fclose(fileID);
    if(length(m)~=2^nb) % something is wrong with the file
        m = randperm(2^nb)-ones(1,2^nb);
        writeRandomLabelingInFile(nb,m);
    end
else
    m = randperm(2^nb)-ones(1,2^nb);
    writeRandomLabelingInFile(nb,m);
end

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

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