% WEBDEMO demo_mutual.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
% modified by: Marvin Geiselhart, geiselhart@inue.uni-stuttgart.de
% update: 05.08.2020 (replaced distrubution estimation by direct formula for Mutual Information)
%
% Written for GNU Octave (works for MATLAB as well)
%
% EXAMPLE:  demo_mutual(1,1,1,3,-10,20,1,0,6,0,10000,2,0,0,0,0,0,'demo_mutual')
% PARAMETERS:
%  nb,              % number of bits per symbol
%  labeling,        % 1: Gray, 2: anti-Gray, 3: natural mapping
%  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
%  mutInf_start,    % range of axis "mutual information"
%  mutInf_end,      % range of axis "mutual information"
%  plotOverEbNo,    % plot the curves over Es/N0 or Eb/N0.
%  nOfSymbols,		  % number of symbols for simulation
%  algorithmSoftDemapping, % algorithmSoftDemapping = 1: probability domain, = 2: Log-Domain, Jacobian Log, = 3: Log-Domain, max-Log
%  cumsum_mutualInf, % ...
%  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_mutual( ...
	nb, ...
	labeling, ...
	mapPamQam, ...
	channel, ...
	EsNo_start, ...
	EsNo_end, ...
	EsNo_inc, ...
	mutInf_start, ...
	mutInf_end, ...
	plotOverEbNo, ...
	nOfSymbols, ...
	algorithmSoftDemapping, ...
	cumsum_mutualInf, ...
	hideLegend, ...
	keepPlot, ...
	svg_en, ...
	webdemo_en, ...
	filename )

% default parameters
switch nargin
case 0
	nb = 2;
	labeling = 1;
	mapPamQam = 1;
	channel = 3;
	EsNo_start = -10;
	EsNo_end = 20;
	EsNo_inc = 1;
	mutInf_start = 0;
	mutInf_end = 6;
	plotOverEbNo = 0;
	nOfSymbols = 10000;
	algorithmSoftDemapping = 2;
	cumsum_mutualInf = 0;
	hideLegend = 1;
	keepPlot = 0;
	svg_en = 0;
	webdemo_en = 0;
	filename = 'demo_mutual';
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_EsNo_values = floor((EsNo_end - EsNo_start)/EsNo_inc) + 1;
nOfBins = 300;

I_s = zeros(bitsPerSymbol,nOf_EsNo_values);
I_h = zeros(bitsPerSymbol,nOf_EsNo_values);
I_s_ray = zeros(bitsPerSymbol,nOf_EsNo_values);
I_h_ray = zeros(bitsPerSymbol,nOf_EsNo_values);
C_ray = zeros(1,nOf_EsNo_values);

EsNo_dB_list = EsNo_start:EsNo_inc:EsNo_end;
EsNo_list = 10.^(EsNo_dB_list/10);
Es = 1; % normalized constellation

EbNo_start = -50;	% for capacity over EbNo
EbNo_end = 25;
EbNo_dB_list = EbNo_start:EsNo_inc:EbNo_end;
EbNo_list = 10.^(EbNo_dB_list/10);

% channel: add AWGN
counter=1;

% hideLegend: hideLegend hideHardMutialInf 
% 0: 0 0
% 1: 1 0
% 2: 0 1
% 3: 1 1
hideHardMutialInf = 0
if ( hideLegend == 2 )
	hideLegend = 0;
	hideHardMutialInf = 1;
elseif ( hideLegend == 3 )
	hideLegend = 1;
	hideHardMutialInf = 1;
end

noise_re = randn(1,nOfSymbols);
noise_im = randn(1,nOfSymbols);

noise_re = noise_re / std(noise_re);
noise_im = noise_im / std(noise_im);

noise_normal = noise_re + 1j*noise_im;

