% tesine: generalizzare a due o più dimensioni
function IXY_DiscreteInputAWGN

	% costellazione
	X = [-3 -1 1 3];
	K = numel(X);
	
	% varianza del rumore gaussiano
	sigma2n = 0.1;
	
	% numero di intervalli da calcolare
	M = 4*K;
	c = max(X)+3*sqrt(sigma2n);
	
	clc
	while true

		% valori della funzione negli estremi degli M intervalli
		f = zeros(1,1+M);

		for m = 2:M
			% estremo destro dell'(m-1)-esimo intervallo
			z = -pi/2+(m-1)*pi/M;

			% sostituzione y -> c tan(z)
			y = c*tan(z);

			% denominatore dell'argomento del logaritmo
			val_K = 0;
			for l = 1:K
				val_K = val_K + ...
					1/K * fN(sigma2n, y - X(l));
			end

			% argomento dell'integrale
			val = 0;
			for i = 1:K
				val_N = fN(sigma2n, y - X(i));
				
				% si devono gestire i NaN nei 
				% singoli addendi e non alla fine!
				temp = val_N * log2(val_N / val_K);
				if isnan(temp)
					temp = 0;
				end

				val = val + 1/K * temp * c/cos(z)^2;
			end
			f(m) = val;
		end

		s = 0;
		S = 0;
		for m = 1:M
			if f(m) < f(m+1)
				s = s + f(m);
				S = S + f(m+1);
			else
				s = s + f(m+1);
				S = S + f(m);
			end
		end
		s = s * pi/M;
		S = S * pi/M;

		sigma2x = sum(X.^2)/K;
		C_AWGN = 0.5*log2(1+sigma2x/sigma2n);

		z = -pi/2:pi/M:pi/2;
		plot(z,f);
	
		if S/s > 1.1
			drawnow
			M = 2*M;
			fprintf('[%g %g %g] M %g->%g...\n',...
				s,(s+S)/2,S,M/2,M);
		else
			break
		end
	end
	fprintf('SNR: %g dB -- C_AWGN: %g -- IXY(est): %g [%g %g]\n',...
		10*log10(sigma2x/sigma2n),C_AWGN,(s+S)/2,s,S);
end

% densità di probabilità gaussiana
function val = fN(sigma2n,n)
	val = exp(-n^2/(2*sigma2n))/sqrt(2*pi*sigma2n);
end
