Sezione 5.2: Filtri digitali Su Capitolo 5: Filtri analogici e numerici Sezione 5.4: Filtraggio polifase 

5.3 Filtri numerici

Affrontiamo[226]  [226] Senza pretesa di rigore e completezza, essendo questi argomenti trattati anche in diversi altri corsi. ora il tema di come passare dai concetti esposti al § 5.2 agli schemi di calcolo da adottare nel caso di una implementazione completamente numerica, associata all’elaborazione di dati campionati. Come descritto al § 5.2.2, i campioni yk del risultato di una convoluzione y(t) = x(t) * h(t) tra segnali limitati in banda x(t) e h(t) si ottengono mediante la convoluzione discreta (10.98)
yk = Tn = −∞hnxk − n
x = randn(1,1000);   # segnale di ingresso
Nx = 1000;           # campioni di ingresso
h = [1 0.5 0.25 0.125 0.0625]; # coeff. filtro
Nc = 5;              # numero coeff.
Ny = Nx+Nc-1;        # campioni di uscita
y = zeros(1,Ny);     # vettore di uscita
for k = 1:Ny         # campioni di uscita
 for n = 1:Nc        # convoluzione
  if ( (k-n+1)<=Nx && n<=k ) # controllo
   y(k) = y(k) + h(n)*x(k-n+1);
  endif
 end
end
equivalente a campionare l’uscita y(t) di un filtro trasversale con coefficienti cn = Thn proporzionali ai campioni hn della risposta impulsiva h(t) del filtro analogico di partenza. Anziché campionare l’uscita del filtro digitale, l’operazione di filtraggio viene implementata per via hw o sw in modo da eseguire direttamente la convoluzione discreta yk = Nn = 0 hnxk − n, come ad es. svolto dal codice Octave mostrato sopra. Notiamo esplicitamente che nel filtraggio numerico viene omesso il fattore T = 1fc, che è inserito autonomamente dal filtro di restituzione (eq. (10.68)) nel momento del passaggio da numerico ad analogico.

5.3.1 Sintesi FIR a partire dalla descrizione tempo continua

Acquisiamo innanzitutto due metodi per ricavare un insieme finito di coefficienti hn per il filtro numerico, a partire dalla descrizione analogica del comportamento desiderato in termini di h(t) o di H(f).
Finestratura della risposta impulsiva
In questo approccio i valori di hn si ottengono campionando una versione finestrata hw(t) = h(t)w(t) della h(t) desiderata, e quindi ritardando gli stessi al fine di ottenere un filtro causale, come discusso all’esempio di pag. 1. Come noto dal § 3.8.4, l’operazione di finestratura produce una risposta in frequenza per il filtro sintetizzato pari a Hw(f) = H(f) * W(f), e dunque va posta particolare attenzione nella scelta della funzione finestra w(t).
Oscillazione uniforme in frequenza
E’ una tecnica iterativa[227]  [227] Vedi https://en.wikipedia.org/wiki/Parks-McClellan_filter_design_algorithm volta a rendere minimo il massimo errore di approssimazione tra l’andamento desiderato per |H(f)|, espresso nella forma di uno schema di tolleranza (vedi pag. 1), e l’andamento ottenibile esprimendo |H(f)| come una combinazione di polinomi di Chebyshev[228]  [228] Vedi http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/00700_OptimalFIR.pdf . Il risultato è una (f) che presenta oscillazioni ridotte ma a tutte le frequenze[229]  [229] Il filtro risultante è detto per questo equiripple., mentre i coefficienti hn si ottengono mediante dft inversa (§ 4.5) della sequenza Hn ottenuta campionando (f).

5.3.2 Trasformata zeta e filtraggio