for EsNo_dB = EsNo_dB_list;
	EsNo = 10^(EsNo_dB/10);
    % Noise variance
    No   = Es/EsNo;
    n = sqrt(No/2)*noise_normal;
    
    if(channel~=2)
    	y = x + n;

    	a = ones(1,nOfSymbols);
    	if ( ~hideHardMutialInf )
    		rx_bits_hard = hardDemapping(bitsPerSymbol,labeling,mapPamQam,y);
    	end
    	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);
        
        if ( ~hideHardMutialInf )
        	rx_bits_hard_ray = hardDemapping(bitsPerSymbol,labeling,mapPamQam,y_ray);
        end
        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);
            hL0_s = hist(rx_bits_soft(tx_bits(1:end,m)==0,m),x_hl_s);
            
            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;
            hL0_s = hL0_s / sum(hL0_s) / delta_x_s;
            
            toSum_s_1 = hL1_s*delta_x_s.*log2(hL1_s./hL_s); 	
            toSum_s_1 = 0.5*hL1_s*delta_x_s.*log2(2*hL1_s./(hL1_s + hL0_s)); 
	    	toSum_s_0 = 0.5*hL0_s*delta_x_s.*log2(2*hL0_s./(hL1_s + hL0_s));
	    	isn_s = isnan(toSum_s_1);
            teilErg_s_1 =  sum( toSum_s_1(~isn_s));
	    	isn_s = isnan(toSum_s_0);
            teilErg_s_0 =  sum( toSum_s_0(~isn_s));
            I_s(m,counter) = I_s(m,counter) + teilErg_s_1 + teilErg_s_0;
	    	soft_mi = mean(1-log2(1+exp(-rx_bits_soft(1:end,m).*(2*tx_bits(1:end,m)-1)))); % Cammerer Formula
            %I_s(m,counter) = I_s(m,counter) + soft_mi; %teilErg_s_1 + teilErg_s_0;
            
            if ( ~hideHardMutialInf )
            % [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_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_h = hL1_h*delta_x_h.*log2(hL1_h./hL_h);
            % isn_h = isnan(toSum_h);
            % teilErg_h =  sum( toSum_h(~isn_h));

            pe = mean(rx_bits_hard(1:end,m) != tx_bits(1:end,m));
            hard_mi = binary_mi(pe);

            I_h(m,counter) = I_h(m,counter) + hard_mi; %teilErg_h;
        end
    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);
            
            %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;
            
            %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));
            soft_mi = mean(1-log2(1+exp(-rx_bits_soft_ray(1:end,m).*(2*tx_bits(1:end,m)-1)))); % Cammerer Formula
            I_s_ray(m,counter) = I_s_ray(m,counter) + soft_mi; %teilErg_s_ray;
            
            if ( ~hideHardMutialInf )
            % [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_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_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));

            pe = mean(rx_bits_hard_ray(1:end,m) != tx_bits(1:end,m));
            hard_mi = binary_mi(pe);

            I_h_ray(m,counter) = I_h_ray(m,counter) + hard_mi; % teilErg_h_ray;
        end
    end
end

counter=counter+1;

end

if(cumsum_mutualInf != 2)
	if(channel~=2)
		I_s = cumsum(I_s,1);
		if ( ~hideHardMutialInf )
			I_h = cumsum(I_h,1);
		end
	end
	if(channel~=1)
		I_s_ray = cumsum(I_s_ray,1);
		if ( ~hideHardMutialInf )
			I_h_ray = cumsum(I_h_ray,1);
		end
	end
end

% Capacity: AWGN
if(channel~=2)
	C = log2( 1 + EsNo_list );
	C_EbNo = log2( 1 + EbNo_list );
end
% Capacity: Rayleigh
if(channel~=1)
	C_ray_dB_list = zeros(1,nOf_EsNo_values);
	if (plotOverEbNo)
		C_ray_dB_list = EbNo_dB_list;
	else
		C_ray_dB_list = EsNo_dB_list;
	end

	counter=1;
	for EsNo_dB = C_ray_dB_list;
		EsNo = 10^(EsNo_dB/10);

        % Rayleigh fading channel coefficients with zero mean and sigma^2 variance:
        a = 1/sqrt(2)*(randn(1,nOfSymbols)+ 1i*randn(1,nOfSymbols));
        C_ray(counter) = sum(log2( 1 + abs(a).^2*EsNo))/nOfSymbols;
        
        counter=counter+1;
    end
end

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

