% WEBDEMO demo_ber_RdLimit.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: Klaus Oestreich, oestreich@inue.uni-stuttgart.de
% date: 13.12.2013
% update: 13.12.2014
%
% Written for GNU Octave (works for MATLAB as well)
%
% EXAMPLE:  demo_ber_RdLimit(2,1,1,0.5,1,0,3,-1,20,1,10000,2,0,0,0,0,0,'demo_ber')
% PARAMETERS:
%  nb,              % number of bits per symbol
%  labeling,        % 1: Gray, 2: anti-Gray, 3: natural labeling
%  mapPamQam,       % 1: QAM, 2: PAM, 3: PSK
%  Rc,              % Coderate for the R(D)-Limit
%  plotOverEbNo,    % plot the curves over Es/N0 or Eb/N0.
%  plotRdLimit,     % If true, the R(D)-Limit is calculated and plotted.
%  channel,         % 1: AWGN, 2: Rayleigh-Fading, 3: simulate both channels (AWGN and Rayleigh)
%	 EbNo_start,		  % Range of Eb/N0, first value
%	 EbNo_end,				% Range of Eb/N0, last value
%	 EbNo_inc,				% Range of Eb/N0, step size
%  nOfSymbols,      % number of symbols for simulation
%  algorithmSoftDemapping, % algorithmSoftDemapping = 1: probability domain, = 2: Log-Domain, Jacobian Log, = 3: Log-Domain, max-Log
%  hideLegend,			% 1: no legend
%  lowBER           % 1: down to 10^-5
%  keepPlot,        % if 1, then the results of the previous plot are shown in the following plot as well.
%  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_ber_RdLimit( ...
    nb, ...
    labeling, ...
    mapPamQam, ...
    Rc, ...
    plotOverEbNo, ...
    plotRdLimit, ...
    channel, ...
    EbNo_start, ...
    EbNo_end, ...
    EbNo_inc, ...
    nOfSymbols, ...
    algorithmSoftDemapping, ...
    hideLegend, ...
    lowBER, ...
    keepPlot, ...
    svg_en, ...
    webdemo_en, ...
    filename )

% default parameters
switch nargin
    case 0
        nb = 2;
        labeling = 1;
        mapPamQam = 1;
        Rc = 0.5;
        plotOverEbNo = 1;
        plotRdLimit = 0;
        channel = 3;
        EbNo_start = -1;
        EbNo_end = 20;
        EbNo_inc = 1;
        nOfSymbols = 10000;
        algorithmSoftDemapping = 2;
        hideLegend = 1;
        lowBER = 0;
        keepPlot = 0;
        svg_en = 0;
        webdemo_en = 0;
        filename = 'demo_ber';
    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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bitsPerSymbol = nb;
l_max =2^bitsPerSymbol-1;

labelingSequence = randi([0,l_max],1,nOfSymbols);

% get signal to transmit
x = mapping(bitsPerSymbol,labeling,mapPamQam,labelingSequence);

nOf_EbNo_values = floor((EbNo_end - EbNo_start)/EbNo_inc) + 1;

BER_h = ones(1,nOf_EbNo_values);
BER_s = ones(1,nOf_EbNo_values);
BER_h_ray = ones(1,nOf_EbNo_values);
BER_s_ray = ones(1,nOf_EbNo_values);
EbNo_dB_list=EbNo_start:EbNo_inc:EbNo_end;
EbNo_list = 10.^(EbNo_dB_list/10);

%%%%% RD
if (plotRdLimit)
    nOfBins = 300;
    I_s = zeros(1,nOf_EbNo_values);
    I_h = zeros(1,nOf_EbNo_values);
    I_s_ray = zeros(1,nOf_EbNo_values);
    I_h_ray = zeros(1,nOf_EbNo_values);
end
%%%%% RD

Es = 1;
Eb = Es/(bitsPerSymbol*Rc);