Al § 4.5.1 abbiamo già discusso come una sequenza di campioni xn =  x(t)|t = nT di un segnale limitato in banda W = 12T possa essere descritta dalla sua trasformata zeta
X(z) = Z{xn} = n = −∞xn z− n
e come calcolando quest’ultima per z = e jω si ottenga una rappresentazione spettrale
X(e jω) =  X(z)|z = e jω = n = −∞xn ejωn
di {xn} che corrisponde alla periodizzazione della X(f) = F {x(t)}, in cui l’intervallo di frequenze  − fc2 < f < fc2 si mappa in quello  − π < ω < π, ovvero f = ω2πT.
Allo stesso modo si possono definire le trasformate zeta della sequenza di uscita da un filtro Y(z) = Z{yn} e dei campioni della risposta impulsiva H(z) = Z{hn} del filtro stesso, e dato che anche ora vale la corrispondenza tra convoluzione nel tempo e prodotto delle trasformate[230]  [230] Vedi ad es. https://it.wikipedia.org/wiki/Trasformata_zeta, si può scrivere
(10.103) Y(z) = H(z) X(z)
in cui H(z) è l’analogo della H(s) (eq. (10.90)) tempo-continuo e per questo prende il nome di funzione di trasferimento tempo-discreto. Sostituendo z = e jω si ottiene la risposta in frequenza tempo-discreto (10.81) H(e jω) = n = −∞hn e− jωn periodica in ω con periodo 2π, relativa alla sequenza risposta impulsiva hn, e che corrisponde alla dtft (§ 4.4) H(f) = n = −∞hn e −j2πfnT della stessa sequenza hn qualora si ponga ω = 2πfT, vedi anche la fig. 4.30 a pag. 1.

5.3.2.1 Filtri a risposta impulsiva finita

Affrontiamo la descrizione di un filtro fir mediante la trasformata zeta rifacendoci dapprima al caso del filtro del primo ordine (§ 5.2.3) descritto da una h(t) = δ(t) + αδ(t − T), e consideriamo al suo ingresso una sequenza impulsiva δn =  1    n = 0 0 n ≠ 0 al posto dell’impulso di Dirac δ(t). Il filtro è ora caratterizzato nei termini della sequenza di risposta all’impulso hn = δn + αδn − 1 la cui trasformata zeta risulta[231]  [231] E’ sufficiente applicare la definizione H(z) = n = −∞hn z − n = n = −∞(δn + αδn − 1)  z − n = 1 + αz−1, dato che Z{δn} = 1 e che un ritardo di m indici ha trasformata
Z{xn − m} = n = −∞xn − m z − n = k = −∞xk z − k − m = zmk = −∞xk z − k = zmX(z)
(si è posto n − m = k). In particolare, un ritardo unitario corrisponde al prodotto per z−1 della sequenza trasformata, e dunque Z{δn − 1} = z−1.
(10.104) H(z) = 1 + α z−1
figure f7.41.png
Allo schema computazionale di pag. 1 si affianca quindi quello algoritmico mostrato a lato, in cui il blocco z−1 rappresenta il ritardo di un campione, mentre il prodotto per α è ora raffigurato come una etichetta posta sull’arco da cui passa la sequenza in transito, cosicché possiamo scrivere
Y(z) = X(z) + α z− 1 X(z) = (1 + αz− 1) X(z)
ossia la (10.104) applicata alla (10.103). Estendendo ora la trattazione ad un generico filtro fir di ordine N descritto da una hn = Nk = 0akδn − k otteniamo
(10.105)
HFIR(z) = Nk = 0akzk = a0Nk = 0ak a0 zk = a0 Nn = 1(1 − cnz−1)
in cui gli N valori cn rappresentano gli zeri[232]  [232] Si dicono zeri di un polinomio P(z) = Nk = 0βkzk di grado N le radici z = cn (n = 1, 2, ⋯, N) dell’equazione P(z) = 0. La (10.105) può essere riscritta come H(z) = zNNk = 0akzN − k = Nk = 0akzN − k zN e dunque si azzera per N(z) = Nk = 0akzN − k = 0, che è un polinomio a potenze positive. Una volta trovate le sue radici cn possiamo scrivere N(z) = a0Nn = 1(z − cn) o equivalentemente H(z) = N(z)zN  = a0Nn = 1(1 − cnz−1). (reali, od in coppie complesse coniugate) di H(z).
Esercizio: media mobile E’ il filtro numerico definito a pag. 1 come un rect discreto, ovvero con hn = 1 per n = 0, 1, ⋯, N − 1 e zero altrimenti.
||sin(ω72)sin(ω2) ||
figure f7.41-medmob.png