if(plotOverEbNo)
	if(channel~=1)
		unreachable_X_ray = [-2, -2, C_ray_dB_list(1)-10*log10(C_ray(1)),C_ray_dB_list-10*log10(C_ray)];
		unreachable_Y_ray = [C_ray(end), 0, 0, C_ray];
		fill(unreachable_X_ray, unreachable_Y_ray, seablue(8,:) ); hold on
	end
	if(channel~=2)
		unreachable_X = [-2, -2, EbNo_dB_list(1)-10*log10(C_EbNo(1)), EbNo_dB_list-10*log10(C_EbNo)];
		unreachable_Y = [C_EbNo(end), 0, 0, C_EbNo];
		fill(unreachable_X, unreachable_Y, seablue(9,:) ); hold on
	end

	x_lim = 0;

	if(channel~=1)
		h1 = plot( EbNo_dB_list-10*log10(C_ray), C_ray, '--', 'color', seablue(1,:) ); hold on
		if ( cumsum_mutualInf > 0 )
			if ( ~hideHardMutialInf )
				h2 = plot( (repmat(EsNo_dB_list,bitsPerSymbol,1)-10*log10(I_h_ray))', I_h_ray', '--', 'color', red(1,:) ); hold on
			end
			h3 = plot( (repmat(EsNo_dB_list,bitsPerSymbol,1)-10*log10(I_s_ray))', I_s_ray', '--', 'color', green(1,:) ); hold on
		else
			if ( ~hideHardMutialInf )
				h2 = plot( EsNo_dB_list-10*log10(I_h_ray(bitsPerSymbol,:)), I_h_ray(bitsPerSymbol,:), '--', 'color', red(1,:) ); hold on
			end
			h3 = plot( EsNo_dB_list-10*log10(I_s_ray(bitsPerSymbol,:)), I_s_ray(bitsPerSymbol,:), '--', 'color', green(1,:) ); hold on
		end
		x_lim = EsNo_dB_list(end)-10*log10(I_s_ray(end));
	end
	if(channel~=2)
		h4 = plot( EbNo_dB_list-10*log10(C_EbNo), C_EbNo, 'color', seablue(1,:) ); hold on
		if ( cumsum_mutualInf > 0 )
			if ( ~hideHardMutialInf )
				h5 = plot( (repmat(EsNo_dB_list,bitsPerSymbol,1)-10*log10(I_h))', I_h', '-', 'color', red(2,:) ); hold on
			end
			h6 = plot( (repmat(EsNo_dB_list,bitsPerSymbol,1)-10*log10(I_s))', I_s', '-', 'color', green(2,:) ); hold on
		else
			if ( ~hideHardMutialInf )
				h5 = plot( EsNo_dB_list-10*log10(I_h(bitsPerSymbol,:)), I_h(bitsPerSymbol,:), '-', 'color', red(2,:) ); hold on
			end
			h6 = plot( EsNo_dB_list-10*log10(I_s(bitsPerSymbol,:)), I_s(bitsPerSymbol,:), '-', 'color', green(2,:) ); hold on
		end
		x_lim = EsNo_dB_list(end)-10*log10(I_s(end));
	end

	if( (EsNo_end-EsNo_start) == 30 && hideLegend~=1)
		textUnachievable = text(2,2.7,'unachievable region','horizontalalignment','center','fontsize',10,'color','black');
		set(textUnachievable, 'rotation', 55)
	end

	axis([-2 x_lim mutInf_start mutInf_end]);
	xlabel('Eb/N0 [dB]');