EsNo_dB_list =  EbNo_dB_list + 10*log10(bitsPerSymbol*Rc);
%%%%% RD
if (plotRdLimit)
    counter=1;
    for EsNo_dB = EsNo_dB_list;
        EsNo = 10^(EsNo_dB/10);
        % Noise variance
        No = Es/EsNo;
        n = sqrt(No/2)*(randn(1,nOfSymbols)+1i*randn(1,nOfSymbols));
        
        if(channel~=2)
            y = x + n;
            a = ones(1,nOfSymbols);
            
            rx_bits_hard = hardDemapping(bitsPerSymbol,labeling,mapPamQam,y);
            rx_bits_soft = softDemapping(bitsPerSymbol,labeling,mapPamQam,y,No,algorithmSoftDemapping,a);
        end
        if(channel~=1)
            % Rayleigh fading channel coefficients with zero mean and sigma^2 variance:
            a = 1/sqrt(2)*(randn(1,nOfSymbols)+ 1i*randn(1,nOfSymbols));
            y_ray = a.*x + n;
            % Optimal channel weights
            channel_weights = conj(a);
            % Rx. Signal Multiplied with the optimal weights of different channels
            y_ray = channel_weights.*y_ray;
            y_ray = y_ray./(channel_weights.*a);
            rx_bits_hard_ray = hardDemapping(bitsPerSymbol,labeling,mapPamQam,y_ray);
            rx_bits_soft_ray = softDemapping(bitsPerSymbol,labeling,mapPamQam,y_ray,No,algorithmSoftDemapping,a);
        end
        
        tx_bits = de2bi(labelingSequence);
        
        for m = 1:1:bitsPerSymbol;
            % AWGN:
            if(channel~=2)
                [hL_s, x_hl_s] = hist(rx_bits_soft(1:end,m),nOfBins);
                hL1_s = hist(rx_bits_soft(tx_bits(1:end,m)==1,m),x_hl_s);
                
                [hL_h, x_hl_h] = hist(rx_bits_hard(1:end,m),nOfBins);
                hL1_h = hist(rx_bits_hard(tx_bits(1:end,m)==1,m),x_hl_h);
                
                delta_x_s = x_hl_s(2) - x_hl_s(1);
                hL_s = hL_s / sum(hL_s) / delta_x_s;
                hL1_s = hL1_s / sum(hL1_s) / delta_x_s;
                
                delta_x_h = x_hl_h(2) - x_hl_h(1);
                hL_h = hL_h / sum(hL_h) / delta_x_h;
                hL1_h = hL1_h / sum(hL1_h) / delta_x_h;
                
                toSum_s = hL1_s*delta_x_s.*log2(hL1_s./hL_s);
                isn_s = isnan(toSum_s);
                teilErg_s =  sum( toSum_s(~isn_s));
                I_s(counter) = I_s(counter) + teilErg_s;
                
                toSum_h = hL1_h*delta_x_h.*log2(hL1_h./hL_h);
                isn_h = isnan(toSum_h);
                teilErg_h =  sum( toSum_h(~isn_h));
                I_h(counter) = I_h(counter) + teilErg_h;
            end
            
            % Rayleigh:
            if(channel~=1)
                [hL_s_ray, x_hl_s_ray] = hist(rx_bits_soft_ray(1:end,m),nOfBins);
                hL1_s_ray = hist(rx_bits_soft_ray(tx_bits(1:end,m)==1,m),x_hl_s_ray);
                
                [hL_h_ray, x_hl_h_ray] = hist(rx_bits_hard_ray(1:end,m),nOfBins);
                hL1_h_ray = hist(rx_bits_hard_ray(tx_bits(1:end,m)==1,m),x_hl_h_ray);
                
                delta_x_s_ray = x_hl_s_ray(2) - x_hl_s_ray(1);
                hL_s_ray = hL_s_ray / sum(hL_s_ray) / delta_x_s_ray;
                hL1_s_ray = hL1_s_ray / sum(hL1_s_ray) / delta_x_s_ray;
                
                delta_x_h_ray = x_hl_h_ray(2) - x_hl_h_ray(1);
                hL_h_ray = hL_h_ray / sum(hL_h_ray) / delta_x_h_ray;
                hL1_h_ray = hL1_h_ray / sum(hL1_h_ray) / delta_x_h_ray;
                
                toSum_s_ray = hL1_s_ray*delta_x_s_ray.*log2(hL1_s_ray./hL_s_ray);
                isn_s_ray = isnan(toSum_s_ray);
                teilErg_s_ray =  sum( toSum_s_ray(~isn_s_ray));
                I_s_ray(counter) = I_s_ray(counter) + teilErg_s_ray;
                
                toSum_h_ray = hL1_h_ray*delta_x_h_ray.*log2(hL1_h_ray./hL_h_ray);
                isn_h_ray = isnan(toSum_h_ray);
                teilErg_h_ray =  sum( toSum_h_ray(~isn_h_ray));
                I_h_ray(counter) = I_h_ray(counter) + teilErg_h_ray;
            end
            
        end
        counter=counter+1;
    end
