% 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: 14.01.2014
% update: 13.12.2014
%
% Written for GNU Octave (works for MATLAB as well)
%
% EXAMPLE:  demo_symbolErrorRate( 1, 1, 3, -1, 30, 2, 10000, 0, 0, 0, 0, 'demo_symbolErrorRate')
% PARAMETERS:
%  nb,              % number of bits per symbol
%  mapPamQam,       % 1: QAM, 2: PAM, 3: PSK
%  channel,			    % 1: AWGN, 2: Rayleigh-Fading, 3: simulate both channels (AWGN and Rayleigh)
%	 EsNo_start,		  % Range of Es/N0, first value
%	 EsNo_end,				% Range of Es/N0, last value
%	 EsNo_inc,				% Range of Es/N0, step size
%  nOfSymbols,		  % number of symbols for simulation
%  hideLegend,			% 1: no legend
%  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_symbolErrorRate( ...
    nb, ...
    mapPamQam, ...
    channel, ...
    EsNo_start, ...
    EsNo_end, ...
    EsNo_inc, ...
    nOfSymbols, ...
    hideLegend, ...
    keepPlot, ...
    svg_en, ...
    webdemo_en, ...
    filename )

% default parameters
switch nargin
    case 0
        nb = 1;
        mapPamQam = 1;
        channel = 3;
        EsNo_start = -1;
        EsNo_end = 30;
        EsNo_inc = 2;
        nOfSymbols = 10000;
        hideLegend = 1;
        keepPlot = 0;
        svg_en = 0;
        webdemo_en = 0;
        filename = 'demo_symbolErrorRate';
    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
labeling = 1;
x = mapping(bitsPerSymbol,labeling,mapPamQam,labelingSequence);

nOf_EsNo_values = floor((EsNo_end - EsNo_start)/EsNo_inc) + 1;

BER_h = ones(1,nOf_EsNo_values);
BER_h_ray = ones(1,nOf_EsNo_values);
EsNo_dB_list = EsNo_start:EsNo_inc:EsNo_end;
EsNo_list = 10.^(EsNo_dB_list/10);
Es = 1; % normalized constellation

% channel: add AWGN
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;
        rx_symbols_hard = hardSymbolDemapping(bitsPerSymbol,labeling,mapPamQam,y);
        
        B_h = sum((labelingSequence - rx_symbols_hard)~=0) / nOfSymbols;
        BER_h(counter) = B_h;
        
    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_symbols_hard_ray = hardSymbolDemapping(bitsPerSymbol,labeling,mapPamQam,y_ray);
        
        B_h_ray = sum((labelingSequence - rx_symbols_hard_ray)~=0) / nOfSymbols;
        BER_h_ray(counter) = B_h_ray;
    end
    
    counter=counter+1;
    
end

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

M =2^bitsPerSymbol;
if(channel~=2)
    if (mapPamQam == 1) % QAM
        if (bitsPerSymbol == 1)
            EsNo_calc = 0.5 * erfc( (1/sqrt(2)) * sqrt( 2*EsNo_list ) );
        elseif (bitsPerSymbol == 2)
            EsNo_calc = 2 * 0.5 * erfc( (1/sqrt(2)) * sqrt( EsNo_list ) );
        else
            EsNo_calc = 3 * 0.5 * erfc( (1/sqrt(2)) * sqrt( 3* EsNo_list / (M-1) ) );
        end
    elseif (mapPamQam == 2) % PAM
        EsNo_calc = 2 * (1 - 1/M) * 0.5 * erfc( (1/sqrt(2)) * sqrt( 6 * EsNo_list / ( M*M - 1) ));
    else % PSK
        if (bitsPerSymbol == 1)
            EsNo_calc = 0.5 * erfc( (1/sqrt(2)) * sqrt( 2*EsNo_list ) );
        elseif (bitsPerSymbol == 2)
            EsNo_calc = 2 * 0.5 * erfc( (1/sqrt(2)) * sqrt( EsNo_list ) );
        else
            EsNo_calc = 2 * 0.5 * erfc( (1/sqrt(2)) * sqrt( 2* EsNo_list ) * sin(pi/M));
        end
    end
end

if(channel~=2)
    semilogy(EsNo_dB_list, EsNo_calc, '-.', 'color', gray(5,:));
    hold on;
    
    semilogy(EsNo_dB_list, BER_h, '+-',  'color', red(5,:));
    hold on;
end
if(channel~=1)
    semilogy(EsNo_dB_list, BER_h_ray, 'o-',  'color', magenta(5,:));
    hold on;
end

grid on;

xlim([EsNo_start EsNo_end]);
ylim([1e-4 1]);
xlabel('Es/N0 [dB]');
ylabel('SER');

if(hideLegend~=1)
    if(channel==1)
        legend('closed-form approx. (AWGN)','hard-demapping (AWGN)','location', 'northeast');
    elseif(channel==2)
        legend('hard-demapping (Rayleigh)','location', 'northeast');
    else
        legend('closed-form approx. (AWGN)','hard-demapping (AWGN)', 'hard-demapping (Rayleigh)', 'location', 'northeast');
    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 hardSymbolDemapper = hardSymbolDemapping(nb,labeling,mapPamQam,rxSig)

% smallest euclidian distance

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

l_max=2^nb;

hardSymbolDemapper = 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
    hardSymbolDemapper(index_min_distance==n) = m(n);
end


% same result but less sufficient:
%for a = 1:length(rxSig)
%dec=m(1);
%dist = realmax;
%for n = 1:l_max
%   distneu = abs(rxSig(a)-l(n));
%   if ( distneu < dist)
%       dec=m(n);
%       dist=distneu;
%   end
%end
%hardSymbolDemapper(a)=dec;
%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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