else
	if(channel~=1)
		unreachable_X_ray = [EsNo_dB_list(1), EsNo_dB_list, EsNo_dB_list(end)];
		unreachable_Y_ray = [C_ray(end)+6, C_ray, C_ray(end)+6];
		fill(unreachable_X_ray, unreachable_Y_ray, seablue(9,:) ); hold on
	end
	if(channel~=2)
		unreachable_X = [EsNo_dB_list(1), EsNo_dB_list, EsNo_dB_list(end)];
		unreachable_Y = [C(end)+6, C, C(end)+6];
		fill(unreachable_X, unreachable_Y, seablue(8,:) ); hold on
	end

	if(channel~=1)
		h1 = plot( EsNo_dB_list, C_ray, '--', 'color', seablue(1,:) ); hold on
		if ( cumsum_mutualInf > 0 )
			if ( ~hideHardMutialInf )
				h2 = plot( EsNo_dB_list, I_h_ray', '--', 'color', red(1,:) ); hold on
			end
			h3 = plot( EsNo_dB_list, I_s_ray', '--', 'color', green(1,:) ); hold on
		else
			if ( ~hideHardMutialInf )
				h2 = plot( EsNo_dB_list, I_h_ray(bitsPerSymbol,:), '--', 'color', red(1,:) ); hold on
			end
			h3 = plot( EsNo_dB_list, I_s_ray(bitsPerSymbol,:), '--', 'color', green(1,:) ); hold on
		end
	end
	if(channel~=2)
		h4 = plot( EsNo_dB_list, C, 'color', seablue(2,:) ); hold on
		if ( cumsum_mutualInf > 0 )
			if ( ~hideHardMutialInf )
				h5 = plot( EsNo_dB_list, I_h', '-', 'color', red(2,:) ); hold on
			end
			h6 = plot( EsNo_dB_list, I_s', '-', 'color', green(2,:) ); hold on
		else
			if ( ~hideHardMutialInf )
				h5 = plot( EsNo_dB_list, I_h(bitsPerSymbol,:), '-', 'color', red(2,:) ); hold on
			end
			h6 = plot( EsNo_dB_list, I_s(bitsPerSymbol,:), '-', 'color', green(2,:) ); hold on
		end
	end

	if( (EsNo_end-EsNo_start) == 30 && hideLegend~=1)
		textUnachievable = text(5,3,'unachievable region','horizontalalignment','center','fontsize',10,'color','black');
		set(textUnachievable, 'rotation', 45)
	end

	axis([EsNo_start EsNo_end mutInf_start mutInf_end]);
	xlabel('Es/N0 [dB]');
end

grid on;

%str = 'mutual information';
%title(str);

ylabel('mutual information [bit/channel use]');

if(hideLegend~=1)
	if(channel==1)
		if ( hideHardMutialInf )
			legend([h4 h6(1)],{'Shannon capacity limit (AWGN)','soft-demapping (AWGN)'},'location', 'northwest');
		else
			legend([h4 h5(1) h6(1)],{'Shannon capacity limit (AWGN)','hard-demapping (AWGN)','soft-demapping (AWGN)'},'location', 'northwest');
		end
	elseif(channel==2)
		if ( hideHardMutialInf )
			legend([h1 h3(1)],{'Shannon capacity limit (Rayleigh)','soft-demapping (Rayleigh)'},'location', 'northwest');
		else
			legend([h1 h2(1) h3(1)],{'Shannon capacity limit (Rayleigh)','hard-demapping (Rayleigh)','soft-demapping (Rayleigh)'},'location', 'northwest');
		end
	else
		if ( hideHardMutialInf )
			legend([h1 h3(1) h4 h6(1)],{'Shannon capacity limit (Rayleigh)','soft-demapping (Rayleigh)','Shannon capacity limit (AWGN)','soft-demapping (AWGN)'},'location', 'northwest');
		else
			legend([h1 h2(1) h3(1) h4 h5(1) h6(1)],{'Shannon capacity limit (Rayleigh)','hard-demapping (Rayleigh)','soft-demapping (Rayleigh)','Shannon capacity limit (AWGN)','hard-demapping (AWGN)','soft-demapping (AWGN)'},'location', 'northwest');
		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 = labArray;

	for n = 1:l_max+1
		mapArray(mapArray==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= [ get_labeling(nb,1,1); % Gray, QAM
             get_labeling(nb,0,1); % anti-Gray
                     0:3;  % natural
                     randperm(2^nb)-ones(1,2^nb); % random
                     get_labeling(nb,1,0); % Gray, PAM
                     get_labeling(nb,0,0); % anti-Gray
                     0:3;  % natural
                     randperm(2^nb)-ones(1,2^nb); % random
                     get_labeling(nb,1,0); % Gray, PSK
                     get_labeling(nb,0,0); % anti-Gray
                     0: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, 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)];

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

                elseif( nb == 4 )
                	mArray = [ get_labeling(nb,1,1); 
                	get_labeling(nb,0,1); 
                	0:15;
                	randperm(2^nb)-ones(1,2^nb);
                          get_labeling(nb,1,0); %0, 1, 3, 2, 6, 7, 5, 4, 12, 13, 15, 14, 10, 11, 9, 8; 
                          get_labeling(nb,0,0); 
                          0:15;
                          randperm(2^nb)-ones(1,2^nb);
                          get_labeling(nb,1,0); 
                          get_labeling(nb,0,0); 
                          0: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);

                      elseif( nb == 6 )

                      	mArray = [ get_labeling(nb,1,1);
                      	get_labeling(nb,0,1); 
                      	0:63;
                      	randperm(2^nb)-ones(1,2^nb);
                      	get_labeling(nb,1,0);
                      	get_labeling(nb,0,0);
                      	0:63;
                      	randperm(2^nb)-ones(1,2^nb); 
                      	get_labeling(nb,1,0);
                      	get_labeling(nb,0,0);
                      	0:63;
                      	randperm(2^nb)-ones(1,2^nb)];

                      	lArray =  [ (1/sqrt(42))*[ reshape(ones(8)*(-9-9i)+cumsum(ones(8)*2i,2)+cumsum(ones(8)*2),1,[]) ]; 
                      	(1/sqrt(1365))*[-63:2:63];
                      	cumprod(ones(1,64)*exp(i*pi/32)) ];

                      	m = mArray(labeling+(mapPamQam-1)*4,1:end);
                      	l = lArray(mapPamQam,1:end); 

 else % nb == 8

 	mArray = [ get_labeling(nb,1,1);
 	get_labeling(nb,0,1);
 	0:255;
 	randperm(2^nb)-ones(1,2^nb);
 	get_labeling(nb,1,0);
 	get_labeling(nb,0,0);
 	0:255;
 	randperm(2^nb)-ones(1,2^nb);
 	get_labeling(nb,1,0);
 	get_labeling(nb,0,0);
 	0:255;
 	randperm(2^nb)-ones(1,2^nb)];

 	lArray =  [ (1/sqrt(170))*[ reshape(ones(16)*(-17-17i)+cumsum(ones(16)*2i,2)+cumsum(ones(16)*2),1,[]) ]; 
 	(1/sqrt(21845))*[-255:2:255];
 	cumprod(ones(1,256)*exp(i*pi/128)) ];

 	m = mArray(labeling+(mapPamQam-1)*4,1:end);
 	l = lArray(mapPamQam,1:end);
 end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get labeling
function labelArray = get_labeling(nb,grey,qam)

	nb_codes = 2^nb;
	Gray_tab = zeros(nb_codes,nb);
	Gray_tab(1,:)=-1;

	if (grey == 1)
		for i = 0:nb-1
			for j = 0:(2^i)-1
				for s = 0:nb-1
					Gray_tab(2^(i+1)-1-j+1,s+1) = Gray_tab(j+1,s+1);
				end
				Gray_tab(2^(i+1)-1-j+1,i+1) = 1;
			end
		end

	else
% anti-grey
for i = 0:nb-2
	for j = 0:(2^i)-1
		for s = 0:nb-2
			Gray_tab(2^(i+1)-1-j+1,s+1) = Gray_tab(j+1,s+1);
		end
		Gray_tab(2^(i+1)-1-j+1,nb) = -1;
		Gray_tab(2^(i+1)-1-j+1,i+1) = 1;
	end
end

% second: rearrange and insert inverted versions (see Hamming-book)
for i = nb_codes/2-1:-1:0
	for s = 0:nb-1
		Gray_tab(2*i+2,s+1) = -Gray_tab(i+1,s+1);
		Gray_tab(2*i+1,s+1) = Gray_tab(i+1,s+1);
	end
end

end

Gray_tab(Gray_tab==-1)=0;

Gray_dec = zeros(nb_codes,nb);
for i=1:nb
	Gray_dec(:,i)=2^(i-1);
end
Gray_dec = cumsum(Gray_dec.*Gray_tab,2);
labelArray = Gray_dec(:,nb)';

if(qam==1&&nb>2)
	labelArray = reshape(labelArray,sqrt(nb_codes),[])';
	labelArray(2:2:end) = fliplr(labelArray(2:2:end,:));
	labelArray = labelArray'(:)';
end

end


function mi = binary_mi(pe)
	mi = (1-pe)*log2(2*(1-pe));
	if (pe>0)
		mi = mi + pe*log2(2*pe);
	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