end
%%%%% RD

counter=1;
for EbNo_dB = EbNo_dB_list;
    EbNo = 10.^(EbNo_dB/10);
    % Noise variance
    No   = Eb/EbNo;
    n = sqrt(No/2)*(randn(1,nOfSymbols)+1i*randn(1,nOfSymbols));
    tx_bits = de2bi(labelingSequence);
    
    if(channel~=2)
        y = x + n;
        
        a = ones(1,nOfSymbols);
        rx_bits_hard = hardDemapping(bitsPerSymbol,labeling,mapPamQam,y);
        rx_bits_soft = softDemapping(bitsPerSymbol,labeling,mapPamQam,y,No,algorithmSoftDemapping,a);
        % VZ von soft-Bits auswerten:
        rx_bits_soft = floor((sign(rx_bits_soft)+1)/2);
        B_h = sum(sum(abs(tx_bits - rx_bits_hard))) / (nOfSymbols*bitsPerSymbol);
        B_s = sum(sum(abs(tx_bits - rx_bits_soft))) / (nOfSymbols*bitsPerSymbol);
        BER_h(counter) = B_h;
        BER_s(counter) = B_s;
    end
    
    if(channel~=1)
        % Rayleigh fading channel coefficients with zero mean and sigma^2 variance:
        a = 1/sqrt(2)*(randn(1,nOfSymbols)+ 1i*randn(1,nOfSymbols));
        y_ray = a.*x + n;
        % Optimal channel weights
        channel_weights = conj(a);
        % Rx. Signal Multiplied with the optimal weights of different channels
        y_ray = channel_weights.*y_ray;
        y_ray = y_ray./(channel_weights.*a);
        rx_bits_hard_ray = hardDemapping(bitsPerSymbol,labeling,mapPamQam,y_ray);
        rx_bits_soft_ray = softDemapping(bitsPerSymbol,labeling,mapPamQam,y_ray,No,algorithmSoftDemapping,a);
        % VZ von soft-Bits auswerten:
        rx_bits_soft_ray = floor((sign(rx_bits_soft_ray)+1)/2);
        
        B_h_ray = sum(sum(abs(tx_bits - rx_bits_hard_ray))) / (nOfSymbols*bitsPerSymbol);
        B_s_ray = sum(sum(abs(tx_bits - rx_bits_soft_ray))) / (nOfSymbols*bitsPerSymbol);
        BER_h_ray(counter) = B_h_ray;
        BER_s_ray(counter) = B_s_ray;
    end
    
    counter=counter+1;
    
end

if (keepPlot)
    hold on;
else
    hold off;
end

if(channel~=2)
    M =2^bitsPerSymbol;
    if (mapPamQam == 1)
        if (bitsPerSymbol == 1)
            EbNo_calc = 0.5 * erfc( (1/sqrt(2)) * sqrt( 2*EbNo_list ) );
        elseif (bitsPerSymbol == 2)
            EbNo_calc = 0.5 * erfc( (1/sqrt(2)) * sqrt( 2*EbNo_list ) );
        else
            EbNo_calc = 4/(log2(M)) * 0.5 * erfc( (1/sqrt(2)) * sqrt( 3* EbNo_list * log2(M) / (M-1) ) );
        end
    elseif (mapPamQam == 2)
        EbNo_calc = 2/(log2(M)) * (1 - 1/M) * 0.5 * erfc( (1/sqrt(2)) * sqrt( 6 * log2(M) * EbNo_list / ( M*M - 1) ));
    else
        if (bitsPerSymbol == 1)
            EbNo_calc = 0.5 * erfc( (1/sqrt(2)) * sqrt( 2*EbNo_list ) );
        elseif (bitsPerSymbol == 2)
            EbNo_calc = 0.5 * erfc( (1/sqrt(2)) * sqrt( 2*EbNo_list ) );
        else
            EbNo_calc = 2/(log2(M)) * 0.5 * erfc( (1/sqrt(2)) * sqrt( 2* EbNo_list * log2(M) ) * sin(pi/M));
        end
    end
end



