% WEBDEMO demo_LLR_QAM.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: L_kaus Oestreich, oestreich@inue.uni-stuttgart.de
% date: 11.11.2013
%
% Written for GNU Octave (works for MATLAB as well)
%
% EXAMPLE:  demo_LLR_QAM(0, 0.5, 2, 1, 1, 1, 0, 0,'demo_LLR_QAM')
% 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
%  L_k,             % number of bit for which the L-value shall be plotted
%  keepRandomLabeling, % if selected, the current random labeling is stored for further investigations
%  LLRampl,			  % maximum magnitude of the LLR colormap
%  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_QAM( ...
    EsN0, ...
    Pa, ...
    nb, ...
    mapping, ...
    L_k, ...
    keepRandomLabeling, ...
    LLRampl, ...
    svg_en, ...
    webdemo_en, ... 
    filename )
  
% default parameters 
switch nargin
    case 0
        EsN0 = 24;
        Pa = [0.5 0.5 0.5 0.5];
        nb = 4;
        mapping = 3;
        L_k = 4;
        keepRandomLabeling = 1;
        LLRampl = 0;
        svg_en = 0;
        webdemo_en = 0;
        LLRampl = 400;
        filename = 'demo_LLR_QAM';
    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;

% still in testing mode: new color map, good for printing (white background), colorful like jet
nb_cols=256;
newmap= jet(nb_cols);
% generate smooth transition of jet() from "black to white" brightness
for (k=1:nb_cols)
newmap(k,1:3) = (newmap(k,1:3)*(k-1) + [nb_cols-k, nb_cols-k, nb_cols-k])/(nb_cols-1);
end
jet_light_inue=newmap;
colormap( jet_light_inue );

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

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

if( L_k > nb )
    L_k = nb;
end

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

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

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

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]';

% compute points
if( length(Pa(Pa >= 0)) == length(Pa) && length(Pa(Pa <= 1)) == length(Pa)  )
    La = log(Pa./(1-Pa));
    
    [x_rrx,x_irx] = meshgrid(xmin:x_step:xmax);
    x_rx = x_rrx + 1i * x_irx;
    
    xaxis = (xmin:x_step:xmax);
    yaxis = (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
    
    y_o = l(bit(:,L_k)==1);
    y_u = l(bit(:,L_k)==0);
    
    % Lch
    c=-1/(2*var);
    
    exp_o(1:(ns/2), 1:length(xaxis), 1:length(yaxis) ) = 0;
    exp_u(1:(ns/2), 1:length(xaxis), 1:length(yaxis) ) = 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
        
        % use max log approximation for SNR > 20
        if EsN0 > 20
            exp_o(n,:,:) =  c*abs(x_rx-y_o(n)).^2 + log(pOtherBitsO);
            exp_u(n,:,:) =  c*abs(x_rx-y_u(n)).^2 + log(pOtherBitsU); 
        else
            exp_o(n,:,:) =  exp(c*abs(x_rx-y_o(n)).^2) .* pOtherBitsO;
            exp_u(n,:,:) =  exp(c*abs(x_rx-y_u(n)).^2) .* pOtherBitsU;
        end
        
%         figure(2)
%         imagesc(exp_o);
%         figure(3)
%         imagesc(exp_u);
      
        
    end

    % use max log approximation for SNR > 20
    if EsN0 > 20
        y_Lch = squeeze(max(exp_o,[],1 ) - max(exp_u,[],1 ));
        y_L = y_Lch + La(L_k);
    else
        y_Lch = squeeze(log( sum(exp_o,1 )./ sum(exp_u,1 ) ) );
        y_L = y_Lch + La(L_k);
    end
    
    imagesc(xaxis,yaxis,(y_L));
    set(gca,'YDir','normal');   % By default, imagesc plots the y-axis from lowest to highest value, top to bottom. 
    
    if ( LLRampl~=0 )
        caxis([-LLRampl LLRampl]);
    end
    
    colorbar;
    hold on;
        
    % plot mapping
    plot(real(l), imag(l), '.', 'color', 'k', 'markersize',10);
    hold on;
    
    if nb == 1
      labelOffset = 0.25;
    else
      labelOffset = 0.5;
    end
      
    for k=1:ns
        str='';
        for n=1:nb
            str=[num2str(bit(k,n)),str];
        end
        text(real(l(k)),imag(l(k))-labelOffset,str,'horizontalalignment','center','fontsize',6,'color','black');
    end
    
    grid on;
    str = ['L-value of bit ' num2str(L_k,'%-.1f ') ', L_A(x_k) = ' num2str(La,'%-.1f ') ', P_{apriori}[x_k=-1] = ' num2str(1-Pa,'%-.1f ')];
    title(str);
    
    xlim([xmin xmax]);
    if nb == 1
        ylim([-0.5 0.5]);
    else
        ylim([xmin xmax]);
    end
    xlabel('Re(y)');
    ylabel('Im(y)');
    
    hold off
    
else
    % Create an invisible figure.
    fh = figure(1);
    set(fh, 'visible', 'off');
    
    str = 'P_{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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = [-1-1i, 1-1i, 1+1i, -1+1i];
elseif( nb == 3 )
    if ( map == 1 ) % Gray
        m = [0, 2, 6, 4, 1, 3, 7, 5];
    elseif ( map == 2 ) % anti-Gray
        m = [0, 6, 4, 2, 1, 7, 5, 3];
    elseif ( map == 3 ) % natural
        m = [0, 2, 4, 6, 1, 3, 5, 7];
    else % random
        if(newRandomLabeling)
            m = randperm(2^nb)-ones(1,2^nb);
            writeRandomLabelingInFile(nb,m);
        else
            m = readRandomLabelingFromFile(nb);
        end
    end
    l = [-1-1i,1-1i,1+1i,-1+1i, -(1+sqrt(3))*1i,1+sqrt(3),(1+sqrt(3))*1i,-1-sqrt(3)];
else
    if ( map == 1 ) % Gray
        m = [0, 1, 3, 2, 4, 5, 7, 6, 12, 13, 15, 14, 8, 9, 11, 10];
    elseif ( map == 2 ) % anti-Gray
        m = [0, 3, 2, 1, 12, 15, 14, 13, 8, 11, 10, 9, 4, 7, 6, 5];
    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 = [-3-3i,-1-3i,1-3i,3-3i, -3-1i,-1-1i,1-1i,3-1i, -3+1i,-1+1i,1+1i,3+1i, -3+3i,-1+3i,1+3i,3+3i];
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ...
function [xmin xmax] = get_xminmax(nb)

if( nb == 1 )
    xmin = -2;
    xmax = 2;
elseif( nb == 2 )
    xmin = -2;
    xmax = 2;
elseif( nb == 3 )
    xmin = -3.5;
    xmax = 3.5;
else
    xmin = -4;
    xmax = 4;
end

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


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

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

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


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

str_file = ['storeRandomLabeling',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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
