function IXY_demo

	% rapporto segnale-rumore
	SNR_dB = -10:0.01:10;
	
	SNR = 10.^(SNR_dB/10);
	p = Q(sqrt(SNR));
	I_BSC = 1-H2(p);
	
	C_AWGN = 0.5*log2(1+SNR);
	
	V = [1 -1];
	dR = 0.1;
	I_BPSK = zeros(1,numel(SNR));
	for k = 1:numel(SNR)
		sigma = 1/sqrt(SNR(k));
		
		R = [-inf,(-1:dR:1)*(1+4*sigma),inf];
		pY_V = zeros(2,numel(R)-1);
		for i = 1:2
			for j = 1:numel(R)-1
				pY_V(i,j) = Q((R(j)-V(i))/sigma) - ...
					Q((R(j+1)-V(i))/sigma);
			end
		end
		pY = 0.5*pY_V(1,:) + 0.5*pY_V(2,:);
		temp = 0;
		for i = 1:2
			for j = 1:numel(R)-1
				temp = temp + 0.5*pY_V(i,j) * ...
					log2(eps+pY_V(i,j)/pY(j));
			end
		end
		I_BPSK(k) = temp;
		
	end
	
	hold off
	plot(SNR_dB,I_BSC);
	hold on
	plot(SNR_dB,C_AWGN,'r')
	plot(SNR_dB,I_BPSK,'m')
	xlabel('SNR[dB]');
	ylabel('I(X;Y)');
	title('Informazione mutua per BPSK su AWGN');
	legend({'d;d_{est}','C_{AWGN}','v;r'},...
		'Location','SouthEast')

end

% calcola l'entropia di una sorgente binaria con
% probabilità p ed 1-p
function val = H2(p)
	val = -p.*log2(p)-(1-p).*log2(1-p);
	val(p == 0 | p == 1) = 0;
end

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