function [IXY,IXY_LLR] = IXY_4PAM_AWGN_Soft(SNR_dB)

% 4-PAM con mapping di Gray:
%
%  -3d/2  -d/2    d/2   3d/2
% ---+------+------+------+---->
%    00     01     11     10  r
%            (b1,b2)


% in modo che la varianza dei simboli 
% informativi sia 1
d = sqrt(4/5);

% il rapporto S/N è 5*d^2/4*sigma_w^2
SNR = 10^(SNR_dB/10);
sigma_w = sqrt(5*d^2 / (4*SNR));

% punti di costellazione
V = [-3*d/2 -d/2 d/2 3*d/2];

% calcola l'integrale
S = -1;
N = 128;
while true
	dz = pi/N;
	z = ((1:N)/N-1/(2*N)-1/2)*pi;
	I = 0*z;
	for n = 1:N
		I(n) = integrando(sigma_w,V,z(n));
	end
	S_old = S;
	S = sum(I)*dz;
	if abs(S-S_old) < 1e-2
		break
	else
		N = N*2;
	end
end
IXY = S;

% numero di punti per la simulazione Monte-Carlo
% per il calcolo degli integrali
N = 8000;

% numero di segmenti in cui quantizzare
% l'intervallo (-pi/2,pi/2)
M = 32;

% estremi degli intervalli di quantizzazione
e = (0:M)/M*pi-pi/2;
% ey = a*tan(e);

% stima le probabilità di transizione di stato per i canali soft
% associati ai singoli bit e rivelazione soft con LLR
y_b1_0 = [V(1)*ones(1,N/2),V(2)*ones(1,N/2)]+sigma_w*randn(1,N);
L_b1_0 = LLR_b1(sigma_w,y_b1_0,V);
y_b1_1 = [V(3)*ones(1,N/2),V(4)*ones(1,N/2)]+sigma_w*randn(1,N);
L_b1_1 = LLR_b1(sigma_w,y_b1_1,V);

y_b2_0 = [V(1)*ones(1,N/2),V(4)*ones(1,N/2)]+sigma_w*randn(1,N);
L_b2_0 = LLR_b2(sigma_w,y_b2_0,V);
y_b2_1 = [V(2)*ones(1,N/2),V(3)*ones(1,N/2)]+sigma_w*randn(1,N);
L_b2_1 = LLR_b2(sigma_w,y_b2_1,V);

% questo fattore evita l'affollamento degli intervalli centrali
% e l'utilizzo degli intervalli più esterni
a = 2^ceil(log2(max(abs([L_b1_0,L_b1_1,L_b2_0,L_b2_1]))/tan(e(end-1))));

% ciascuno dei rivelatori soft LLR è equivalente
% ad un canale DMC con probabilità di transizione date da
p_b1_0 = histc(L_b1_0,a*tan(e))/N;
p_b1_1 = histc(L_b1_1,a*tan(e))/N;
T1 = [p_b1_0;p_b1_1].';

p_b2_0 = histc(L_b2_0,a*tan(e))/N;
p_b2_1 = histc(L_b2_1,a*tan(e))/N;
T2 = [p_b2_0;p_b2_1].';

% calcola l'informazione mutua dei canali con ingresso
% binario e uscita continua associati ai rapporti di 
% verosimiglianza dei bit
IXY_LLR = 0;
for n = 1:2
    for m = 1:M
        if T1(m,n)
		IXY_LLR = IXY_LLR + T1(m,n)/2 ...
			*log2(T1(m,n)/(1/2*sum(T1(m,:))));
        end
        if T2(m,n)
		IXY_LLR = IXY_LLR + T2(m,n)/2 ...
			*log2(T2(m,n)/(1/2*sum(T2(m,:))));
        end
    end
end

% \sum_n P_n*f_w(tan(z)-X(n))/cos^2(z) * ...
% log2( f_w(tan(z)-X(n))/(\sum_m f_w(tan(z)-X(n))*P_m) )
function val = integrando(sigma_w,X,z,p)
p = 0*X + 1/numel(X);
s = f_w(sigma_w,tan(z)-X);
if cos(z)==0
	val = 0;
else
	s_avg = sum(s .* p);
	val = 0;
	for n = 1:numel(X)
		if s(n)
			val = val + s(n) * p(n) * log2(s(n)/s_avg);
		end
	end
	val = val/cos(z)^2;
end

% densità di probabilità Gaussiana
function val = f_w(sigma_w,w)
val = exp(-w.^2/(2*sigma_w^2))/sqrt(2*pi*sigma_w^2);

% funzioni semplici per il calcolo degli LLR del primo e del secondo
% bit con la modulazione 4-PAM ed il mapping indicato all'inizio
% dello script
function val = LLR_b1(sigma_w,y,V)

val = log((f_w(sigma_w,y-V(1))+f_w(sigma_w,y-V(2)))...
	./(f_w(sigma_w,y-V(3))+f_w(sigma_w,y-V(4))));

for i = 1:numel(y)
	if val(i) > 36
		val(i) = 36;
	elseif val(i) < -36
		val(i) = -36;
	end
end

function val = LLR_b2(sigma_w,y,V)

val = log((f_w(sigma_w,y-V(1))+f_w(sigma_w,y-V(4)))...
	./(f_w(sigma_w,y-V(2))+f_w(sigma_w,y-V(3))));

for i = 1:numel(y)
	if val(i) > 36
		val(i) = 36;
	elseif val(i) < -36
		val(i) = -36;
	end
end