%%%%% RD
if (plotRdLimit)
    EbNo_i = EbNo_start:0.1:EbNo_end;
    REbNo_i = EbNo_start*Rc*bitsPerSymbol:0.1*Rc*bitsPerSymbol:EbNo_end*Rc*bitsPerSymbol;
    nOfRDvalues = length(REbNo_i);
    p_s = ones(1,nOfRDvalues)*1e-15;
    p_h = ones(1,nOfRDvalues)*1e-15;
    if(channel~=2)
        I_s_i = interp1(EbNo_dB_list*Rc*bitsPerSymbol,I_s,REbNo_i);
        I_h_i = interp1(EbNo_dB_list*Rc*bitsPerSymbol,I_h,REbNo_i);
    end
    
    p_s_ray = ones(1,nOfRDvalues)*1e-15;
    p_h_ray = ones(1,nOfRDvalues)*1e-15;
    if(channel~=1)
        I_s_i_ray = interp1(EbNo_dB_list*Rc*bitsPerSymbol,I_s_ray,REbNo_i);
        I_h_i_ray = interp1(EbNo_dB_list*Rc*bitsPerSymbol,I_h_ray,REbNo_i);
    end
    
    if(channel~=2)
        computeNextRd = 1;
        p0=[1e-5 0.5];
        n = 1;
        while computeNextRd
            fun = @(p) -p*log2(p) - (1-p)*log2(1-p) - (1 - I_s_i(n)/(nb * Rc));
            try
                p_res = fzero(fun,p0);
                if( isnan(p_res) )
                    computeNextRd = 0;
                else
                    %p0 = p_res*0.1;
                    p_s(n) = p_res;
                    if (p_res < 1e-4 || n==nOfRDvalues)
                        computeNextRd = 0;
                    else
                        n = n + 1;
                    end
                end
            catch
                computeNextRd = 0;
            end
        end
        
        computeNextRd = 1;
        p0=[1e-5 0.5];
        n_h = 1;
        while computeNextRd
            fun = @(p) -p*log2(p) - (1-p)*log2(1-p) - (1 - I_h_i(n_h)/(nb * Rc));
            try
                p_res = fzero(fun,p0);
                if( isnan(p_res) )
                    computeNextRd = 0;
                else
                    %p0 = p_res*0.1;
                    p_h(n_h) = p_res;
                    if (p_res < 1e-4 || n_h==nOfRDvalues)
                        computeNextRd = 0;
                    else
                        n_h = n_h + 1;
                    end
                end
            catch
                computeNextRd = 0;
            end
        end
        
    end
    
    
    %%%%%%%% Rayleigh
    if(channel~=1)
        
        computeNextRd = 1;
        p0=[1e-5 0.5];
        n_ray = 1;
        while computeNextRd
            fun = @(p) -p*log2(p) - (1-p)*log2(1-p) - (1 - I_s_i_ray(n_ray)/(nb * Rc));
            try
                p_res = fzero(fun,p0);
                if( isnan(p_res) )
                    computeNextRd = 0;
                else
                    %p0 = p_res*0.1;
                    p_s_ray(n_ray) = p_res;
                    if (p_res < 1e-4 || n_ray==nOfRDvalues)
                        computeNextRd = 0;
                    else
                        n_ray = n_ray + 1;
                    end
                end
            catch
                computeNextRd = 0;
            end
        end
        
        computeNextRd = 1;
        p0=[1e-5 0.5];
        n_h_ray = 1;
        while computeNextRd
            fun = @(p) -p*log2(p) - (1-p)*log2(1-p) - (1 - I_h_i_ray(n_h_ray)/(nb * Rc));
            try
                p_res = fzero(fun,p0);
                if( isnan(p_res) )
                    computeNextRd = 0;
                else
                    %p0 = p_res*0.1;
                    p_h_ray(n_h_ray) = p_res;
                    if (p_res < 1e-4 || n_h_ray==nOfRDvalues)
                        computeNextRd = 0;
                    else
                        n_h_ray = n_h_ray + 1;
                    end
                end
            catch
                computeNextRd = 0;
            end
        end
        
        %%%%%%%% Rayleigh
    end
    
end
%%%%% RD

shiftForEsNo = 0;
strXlabel = 'Eb/N0 [dB]';
if ( ~plotOverEbNo )
    shiftForEsNo = +10*log10(bitsPerSymbol*Rc);
    strXlabel = 'Es/N0 [dB]';
end

