function image_filt2
	
	% ritardo massimo introdotto dal filtraggio inverso (per
	% immagini piccole o se si vogliono limitare gli effetti di
	% bordo questo deve valere poche unita'. se invece
	% interessa l'accuratezza di ricostruzione occorre
	% sceglierlo pari ad alcune volte il massimo ritardo
	% introdotto dal filtro di distorsione)
% 	D = 4;
	D = 100;
	
	% taglia i transitori dall'immagine filtrata
	do_TRUNC = false;
% 	do_TRUNC = true;
	
	% quantizza l'immagine filtrata
	do_QUANT = false;
% 	do_QUANT = true;

	% pseudoinversa con o senza regolarizzazione
	try_REGULARIZATION = false;
% 	try_REGULARIZATION = true;
	REGULARIZATION_LAMBDA = 0.0015;

	% risposta all'impulso del filtro
	L = 14;
	h = exp(-5*(-L/2:L/2).^2/(L/2)^2); % Gaussiano
% 	h = 1.05.^(0:L); % esponenziale crescente (troncato)
% 	h = 0.9.^(0:L); % esponenziale decrescente (troncato)
% 	h = [1./0.8.^(-floor(L/2):0),0.8.^(1:L-floor(L/2))];
% 	h = ones(1,1+L);
	h = h/sum(h);
	if do_TRUNC
		h_delay = round(sum(abs(h).^2.*(0:L))...
			/sum(abs(h).^2));
	end
	
	% visualizza la risposta del filtro
	figure(2)
	subplot(2,2,1)
	stem(-3:L+3,[0 0 0 h 0 0 0])
	title('risposta all''impulso del filtro')
	
	% carica l'immagine (a livelli di grigio)
% 	a = double(imread('imm1.png'));
% 	a = double(imread('imm2.png'));
 	a = double(imread('lena.png'));
%  	a = double(imread('audi.png'));
	if size(a,3) > 1
		a = 0.299*a(:,:,1)+0.114*a(:,:,2)+0.587*a(:,:,3);
	end
	[nr,nc]=size(a);
	
	% mostra l'immagine originale
	figure(1)
	subplot(2,2,1)
	imagesc(a,[0 255])
	colormap gray
	axis image
	title('immagine originale')
	
	% filtra l'immagine
	a_blur = conv2(conv2(a,h.'),h);
	
	% troncamento ed errore di quantizzazione
	if do_TRUNC
		a_blur = a_blur(h_delay+(1:nr),h_delay+(1:nc));
	end
	if do_QUANT
		a_blur = double(uint8(a_blur));
	end
	
	% mostra l'immagine distorta
	figure(1)
	subplot(2,2,2)
	imagesc(a_blur,[0 255])
	colormap gray
	axis image
	switch do_TRUNC+2*do_QUANT
		case 0
			title('immagine filtrata (non troncata ne'' quantizzata)')
		case 1
			title('immagine filtrata e troncata (non quantizzata)')
		case 2
			title('immagine filtrata e quantizzata (non troncata)')
		case 3
			title('immagine filtrata, troncata e quantizzata')
	end
	
	% costruisci la matrice di Toeplitz che rappresenta il
	% filtraggio
	H = toeplitz([h(end);zeros(D,1)],...
		[h(end:-1:1),zeros(1,D)]);
	
	% pseudoinversa (eventualmente regolarizzata)
	if try_REGULARIZATION
		W = (H'*H+REGULARIZATION_LAMBDA*eye(L+1+D))\H';
	else
		W = pinv(H);
	end
	
	% trova l'offset corrispondente alla distanza minima
	% da una delta di Kronecker
	err = zeros(1,L+1+D);
	for i=0:L+D
		d = conv(W(1+i,end:-1:1),h);
		d0 = d(L+1+D-i);
		d(L+1+D-i) = 0;
		err(1+i) = sum(abs(d).^2)/d0^2;
	end
	
	% visualizza la metrica di costo in funzione del ritardo
	% del filtro
	figure(2)
	subplot(2,2,2)
	semilogy(0:L+D,err)
	title('costo in funzione del ritardo complessivo')
	
	% scegli il ritardo in base alla metrica scelta
	[ignored,best_i]=min(err); best_i = best_i - 1;
	w = W(1+best_i,end:-1:1);
	w = w/sum(w);
	
	% visualizza la risposta del filtro inverso
	figure(2)
	subplot(2,2,3)
	stem(w)
	title('risposta del filtro inverso')
	
	% visualizza il prodotto di convoluzione tra filtro
	% e filtro inverso
	figure(2)
	subplot(2,2,4)
	d = zeros(1,L+D+1);
	d(1+L+D-best_i) = 1;
	hold off
	stem(0:L+D,d,'r:')
	hold on
	stem(-3:L+D+3,[0 0 0 conv(w,h) 0 0 0])
	title('risposta della cascata dei filtri')
	
	% calcola la posizione dell'immagine ricostruita
	if do_TRUNC
		w_delay = L+D-best_i-h_delay;
	else
		w_delay = L+D-best_i;
	end
	% equalizza l'immagine
	a_sharp = conv2(conv2(a_blur,w.'),w);
	if w_delay < 0
		a_sharp = [zeros(-w_delay+size(a_sharp,1),-w_delay),...
			[zeros(-w_delay,size(a_sharp,2));a_sharp]];
		w_delay = 0;
	elseif nr+w_delay > size(a_sharp,1) || ...
			nc+w_delay > size(a_sharp,2)
		a_sharp = [[a_sharp;...
			zeros(nr+w_delay-size(a_sharp,1),size(a_sharp,2))],...
			zeros(nr+w_delay,nc+w_delay-size(a_sharp,1))];
	end
	a_sharp = uint8(a_sharp(w_delay+(1:nr),w_delay+(1:nc)));
	
	% visualizza l'immagine equalizzata
	figure(1)
	subplot(2,2,3)
	imagesc(a_sharp,[0 255])
	colormap gray
	axis image
	title('immagine ricostruita')
	
	% visualizza l'errore di ricostruzione
	subplot(2,2,4)
	imagesc(abs(a-double(a_sharp)),[0 127])
	colormap gray
	axis image
	title('modulo dell''errore di ricostruzione')
	
end
