function [IXY,IXY_bit] = IXY_4PAM_AWGN_Hard(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];

% estremi sinistri e destri 
% delle regioni di decisione
L = [-inf  -d   0   d  ];
R = [  -d   0   d   inf];

% probabilità di transizione del canale
T = zeros(4,4);
for n = 1:4
    for m = 1:4
        T(m,n) = Q((L(m)+V(n))/sigma_w) ...
               - Q((R(m)+V(n))/sigma_w);
    end
end

% calcola l'informazione mutua
IXY = 0;
for n = 1:4
    for m = 1:4
        if T(m,n)
            IXY = IXY + T(m,n)/4*log2(T(m,n)/...
                (1/4*sum(T(m,:))));
        end
    end
end

% alternativamente, considerando un ricevitore
% Hard bit per bit, il modello diventa due canali
% paralleli con probabilità di transizione che si
% ottengono combinando le probabilità di 
% transizione delle regioni di decisione:
temp=[0.5*(T(:,1)+T(:,2)),0.5*(T(:,3)+T(:,4))];
T_b1=[temp(1,:)+temp(2,:);temp(3,:)+temp(4,:)];
temp=[0.5*(T(:,1)+T(:,4)),0.5*(T(:,2)+T(:,3))];
T_b2=[temp(1,:)+temp(4,:);temp(2,:)+temp(3,:)];

% calcola l'informazione mutua dei due canali bit per bit
IXY_bit = 0;
for n = 1:2
    for m = 1:2
        if T_b1(m,n)
			IXY_bit = IXY_bit + T_b1(m,n)/2 ...
				*log2(T_b1(m,n)/(1/2*sum(T_b1(m,:))));
        end
        if T_b2(m,n)
			IXY_bit = IXY_bit + T_b2(m,n)/2 ...
				*log2(T_b2(m,n)/(1/2*sum(T_b2(m,:))));
        end
    end
end

% Q(x) = Pr{N(0,1) > x}
function val = Q(x)
val = 0.5*erfc(x/sqrt(2));