posizione degli zeri figure f7.41-medmob-zeri.png
Tenendo conto che   N− 1n = 0 an = 1 − aN 1 − a ([233]  [233] Vedi ad es. https://it.wikipedia.org/wiki/Serie_geometrica), si ottiene
HMA(z) = N− 1n = 0zn = 1 − zN 1 − z−1
che calcolata per z =  e jω fornisce 
HMA(ω) = 1 − e −jωN 1 −  e −jω  =  e −jωN2  e −jω2 e jωN2 − e −jωN2  e jω2 −  e −jω2  = e −jω(N − 1)2 sin(ωN2) sin(ω2)
di cui il primo fattore è un termine di fase lineare mentre dal secondo se ne ricava il modulo
|HMA(ω)| = |||sin(ωN2) sin(ω2) |||
rappresentato in figura per N = 7, e che si dimostra periodico, con sei zeri equidistribuiti nell’intervallo 0 < ω < 2π. Infatti dall’espressione di H(z) = 1 − zN 1 − z−1 si ottengono sette zeri sul cerchio unitario in posizione cn =  e j2πn(N − 1) di cui il primo in z = 1 si cancella con l’unico polo nella stessa posizione.

5.3.2.2 Risposta impulsiva infinita

Partendo anche questa volta dal caso del primo ordine (§ 5.2.4) notiamo che campionando l’uscita si osserva la sequenza yn = xn + αyn − 1 da cui, trasformando ambo i membri, si ottiene Y(z) = X(z) + αz− 1Y(z) ovvero Y(z)(1 − αz− 1) = X(z) e dunque
H(z) = Y(z)X(z) = 1 1 − αz− 1
con un polo[234]  [234] Ossia una radice del denominatore D(z) di H(z), che scritta come H(z) = zz − α vale D(z) = z − α, e si azzera in z = α. in z = α. Ad un generico filtro iir di ordine N corrisponde una
(10.106)
HIIR(z) = Mk = 0akz− k 1 − Nk = 1bkz− k  =  zNMk = 0akzM − k zM(zN − Nk = 1bkzN − k) = a0Mn = 1(1 − cnz− 1)Nn = 1(1 − dnz− 1)
in cui (ak, bk) sono i coefficienti di numeratore e denominatore[235]  [235] Il rapporto tra polinomi viene normalizzato in modo da far risultare b0 = 1., e (cn, dn) le rispettive radici. Mentre il numeratore di (10.106) esprime la componente fir, la presenza del denominatore che conferisce alla H(z) poli diversi da z = 0 determina la componente con risposta infinita. Scrivendo infatti (10.106) come rapporto tra polinomi H(z) = A(z)B(z) la relazione Y(z) = H(z)X(z) diviene Y(z)B(z) = X(z)A(z) che, antitrasformata[236]  [236] I passaggi iniziano con lo scrivere Y(z)(1 − Nk = 1bkzk) = X(z)Mk = 0akzk ovvero Y(z) − Nk = 1bkzkY(z) = Mk = 0akzkX(z); dato ora che Z− 1{zkX(z)} = xn − k, si arriva al risultato mostrato., dà origine all’equazione alle differenze finite yn − Nk = 1bkyn − k = Mk = 0akxn − k ovvero
(10.107)
yn = Nk = 1bkyn − k + Mk = 0akxn − k
Questa espressione permette di descrivere il funzionamento del filtro in base al diagramma di fig. 5.19. Il secondo termine della (10.107) individua la sequenza intermedia un = Mk = 0akxn − k in modo da poter scrivere yn = Nk = 1bkyn − k + un. La forma diretta di fig. 5.19-a) altro non è che lo schema di un filtro trasversale ruotato in verticale e che calcola un a partire da xn, seguito dal blocco che calcola yn a partire da se stesso e da un. Dato che i due blocchi esprimono relazioni lineari tempo invarianti sussiste per essi la proprietà commutativa, espressa dalla forma canonica di fig. 5.19-b), in cui gli elementi di ritardo sono stati messi in comune, a tutto vantaggio della memoria necessaria ad implementare il filtro numerico.
a) figure f7.42.png
b) figure f7.42b.png
Figure 5.19 Architettura di filtri iir: a) forma diretta; b) forma canonica
Stabilità
I poli del denominatore della (10.106) devono giacere tutti all’interno del cerchio unitario |z| < 1 pena l’instabilità del filtro, sebbene siano ammessi poli con modulo unitario qualora posti in corrispondenza di uno zero, in modo da cancellarne l’effetto.
Sensibilità alla quantizzazione
La precedente considerazione rende evidente il problema legato alla realizzazione del filtro mediante operazioni a precisione finita: l’effetto di quantizzazione subìto dai coefficienti porta a variazioni delle posizioni di poli e zeri che possono determinare effetti indesiderati. Anche per questo motivo esistono architetture alternative[237]  [237] Vedi ad es. http://raffaeleparisi.site.uniroma1.it/didattica/circuiti-a-tempo-discreto/dispense-di-circuiti-a-tempo-discreto, capitolo 8. Alcune architetture si dimostrano (a parità di precisione) migliori di altre nel definire la posizione di zeri e poli nel piano z. a quella canonica, legate a modi di diversi di scrivere la (10.106) come ad esempio nel prodotto di fattori, che dà luogo ad una architettura di celle in cascata.
Comportamento della fase e complessità
Mentre per i filtri fir esiste la possibilità di ottenere una fase lineare (vedi pag. 1) e dunque non distorcente, per i filtri iir questo non è possibile; d’altra parte i secondi permettono di ottenere una buona rapidità di variazione della risposta in frequenza pur mantenendo basso l’ordine, e dunque il carico computazionale.

