% Progetto di un filtro di IIR (Butterworth) col Metodo dell'Invarianza Impulsiva
% partiamo da specifiche assegnate direttamente sulla frequenza normalizzata

clear all;
% Specifiche
deltap= 1-0.89125;		% 1=>|H|=>0.89125
Attp=  -20*log10(1-deltap);	% massima attenuazione
% OCCHIO: è diverso dal parametro Ap= 10*log10(1+deltap) che si usa di solito
fp= 0.2;	% *Fs/2		(per       0<=|f|<=0.1*Fs)
deltas= 0.17783;		%    |H|<=0.0.17783
As= -20*log10(deltas);
fs= 0.3;	% *Fs/2		(per 0.15*Fs<=|f|<=Fs/2)
% calcolo l'ordine del filtro di butterworth che rispetta le specifiche date
% e la freq. di taglio
% ATTENZIONE: le funzioni "buttord", "cheby1ord", "cheby2ord", "ellipord"
%           assumono le freq. di taglio normalizzate a Fs/2 (Nyquist) e non a Fs (freq. di campionamento)
% ATTENZIONE ancora: bisogna chiamare la funzione col parametro 's' per farla lavorare sul filtro Analogico
%      senza 's' opera direttamente sul filtro digitale e il risultato non è corretto
[n,fN0]= buttord(fp,fs,Attp,As,'s');
f0= fN0/2;	% così è normalizzata a Fs e non a Fs/2

% Calcolo la funzione di trasferimento (Laplace) del "prototipo analogico" Butterworth con w0= 1 [rad/s] (f0= 1/(2*pi) [Hz])
[zer,pol,k]= buttap(n);
% e visualizzo i poli, che sono SUL cerchio unitario
figure(1); clf;
zplane(zer,pol);
title('Poli del prototipo passabasso analogico normalizzato');

% trovo i coeff. della funzione di trasferimento, che occorrono per la trasformazione di frequenza (lp2lp)
bnorm= k*poly(zer);
anorm=   poly(pol);		% alternativa: [bnorm,anorm]= zp2tf(zer,pol,k);

% sostituisco s con s/(2*pi*f0) per avere un passabasso "denormalizzato"
[b,a]= lp2lp(bnorm,anorm,2*pi*f0);

% e visualizzo i poli DENTRO il cerchio unitario
figure(2); clf;
zplane(b,a);	% versione di "zplane" applicata alla FdT (x filtri numerici)
% difatti compare implicitamente lo zero di ord. 6 in z=0 poich? il rapporto di polinomi si assume in z^-1
% zplane(roots(b),roots(a)); avrebbe dato solo i poli (assumendo il rapporto di polinomi in "s")
title('Poli del filtro analogico');

% Applico il metodo di invarianza impulsiva per digitalizzare ilfiltro analogico
[bNum, aNum]= impinvar(b,a);	% FdT del filtro numerico
% "impinvar" segue i passi algoritmici previsti dal Metodo dell'Invarianza Impulsiva
% (espansione in frazioni parziali e sostituzione dei poli con esponenziale)
% vedi anche la pagina di Help in linea su "impinvar"

% Si pu? constatare che i poli del filtro numerico non sono pi? sul cerchio unitario...
figure(3); clf;
zplane(bNum,aNum);
% gli zeri hanno un modulo molto grande...
axis([-2 2 -2 2]);	% cos? ne visualizziamo solo alcuni
title('Poli del filtro numerico');

% visualizzazione modulo in scala log e fase in gradi
Npunti= 128;	% numero di punti su cui calcolare H(f)
Fs= 1;	% assumo una freq. di Nyquist normalizzata
figure(4); clf;
freqz(bNum,aNum,Npunti,Fs);
title('Risposte in ampiezza e fase del filtro numerico');


% se voglio controllare la correttezza delle specifiche, preparo un plot
% in scala lineare e poi zoomo sui valori critici

% Ricavo la funzione di trasferimento del filtro numerico
[H,fvec]= freqz(bNum,aNum,Npunti,Fs);
% e la visualizzo assieme al ritardo di gruppo
figure(5); clf;
subplot(2,1,1);
plot(fvec,abs(H),'b');
grid; xlabel('f/Fs'); ylabel('|H|');
subplot(2,1,2);
plot(fvec,unwrap(angle(H)),'r--');
grid; xlabel('f/Fs'); ylabel('arg[H] (rad.)');
title('Risposta in ampiezza in scala lineare: controllare i margini');
figure(6); clf;
grpdelay(bNum,aNum,Npunti,Fs);