function iir_demo
	
	type = 'butter';
% 	type = 'cheby';

	transform = 'imp_inv';
% 	transform = 'bilin';
	
	% frequenza di taglio
	fc = 4000;
	
	% ordine del filtro (di Butterworth o di Chebychev)
	N = 13;
	% tolleranza per il filtro di Chebychev
	tol_dB=-3;
	e = sqrt(10^(-tol_dB/10)-1);
	
	% frequenza di campionamento
	fs = 44100;
	
	% --------------------------------------------------------------------------
	
	switch type
		case 'butter'
			% poli del filtro prototipo di Butterworth
			p = 2*pi*fc*exp(1i*((pi+mod(N,2)*pi)/(2*N)+(0:2*N-1)*pi/N));

			% per la trasformazione bilineare, a causa della distorsione
			% dell'asse delle frequenze, occorre calcolare un altro polinomio
			if strcmp(transform,'bilin')
				fc_bilin = fs/(pi)*tan(pi*fc/fs);
				p_bilin = 2*pi*fc_bilin*exp(1i*((pi+mod(N,2)*pi)/(2*N)+...
					(0:2*N-1)*pi/N));
			end
			
		case 'cheby'
			
			% polinomio di Chebychev
			% T0(x) = 1
			T_prev = [zeros(1,N),1];
			% T1(x) = x/(1i*omegac)
			T      = [zeros(1,N-1),1/(1i*2*pi*fc),0];
			for i = 1:N-1
				temp = T;
				% TN(x) = 2*x*T{N-1}(x) - T{N-2}(x)
				T = 2/(1i*2*pi*fc)*[temp(2:end),0] - T_prev;
				T_prev = temp;
			end
			poly = e^2*conv(T,T);
			poly(end) = poly(end)+1;
			p = roots(poly);
			
			% per la trasformazione bilineare, a causa della distrosione
			% dell'asse delle frequenze, occorre calcolare un altro polinomio
			if strcmp(transform,'bilin')
				fc_bilin = fs/(pi)*tan(pi*fc/fs);
				% T0(x) = 1
				T_prev = [zeros(1,N),1];
				% T1(x) = x/(1i*omegac)
				T      = [zeros(1,N-1),1/(1i*2*pi*fc_bilin),0];
				for i = 1:N-1
					temp = T;
					% TN(x) = 2*x*T{N-1}(x) - T{N-2}(x)
					T = 2/(1i*2*pi*fc_bilin)*[temp(2:end),0] - T_prev;
					T_prev = temp;
				end
				poly = e^2*conv(T,T);
				poly(end) = poly(end)+1;
				p_bilin = roots(poly);
			end
	end
	
	figure(1)
	subplot(4,1,1)
	hold off
	plot(real(p),imag(p),'rx','MarkerSize',7)
	hold on
	plot([0,0],2.4*pi*[-fc,fc],'r--')
	axis equal
	
	% seleziona i poli a parte reale negativa
	ps = p(real(p)<0);
	if strcmp(transform,'bilin')
		ps_bilin = p_bilin(real(p_bilin)<0);
	end
	
	% denominatore della funzione di trasferimento in s
	temp = 1;
	for n = 1:N
		temp = conv(temp,[1,-1/ps(n)]);
	end
	Ds = temp;
	if strcmp(type,'cheby') && mod(N,2)==0
		Ds = Ds*sqrt(1+e^2);
	end
	[omega_s,Hs] = freq_plot(1,Ds,fs/2);

	if strcmp(transform,'bilin')
		temp = 1;
		for n = 1:N
			temp = conv(temp,[1,-1/ps_bilin(n)]);
		end
		Ds_bilin = temp;
		if strcmp(type,'cheby') && mod(N,2)==0
			Ds_bilin = Ds_bilin*sqrt(1+e^2);
		end
	end
	
	switch transform
		case 'imp_inv'
			% posizione dei poli a tempo discreto
			pz = exp(ps/fs);
			zz = [];
	
			% denominatore della funzione di trasferimento in z
			temp = 1;
			for n = 1:N
				temp = conv(temp,[1,-pz(n)]);
			end
			Dz = real(temp);
			if strcmp(type,'cheby') && mod(N,2)==0
				Dz = Dz*sqrt(1+e^2);
			end
			Nz = real(prod(1-pz));
			[omega_z,Hz] = freq_plot(Nz,Dz);
		case 'bilin'
			p1=[1 +1];
			m1=[1 -1];
			base=zeros(1+N);
			base(1,end)=1;
			for n = 1:N
				base(1+n,end-n:end)=conv(base(n,end-n+1:end),p1);
			end
			temp = m1;
			for n = N-1:-1:0
				base(n+1,:)=conv(base(n+1,end-n:end),temp);
				temp=conv(temp,m1);
			end
			% base(1,:) == (1-z^-1)^N*(1+z^-1)^0
			% base(2,:) == (1-z^-1)^(N-1)*(1+z^-1)
			% ...
			% base(N+1,:) == (1-z^-1)^0*(1+z^-1)^N
			Dz = zeros(1,1+N);
			for i = 0:N
				Dz = Dz + Ds_bilin(1+i)*(2*fs)^i * base(N+1-i,:);
			end
			Nz = base(1+N,:);
			[omega_z,Hz] = freq_plot(Nz,Dz);
			pz = roots(Dz);
			zz = roots(Nz);
	end
	
	subplot(4,1,2)
	hold off
	plot(real(pz),imag(pz),'x','MarkerSize',7)
	hold on
	if numel(zz)>0
		plot(real(zz),imag(zz),'o','MarkerSize',7)
	end
	plot(exp(1i*(-pi:pi/32:pi)),'--');
	axis equal
	
	subplot(4,1,3)
	hold off
	plot(omega_z/pi*fs/2,20*log10(abs(Hz)))
	hold on
	plot(omega_s/(2*pi),20*log10(abs(Hs)),'r')
	set(gca,'YLim',[20*log10(min(abs(Hs)))-5,5]);
	set(gca,'XLim',[-fs/2,fs/2]);
	grid on
	
	subplot(4,1,4)
	hold off
	temp = unwrap(angle(Hz));
	plot(omega_z/pi*fs/2,temp-temp(numel(temp)/2+1))
	hold on
	temp = unwrap(angle(Hs));
	plot(omega_s/(2*pi),temp-temp(numel(temp)/2+1),'r')
	set(gca,'XLim',[-fs/2,fs/2]);
	grid on
	
% 	if ~exist('a.wav','file')
% 		system(['mpg123 /home/ftp/audio/mp3/Muse/Origin\ of\ '...
% 			'Simmetry/09.\ Dark\ shines.mp3 -w a.wav']);
% 		x = wavread('a.wav',44100*30);
% 		wavwrite(x,fs,16,'a.wav');
% 	else
% 		x = wavread('a.wav',44100*30);
% 	end
% % 	system('play a.wav');
% 	figure(2)
% 	subplot(2,1,1)
% 	pwelch(x(:,1));
% 	y = filter(real(Nz),real(Dz),x);
% 	subplot(2,1,2)
% 	pwelch(y(:,1));
% 	wavwrite(y,fs,16,'afilt.wav');
% 	system('play afilt.wav');

    M = fs/10;
	ts = (0:M-1)/fs;
	x = 2*(cos(2*pi*50*ts)>0)-1;
	y = filter(real(Nz),real(Dz),x);
	
	figure(2)
	hold off
	plot(ts,x,'--')
	hold on
	plot(ts,y,'r')

end
