% Filtro IIR passabanda (da Butterworth) di ordine N=M=6:   2/8<f/Fs>3/8
% coefficienti numeratore
b= [0.03 0 -0.09 0 0.09 0 -0.03];	% M+1=7 elementi
% coefficienti denominatore
a= [1 1.65 2.42 1.96 1.47 0.57 0.21];	% N+1=7 elementi

% visualizzo risposta in ampiezza e in fase del filtro
Npunti= 512;
Fs= 10000;	% [Hz] freq. di campionamento
figure(1); clf;
freqz(b,a,Npunti,Fs);
% visualizzo il ritardo di gruppo del filtro
figure(2); clf;
grpdelay(b,a,Npunti,Fs);
axis([0 Fs/2 0 10]);	% trascura la discontinuitā in f=0

% visualizzazione risposte impulsiva e indiciale (al gradino)
% Durata segnale
D= 50; % [# campioni]
% impulso di Dirac (di durata 50 campioni
imp= [1; zeros(D-1,1)];
% risposta impulsiva
h= filter(b,a,imp);
figure(3); clf;
% ascisse da 0 a D-1, altrimenti Matlab considera il vettore con indici da 1 a D
stem((0:D-1),h);
title('Risposta Impulsiva'); xlabel('n (# campione)'); ylabel('h(n)');
% gradino unitario
u= ones(D,1);
% risposta indiciale
g= filter(b,a,u);
figure(4); clf;
% ascisse da 0 a D-1, altrimenti Matlab considera il vettore con indici da 1 a D
stem((0:D-1),g);
title('Risposta Indiciale'); xlabel('n (# campione)'); ylabel('g(n)');

% plot zeri-poli
figure(5); clf;
zplane(b,a);
% calcolo esplicito zeri-poli
zer= roots(b);	% roots calcola le radici del polinomio b(1)z^M+...+b(M)*z+b(M+1)
pol= roots(a);	% qui č a(1)z^N+...+a(N)*z+a(N+1)
figure(6); clf;
% zer e pol sono vettori colonna: zplane se ne accorge e plotta direttamente gli elementi dei vettori
zplane(zer,pol);
% la fig. 6 dovrebbe coincidere con la 5

% calcolo funzione di trasferimento (vettori b e a) a partire da zeri e poli
% occorre il fattore di guadagno
k= b(1)/a(1);	% fattore di guadagno asintotico (per z^-1=0)
[b1, a1]= zp2tf(zer,pol,k);
figure(7); clf;
freqz(b1,a1,Npunti,Fs);
% da confrontare con fig.1

% metodo alternativo
% costruzione dei coefficienti del polinomio con radici assegnate:
b2= k*poly(zer);
a2= poly(pol);	% assegna di default a2(1)=1 (coefficiente di z^N
figure(8); clf;
freqz(b1,a1,Npunti,Fs);
% dovrebbe essere uguale alla fig. 1

% Tutte le grandezze viste sono riassunte da
fvtool(b,a);

% Forma in Serie

% Forma in parallelo