if(channel~=2)
    if (plotRdLimit)
        fill([EbNo_i(1:1) EbNo_i(1:n_h+1)]+shiftForEsNo,[1e-15 p_h(1:n_h+1)], lightgreen(9,:));
        hold on;
        h10 = fill([EbNo_i(1:1) EbNo_i(1:n+1)]+shiftForEsNo,[1e-15 p_s(1:n+1)], yellow(9,:));
        hold on;
        
        h6 = semilogy(EbNo_i(1:n_h)+shiftForEsNo,p_h(1:n_h), '-', 'color', green(3,:));
        hold on;
        h7 = semilogy(EbNo_i(1:n)+shiftForEsNo,p_s(1:n), '-', 'color', amber(3,:));
        hold on;
        %textUnachievable = text(-0.35+shiftForEsNo,1.1e-4,'unachievable region','horizontalalignment','left','VerticalAlignment','bottom','fontsize',10,'color','black');
        %set(textUnachievable, 'rotation', 90);
    end
    h1 = semilogy(EbNo_dB_list-10*log10(Rc)+shiftForEsNo, EbNo_calc, '-.', 'color', gray(5,:));
    hold on;
    
    h2 = semilogy(EbNo_dB_list+shiftForEsNo,BER_h, '+-',  'color', red(5,:));
    hold on;
    h3 = semilogy(EbNo_dB_list+shiftForEsNo,BER_s, '--', 'color', blue(5,:));
    hold on;
end

if(channel~=1)
    if (plotRdLimit)
        % Rayleigh RD-Limit:
        h8 = semilogy(EbNo_i(1:n_h_ray)+shiftForEsNo,p_h_ray(1:n_h_ray), '--', 'color', green(3,:));
        hold on;
        h9 = semilogy(EbNo_i(1:n_ray)+shiftForEsNo,p_s_ray(1:n_ray), '--', 'color', amber(3,:));
        hold on;
    end
    h4 = semilogy(EbNo_dB_list+shiftForEsNo,BER_h_ray, '+--',  'color', magenta(5,:));
    hold on;
    h5 = semilogy(EbNo_dB_list+shiftForEsNo,BER_s_ray, '--', 'color', cyan(5,:));
    hold on;
end


grid on;

%str = 'BER with hard- and soft-demapping';
%title(str);

xlim([EbNo_start+shiftForEsNo EbNo_end+shiftForEsNo]);
if (lowBER==1)
    ylim([1e-5 1]);
else
    ylim([1e-4 1]);
end
xlabel(strXlabel);
ylabel('BER');


if(hideLegend~=1)
    if(channel==1)
        if (plotRdLimit)
            legend([h1 h2 h3 h6 h7 h10],{'closed-form approx. (Gray, AWGN)','hard-demapping (AWGN)','soft-demapping (AWGN)', 'R(D)-limit hard-input FEC decoding (AWGN)', 'R(D)-limit soft-input FEC decoding (AWGN)', 'unachievable region (AWGN, soft-demapping)'},'location', 'northeast');
        else
            legend([h1 h2 h3],{'closed-form approx. (Gray, AWGN)','hard-demapping (AWGN)','soft-demapping (AWGN)'},'location', 'northeast');
        end
    elseif(channel==2)
        if (plotRdLimit)
            legend([h4 h5 h8 h9],{'hard-demapping (Rayleigh)','soft-demapping (Rayleigh)', 'R(D)-limit hard-input FEC decoding (Rayleigh)', 'R(D)-limit soft-input FEC decoding (Rayleigh)'},'location', 'northeast');
        else
            legend([h4 h5],{'hard-demapping (Rayleigh)','soft-demapping (Rayleigh)'}, 'location', 'northeast');
        end
    else
        if (plotRdLimit)
            legend([h1 h2 h3 h6 h7 h4 h5 h8 h9 h10],{'closed-form approx. (Gray, AWGN)','hard-demapping (AWGN)','soft-demapping (AWGN)', 'R(D)-limit hard-input FEC decoding (AWGN)', 'R(D)-limit soft-input FEC decoding (AWGN)','hard-demapping (Rayleigh)','soft-demapping (Rayleigh)', 'R(D)-limit hard-input FEC decoding (Rayleigh)', 'R(D)-limit soft-input FEC decoding (Rayleigh)', 'unachievable region (AWGN, soft-demapping)'},'location', 'northeast');
        else
            legend([h1 h2 h3 h4 h5],{'closed-form approx. (Gray, AWGN)','hard-demapping (AWGN)','soft-demapping (AWGN)','hard-demapping (Rayleigh)','soft-demapping (Rayleigh)'}, 'location', 'northeast');
        end
    end
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

