%------------------------------------------------------------------------
%               DAI FILTRI ANALOGICO AI FILTRI DIGITALI
% Esempio per progettare l'analogo discreto di un filtro RC
% Metodi:  - Inv. Risp. Impulsiva (IIR e FIR)
%          - IIR con Trasformata Bilineare
%------------------------------------------------------------------------

clear all  % cancella tutte le variabili dallo spazio di memoria
close all  % chiude tutte le figure aperte

f_cut         = 1e6;          % Frequenza di taglio
Fc            = 5*f_cut ;     % Frequenza di campionamento
DurataFIRnorm = 4        ;    % Durata Finetra FIR norm. (multiplo di RC)
Nf            = 128;          % Numero punti in [0,Fc] su cui si valuta lo spettro

% ----- calcolo parametri ausiliari -------------------------------------
Tc    = 1/Fc ;                % Intervallo di campionamento
RC    = 1/(2*pi*f_cut);       % Costante di tempo
alfa  = 1/RC;                 % p_s = -alfa  polo Tr. Laplace

f     = [0: Fc/Nf : 2*Fc];            % Asse delle frequenze (continue)
H_f   = alfa./(alfa+sqrt(-1)*2*pi*f); % H(f) analogica filtro RC

t     = [0:Tc/10:10*RC];       % Asse dei tempi continuo fino a 10*(RC)
h_t   = alfa*exp(-alfa.*t);    % Risposta impulsiva tempo-continua


% --- FILTRO FIR per finestramento rettangolare di h[n] IIR --------------
DeltaT     = DurataFIRnorm*RC ;    % Durata finestra temporale
N          = ceil(DeltaT/Tc)+1;    % Campioni risp. impulsiva discreta finestrata.

n          = [0:N-1];              % Supporto temporale discreto
h_n        = alfa*exp(-alfa*n*Tc); % Campionamento risposta impulsiva
h_n(1)     = h_n(1)/2;             % Compensazione discontinuita in t =0; 
H_k        = (1/N)*fft(h_n,Nf);    % DFT (Matlab senza 1/N) fatta su Nf punti
F_k        = [0:Nf-1]*Fc/Nf;       % Frequenze discrete che corrispondono alla DFT
                                   % dove si campiona la DTFT della sequenza


% --- FILTRO IIR a un polo INV. RISP. IMP. (con zero comp. discontinuità in t=0)
a     = exp(-alfa*Tc); % polo Tr. ZETA in funzione del polo (-alfa) Tr. Laplace

% h(n) = alfa*( a^n u[n] - 1/2 delta[n] )
% H(z) = alfa * 0.5( 1 + a z^-1)/(1 - a z^-1) = B(z)/A(z)

B = 0.5* alfa *[1,a]; % Coefficienti polinomio B(z)
A = [1,-a];           % Coefficienti polinomio A(z)

% Risposta in frequenza per camp. sul cerchio unitario della trasformata ZETA
[H,F] = freqz(B,A,Nf,'whole',Fc); 

% --- FILTRO IIR a un polo con TRASF. BILINEARE (senza compensazione
% f_cut)  --

pz = (1-alfa*Tc/2)/(1+alfa*Tc/2) ; % Polo in tr.ZETA da polo (-alfa) Tr. Laplace

B = alfa/(alfa+2/Tc)*[1,1 ];       %  Coefficienti polinomio B(z)
A = [1,-pz];                       %  Coefficienti polinomio A(z)

% Risposta in frequenza per camp. sul cerchio unitario della trasformata ZETA
[Hb,F] = freqz(B,A,Nf,'whole',Fc);       
                                  
                                  
% --- FILTRO IIR a un polo con TRASF. BILINEARE (con compensazione f_cut) ----

f_cut_comp = tan(pi*f_cut/Fc)/(pi*Tc); % Compensazione distorsione in frequenza
RC2        = 1/(2*pi*f_cut_comp);
alfa       = 1/RC2 ;

pz = (1-alfa*Tc/2)/(1+alfa*Tc/2) ;     % Polo in tr.ZETA da polo (-alfa) Tr. Laplace