5.3.3 Sintesi di un filtro IIR a partire da un filtro analogico

Il progetto di filtri analogici si basa su metodi consolidati ed efficienti (§ 5.1), che danno luogo ad una rappresentazione nella variabile s del tipo della (10.90)[238]  [238] Per ridurre le possibilità di confusione, adottiamo il pedice a per riferirci al mondo analogico.
Ha(s) = Mm = 0amsm Nn = 0bnsn
a cui corrisponde una risposta impulsiva ha(t), ed una equazione differenziale
(10.108)
Nn = 0 bn dny(t) dtn  = Mm = 0 am dmx(t) dtm
A partire da queste grandezze sono stati individuati metodi che consentono di definire un filtro numerico più o meno equivalente[239]  [239] Nel senso che alimentando il filtro numerico con i campioni di un segnale si ottiene circa lo stesso risultato che campionando l’uscita del filtro analogico di partenza. ad uno analogico, alcuni dei quali illustriamo appresso.

5.3.3.1 Invarianza della risposta impulsiva

Questo approccio ottiene i campioni hn per il filtro numerico come hn = ha(nT), ovvero campionando ha(t) = F −1{Ha(f)}. Qualora Ha(s) abbia N poli in s = dk, per essa sussiste le sviluppo in frazioni parziali[240]  [240] Vedi ad es. https://it.wikipedia.org/wiki/Decomposizione_in_fratti_semplici Ha(s) = Nk = 1 Ak s − dk , a cui corrisponde una ha(t) = Nk = 1 Ak edkt e dunque una
hn = h(nT) = Nk = 1 Ak edkTn
la cui trasformata zeta vale[241]  [241] Vedi esempio a pag. 1
H(z) = Nk = 1 Ak 1 − edkTz−1
ovvero con poli in z = edkT = esT|s = dk ([242]  [242] Infatti con alcuni passaggi
H(z) = Nk = 1 Ak 1 − edkTz−1  = Nk = 1zAk z − edkT
può essere riscritta come
H(z) = k Mn = 1(z − cn)Nn = 1(z − ednT)
in cui però gli zeri cn vanno calcolati.
). L’ultima osservazione comporta l’aver fatto uso della trasformazione[243]  [243] Seppur limitata alla sola corrispondenza per la posizione dei poli, in quanto gli zeri di H(s) non si mappano nel piano z allo stesso modo di come fanno i poli, vedi nota precedente. z = esT, che garantisce il mantenimento della stabilità[244]  [244] Infatti scrivendo s come s = σ + jΩ si ottiene z = esT = eσTe jΩT, e dunque i poli dk = σk + jΩk di Ha(s) che giacciono nel semipiano negativo del piano s, ovvero con σk < 0, vengono mappati all’interno del cerchio unitario nel piano z, in quanto ad essi corrispondono poli per H(z) in z = zk = eσkTe jΩkT, per i quali |zk| = eσkT < 1..
La trasformazione z = esT mette inoltre in corrispondenza periodica l’asse immaginario nel piano s (s = jΩ = j2πf) con la circonferenza del cerchio unitario nel piano z (z = e jω), nel senso che ad ogni  − π < ω < π corrisponde una pulsazione Ω = 2ω  1 2T  =  ω T ; pertanto la fascia del piano s compresa tra  − fc 2  < f <  fc 2 viene mappata all’interno del cerchio unitario del piano z (vedi tab. 5.1), ed anche le fasce superiori ed inferiori subiscono la medesima sorte. Questo aspetto è una manifestazione del fenomeno dell’aliasing che si manifesta per segnali campionati, ed avviene qualora Ha(f) non soddisfi le condizioni di stretta limitazione in banda, rendendo il metodo idoneo solo alla progettazione di filtri passa basso o passa banda.
Un modo alternativo di vedere il problema si basa sulla considerazione che al filtro numerico corrisponde una a(t) = n = 0hnδ(t − nT), a cui come noto (§ 4.1) a sua volta corrisponde una Ha(f) = 1T n = −∞Haf − n T ; dunque se Ha(f) non è limitata in banda tra ± 12T si verifica aliasing.