function hardDemapper = hardDemapping(nb,labeling,mapPamQam,rxSig)

% smallest euclidian distance

[m, l] = get_mapping_and_labeling(nb,labeling,mapPamQam);

l_max=2^nb;
hardDecision = zeros(1,length(rxSig));

distance = zeros(l_max,length(rxSig));
for n = 1:l_max
    distance(n,:) = abs(rxSig-l(n));
end
[~,index_min_distance] = min(distance);
for n = 1:l_max
    hardDecision(index_min_distance==n) = m(n);
end

% same result but less sufficient:
%for a = 1:length(rxSig)
%dec=l(1);
%dist = realmax;
%for n = 1:l_max
%   distneu = abs(rxSig(a)-l(n));
%   if ( distneu < dist)
%       dec=m(n);
%       dist=distneu;
%   end
%end
%hardDecision(a)=dec;
%end

hardDemapper = de2bi(hardDecision);

end


function softDemapper = softDemapping(nb,labeling,mapPamQam,rxSig,No,algorithm,a)

% consider no a-piori knowledge
La = 0.0;

[m, l] = get_mapping_and_labeling(nb,labeling,mapPamQam);
ns=2^nb; % number of symbols

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

% compute L-values
softDemapper = zeros(length(rxSig),nb);

% kL: Nummer von L-Wert bei QPSK z.B. 2 L-Werte zu berechnen.
for kL=1:nb
    
    l_o = l(bit(:,kL)==1);
    l_u = l(bit(:,kL)==0);
    
    % Lch
    variance = No/2;
    c=-1./(2*variance*abs(1./(a.^2)));
    
    if(algorithm==1)
        exp_o = 0;
        exp_u = 0;
        for n=1:(ns/2)
            exp_o = exp_o + exp(c.*(abs(rxSig-l_o(n))).^2);
            exp_u = exp_u + exp(c.*(abs(rxSig-l_u(n))).^2);
        end
        y_Lch = log( exp_o ./ exp_u );
        
    else
        
        max_o = -1000;
        max_u = -1000;
        for n=1:(ns/2)
            if(algorithm==2)
                r_o = log(1 + exp(-abs( max_o - c.*(rxSig-l_o(n)).^2 )));
                r_u = log(1 + exp(-abs( max_u - c.*(rxSig-l_u(n)).^2 )));
            end
            max_o = max ( max_o, c.*abs(rxSig-l_o(n)).^2 );
            max_u = max ( max_u, c.*abs(rxSig-l_u(n)).^2 );
            if(algorithm==2)
                max_o = max_o + r_o;
                max_u = max_u + r_u;
            end
        end
        y_Lch = max_o - max_u;
    end
    
    y_L = y_Lch + La;
    
    softDemapper(1:end,kL) = y_L;
    
end

end



function mapArray =  mapping(nb,labeling,mapPamQam,labArray)

l_max=2^nb-1;
[m, l] = get_mapping_and_labeling(nb, labeling, mapPamQam);

mapArray = zeros(size(labArray));

for n = 1:l_max+1
    mapArray(labArray==m(n)) = l(n);
end


end

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

if( nb == 1 )
    m = [0, 1]; % Gray, anti-Gray, natural, random
    l = [-1, 1];
elseif( nb == 2 )
    
    mArray= [ 0, 1, 3, 2; % Gray, QAM
        0, 3, 2, 1; % anti-Gray
        0, 1, 2, 3;  % natural
        randperm(2^nb)-ones(1,2^nb); % random
        0, 1, 3, 2; % Gray, PAM
        0, 3, 2, 1; % anti-Gray
        0, 1, 2, 3;  % natural
        randperm(2^nb)-ones(1,2^nb); % random
        0, 1, 3, 2; % Gray, PSK
        0, 3, 2, 1; % anti-Gray
        0, 1, 2, 3;  % natural
        randperm(2^nb)-ones(1,2^nb)]; % random
    
    lArray =  [ (1/sqrt(2))*[-1-1i, 1-1i, 1+1i, -1+1i ]; ... % QAM
        (1/sqrt(5))*[-3+0i, -1+0i, 1+0i, 3+0]; ...  % PAM
        (1/sqrt(2))*[-1-1i, 1-1i, 1+1i, -1+1i]]; % PSK
    
    m = mArray(labeling+(mapPamQam-1)*4,1:end);
    l = lArray(mapPamQam,1:end);
    