B  = alfa/(alfa+2/Tc)*[1,1];           %  Coefficienti polinomio B(z)
A  = [1,-pz];                          %  Coefficienti polinomio A(z)

% Risposta in frequenza per camp. sul cerchio unitario della trasformata ZETA
[Hbc,F] = freqz(B,A,Nf,'whole',Fc);    

                                  
% ----  Grafici --------------------------------------------------------------                                
                                 
figure;                                % crea nuova figura vuota
plot( t , h_t ,'b'); hold on; grid on; % Grafico risposta impulsiva tempo-continua
stem( n*Tc , h_n, '-or' );             % Grafico risposta impulsiva discretizzata

xlabel('time (sec)'); ylabel('RC impulse response') % legenda assi
RC_str = num2str(RC,2) ;      % costante tempo in formato carattere
TC_str = num2str(Tc,2);       % Intervallo campion. in formato stringa caratteri

% stringa del titolo
string_title = ['RC Filter imp. response (RC = ', RC_str,' sec ; TC = ',TC_str,' sec)']; 
title( string_title );                        % creazione titolo figura
legend( 'h(t)' , 'h[nT_c] Inv. Risp. Imp.');  % legenda grafico


figure;                               % Crea nuova figura vuota
plot( f , 20*log10(abs(H_f)));        % |H(f)|^2 analogica (senza aliasing) in dB 
hold on; grid on ;                    % mantiene grafico e inserisce  griglia
plot( F , 20*log10(Tc*abs(H)),'-.xk');% |H[k]|^2  filtro IIR Inv.R.Imp. in dB
plot( F_k , 20*log10((N*Tc)*abs(H_k)),'-or');  % |H[k]|^2 finestramento FIR  in dB
plot( F , 20*log10(abs(Hb)),':sg')    % |H[k]|^2  filtro IIR T.Bil. non comp. in dB
plot( F , 20*log10(abs(Hbc)),'-^m')   % |H[k]|^2  filtro IIR T.Bil compensata in dB

axis([0,Fc/2,-30,1]); % Zoom del grafico [fmin, fmax, Hmin_dB , Hmax_dB] 



xlabel('frequency (Hz)'); ylabel('|H(f)|^2  (dB)')  % legenda assi
fcut_str = num2str(f_cut,2);    % Freq. di taglio in formato stringa di caratteri
Fc_str   = num2str(Fc,2);       % Freq. campionamento in formato stringa caratteri
% stringa del titolo
string_title = ['RC Filter Frequency Response ( f_{cut} = ',fcut_str,' Hz ; F_c = ',Fc_str,' Hz)']; 
title( string_title );          % creazione titolo figura
% legenda grafico
legend( '| H(f) | analogica','| H[k] | IIR InvRisImp','| H[k] | FIR','| H[k] | IIR TrBil','| H[k] | IIR TrBil comp'); 


figure;                        % crea nuova figura vuota
plot( f , abs(H_f));           % Grafico di |H(f)|^2 analogica (senza aliasing)
hold on; grid on;              % mantiene grafico e inserisce griglia

stem( F , Tc*abs(H),'-.xk');       % |H[k]| IIR con invrianza risposta impulsiva
stem( F_k ,(N*Tc)*abs(H_k),'-or'); % |H[k]| FIR normalizz. alla durata N*Tc di h[n]
stem( F , abs(Hb),':sg');          % |H[k]| IIR con Trasf. Bilineare non compensata
stem( F , abs(Hbc),'-^m');         % |H[k]| IIR con Trasf. Bilineare compensata

axis([0,Fc/2,0,1.2]); % Zoom del grafico Zoom del grafico [fmin, fmax, Hmin , Hmax]

xlabel('frequency (Hz)'); ylabel('|H(f)| ')  % Legenda assi
fcut_str = num2str(f_cut,2) ;      % Freq. di taglio in formato stringa caratteri
title( string_title );             % Creazione titolo figura
% legenda grafico
legend( '| H(f) | analogica','| H[k] | IIR InvRisImp','| H[k] | FIR','| H[k] | IIR TrBil','| H[k] | IIR TrBil comp');

                                     
                                     
                                     
                                     
                                     






