function demo_quantizer

	% numero di bit della legge di quantizzazione
	m = 4;
% 	M = 2^m; % non silenziata
	M = 2^m-1; % silenziata
	
	% caso
% 	demo = 'uniform';
% 	demo = 'exp';
% 	demo = 'laplace';
	demo = 'gaussian';
	
	switch demo
		case 'uniform'
			% fX(x) = rect(x/2)/2
			X_min = -1;
			X_max = +1;
			fX = @(x)uniform_pdf(x);
			Ex_fun = @(x)Ex_int_uniform(x);
			
		case 'exp'
			% fX(x) = exp(-x)*u(x)
			X_min = 0;
			X_max = 10;
			fX = @(x)exp_pdf(x);
			Ex_fun = @(x)Ex_int_exp(x);
			
		case 'laplace'
			% fX(x) = exp(-abs(x))/2 
			X_min = -10;
			X_max = +10;
			fX = @(x)laplace_pdf(x);
			Ex_fun = @(x)Ex_int_laplace(x);
			
		case 'gaussian'
			% fX(x) = exp(-x^2/2)/sqrt(2*pi) 
			X_min = -6;
			X_max = +6;
			fX = @(x)gaussian_pdf(x);
			Ex_fun = @(x)Ex_int_gaussian(x);
	end
	
	% inizializzazione delle variabili
	x = X_min+(X_max-X_min)/M*(0:M);
	first_loop = 1;
	
	% loop di aggiornamento iterativo delle soglie e dei livelli di
	% quantizzazione
	while first_loop || max(abs(x-x_old)) > 1e-8
		first_loop = 0;
		
		q = Ex_fun(x);
		x_old = x;
		x = [X_min,0.5*(q(1:end-1)+q(2:end)),X_max];
	end
	
	% plot della legge di quantizzazione
	figure(1)
	plot_Qlaw(x,q);
	
	% plot della densita' di probabilità
	figure(2)
	plot_fX(fX,X_min,X_max,q)
end

function plot_Qlaw(x,q)
	ex = 0.05*(x(end)-x(1));
	M = numel(q);
	
	hold off
	plot(x(round(0.5+(1:2*M)/2)),q(round((1:2*M)/2)))
	set(gca,'XLim',[x(1)-ex,x(end)+ex])
	xlabel('x')
	ylabel('Q(x)')
	
% 	lbl=dec2bin(0:M-1);
% 	for i = 1:M
% 		text(0.5*(x(i)+x(i+1)),q(i),lbl(i,:),...
% 			'FontSize',8,...
% 			'HorizontalAlignment','center',...
% 			'VerticalAlignment','baseline');%,'BackgroundColor',[1 1 1]);
% 	end
	
    hold on
    plot([x(1),x(end)],[x(1),x(end)],'r--')
    axis equal
end

function val = uniform_pdf(x)
	val = 0*x;
	val(abs(x)<=1) = 0.5;
end

function q = Ex_int_uniform(x)
	q = (x(2:end)+x(1:end-1))/2;
end

function val = exp_pdf(x)
	val = exp(-x);
	val(x<0) = 0;
end

function q = Ex_int_exp(x)
	p = exp(-x);
	xp = (x+1).*exp(-x);
	q = (xp(1:end-1)-xp(2:end))./(p(1:end-1)-p(2:end));
end

function val = laplace_pdf(x)
	val = exp(-abs(x))/2;
end

function q = Ex_int_laplace(x)
	M = numel(x)-1;
	q = zeros(1,M);
	for i = 1:M
		if x(i+1) < 0
			EX = (x(i+1)-1)*exp(x(i+1))-(x(i)-1)*exp(x(i));
			pX = (exp(x(i+1))-exp(x(i)));
		elseif x(i) > 0
			EX = (x(i)+1)*exp(-x(i))-(x(i+1)+1)*exp(-x(i+1));
			pX = (exp(-x(i))-exp(-x(i+1)));
		else
			EX = (1-x(i))*exp(x(i)) - (1+x(i+1))*exp(-x(i+1));
			pX = 2-exp(-x(i+1))-exp(x(i));
		end
		q(i) = EX/pX;
	end
end

function val = gaussian_pdf(x)
	val = exp(-x.^2/2)/sqrt(2*pi);
end

function q = Ex_int_gaussian(x)
	p = 0.5*erfc(x/sqrt(2));
	q = (exp(-x(1:end-1).^2/2)-exp(-x(2:end).^2/2)) ...
		./ (sqrt(2*pi)*(p(1:end-1)-p(2:end)));
end

function plot_fX(fX,X_min,X_max,q)
	ex = 0.05*(X_max-X_min);
	x = X_min-ex:(X_max-X_min)/1024:X_max+ex;
	hold off
	plot(x,fX(x),'r')
	hold on
	plot(q,fX(q),'rx')
	
	xlabel('x')
	ylabel('f_X(x)')
	set(gca,'XLim',[x(1),x(end)])
end