5.3.3.2 Corrispondenza di poli e zeri

In questo metodo la relazione z = esT si estende[245]  [245] Vedi https://en.wikipedia.org/wiki/Matched_Z-transform_method. Che ci sia qualcosa di sensato nello scrivere z = esT è motivato anche dal fatto che applicando la definizione di trasformata di Laplace ad ha(t) otteniamo
Ha(s) = L{ha(t)} = −∞ha(t)e − stdt = −∞[n = 0hnδ(t − nT)]e − stdt = 
 = n = 0hn −∞δ(t − nT)e − stdt = n = 0hne − snT
e, dato che H(z) = n = −∞hn z − n, ne discende che Ha(s) =  H(z)|z = esT.
anche alla trasformazione degli zeri di Ha(s) in quelli di H(z): in pratica, dopo essere arrivati alla forma fattorizzata
Ha(s) =  Mn = 1(s − cn)Nn = 1(s − dn)
ogni zero cn (o polo dn) si trasforma in ecnT (o ednT) dando luogo a
H(z) =  Mn = 1(z − ecnT)Nn = 1(z − ednT) = k zM zN Mn = 1(1 − ecnTz−1)Nn = 1(1 − ednTz−1)
Sono mantenute le stesse corrispondenze tra piano s e piano z del § 5.3.3.1, così come la possibilità di aliasing. D’altra parte, i valori della risposta impulsiva hn non corrispondono più ai campioni di h(t).

5.3.3.3 Equazioni alle differenze

Qui si parte dalle equazioni differenziali (10.108) che sono a monte della H(s), e che vengono approssimate come equazioni alle differenze finite (10.107). Ciò porta alla equivalenza
s = 1 − z−1 T
che permette di ottenere la H(z) a partire dalla H(s) mediante un cambio di variabile, e che dà luogo alla corrispondenza mostrata al centro di tab. 5.1, in cui i punti dell’intero semipiano negativo del dominio s si mappano in punti interni alla circonferenza nel piano z di raggio 0.5 e centrata in z = 0.5. Pertanto ora sono evitati i problemi legati all’aliasing, ma la presenza dei poli solo nel semipiano destro del piano z impedisce l’uso del metodo per la progettazione di filtri passa-alto.
metodo invarianza h(t),
poli e zeri
eq. alle differenze trasf. bilineare
relazione z = esT s =  1 − z−1 T s =  2 T 1 − z−1 1 + z−1
mapping figure f7.43.png figure f7.43b.png figure f7.43c.png
debolezze aliasing si può agire solo a basse frequenze distorsione dell’asse della frequenza
Table 5.1 Aspetti peculiari delle trasformazioni H(s) ⇒ H(z)

5.3.3.4 Trasformazione bilineare

Anche quest’ultimo metodo, detto di Tustin, trae origine dalla approssimazione numerica di una equazione differenziale, ma può essere anche visto come una approssimazione di primo ordine[246]  [246] https://en.wikipedia.org/wiki/Bilinear_transform della relazione z = esT analizzata in precedenza. Il metodo di Tustin si basa sul cambio di variabile
s = 2T 1 − z−1 1 + z−1
che corrisponde a mappare l’intero semipiano destro del piano s all’interno della circonferenza di raggio unitario del piano z, come mostrato alla destra di tab. 5.1.
Anche in questo caso il risultato è stabile e privo di aliasing; d’altra parte si verifica invece una distorsione dell’asse delle frequenze, dato che la fase ω di z = e jω è ora legata alla frequenza f del filtro analogico tramite la relazione ω = 2 arctanπf fc per cui va adottata una fc = 1T abbastanza più elevata delle frequenze di interesse del filtro[247]  [247] In tal caso infatti si può ritenere l’arcotangente approssimativamente lineare., oppure si progetta un filtro analogico che tenga conto in partenza della distorsione a cui verrà sottoposta la sua risposta in frequenza.
 Sezione 5.2: Filtri digitali Su Capitolo 5: Filtri analogici e numerici Sezione 5.4: Filtraggio polifase