elseif( nb == 3 )
    
    mArray = [ 0, 2, 6, 4, 1, 3, 7, 5; ...
        0, 6, 4, 2, 1, 7, 5, 3; ...
        0, 2, 4, 6, 1, 3, 5, 7;
        randperm(2^nb)-ones(1,2^nb);
        0, 1, 3, 2, 6, 7, 5, 4;
        0, 7, 1, 6, 3, 4, 2, 5;
        0, 1, 2, 3, 4, 5, 6, 7;
        randperm(2^nb)-ones(1,2^nb);
        0, 1, 3, 2, 7, 6, 5, 4;
        0, 7, 1, 6, 3, 4, 2, 5;
        0, 1, 2, 3, 4, 5, 6, 7;
        randperm(2^nb)-ones(1,2^nb)];
    
    lArray = [ (0.4597)*[-1-1i,1-1i,1+1i,-1+1i, -(1+sqrt(3))*1i,1+sqrt(3),(1+sqrt(3))*1i,-1-sqrt(3)];
        (1/sqrt(21))*[-7, -5, -3, -1, 1, 3, 5, 7];
        (1/sqrt(2))*(-1-1i), -1i, (1/sqrt(2))*(1-1i), 1, (1/sqrt(2))*(1+1i), 1i,(1/sqrt(2))*( -1+1i), -1 ];
    
    m = mArray(labeling+(mapPamQam-1)*4,1:end);
    l = lArray(mapPamQam,1:end);
    
else
    
    mArray = [ 0, 1, 3, 2, 4, 5, 7, 6, 12, 13, 15, 14, 8, 9, 11, 10;
        0, 3, 2, 1, 12, 15, 14, 13, 8, 11, 10, 9, 4, 7, 6, 5;
        0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15;
        randperm(2^nb)-ones(1,2^nb);
        0, 1, 3, 2, 6, 7, 5, 4, 12, 13, 15, 14, 10, 11, 9, 8;
        0, 15, 1, 14, 3, 12, 7, 8, 6, 9, 4, 11, 5, 10, 13, 2;
        0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15;
        randperm(2^nb)-ones(1,2^nb);
        0, 1, 3, 2, 6, 7, 5, 4, 12, 13, 15, 14, 10, 11, 9, 8;
        0, 15, 1, 14, 3, 12, 7, 8, 6, 9, 4, 11, 5, 10, 13, 2;
        0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15;
        randperm(2^nb)-ones(1,2^nb)];
    
    lArray =  [ (1/sqrt(10))*[-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];
        (1/sqrt(85))*[-15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15];
        exp(10i*pi/8), exp(11i*pi/8), -1i, exp(13i*pi/8), exp(14i*pi/8), exp(15i*pi/8), 1, exp(1i*pi/8), exp(2i*pi/8), exp(3i*pi/8), 1i, exp(5i*pi/8), exp(6i*pi/8), exp(7i*pi/8), -1, exp(9i*pi/8)];
    
    m = mArray(labeling+(mapPamQam-1)*4,1:end);
    l = lArray(mapPamQam,1:end);
    
end

end

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Copyright (C) 2001 Laurent Mazet
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% usage: b = de2bi(d, [n, [p]])
%%
%% Convert a non-negative integer to bit vector.
%%
%% d: positive integer
%% n: number of rows of the ouput bit vector (default is max. size).
%% p: base of decomposition (default is 2).
%%
%% b : bit vector.

%% 2001-02-02
%%   initial release

function b = de2bi(d, n, p)

if (nargin == 1)
    p = 2;
    n = floor ( log (max (max (d), 1)) ./ log (p) ) + 1;
elseif (nargin == 2)
    p = 2;
elseif (nargin ~= 3)
    error ('usage: b = de2bi (d [, n [, p]])');
end

d = d(:);
if ( any (d < 0) || any (d ~= floor (d)) )
    error ('de2bi: only handles non-negative integers');
end

power = ones (length (d), 1) * (p .^ (0 : n-1) );
d = d * ones (1, n);
b = floor (rem (d, p*power) ./ power);

end