Trasformazione bilineare

Un filtro IIR si può progettare partendo dalla grande raccolta di soluzioni analogiche già note (Butterworth, Chebyshev, ellittici: Filtri di ButterworthIl filtro di Butterworth analogico di ordine $N$ e pulsazione di taglio a $-3$ dB $\Omega_0$ ha $|H_a(j\Omega)|^2=\frac1{1+(\Omega/\Omega_0)^{2N}}$: massimamente piatto in $\Omega=0$, monotono decrescente, valore $\frac1{\sqrt2}$ in $\Omega_0$. È un filtro tutti-poli: i $2N$ poli di $H_a(s)H_a(-s)$ stanno su una circonferenza di raggio $\Omega_0$, $H_a(s)$ prende i $N$ poli del semipiano sinistro $s_k=\Omega_0e^{j\pi(\frac12+\frac{2k-1}{2N})}$, $H_a(s)=\Omega_0^N/\prod(s-s_k)$. Con la bilineare $H(z)=A\frac{(1+z^{-1})^N}{\prod(1-p_kz^{-1})}$, $p_k=\frac{1+s_k}{1-s_k}$, zero di molteplicità $N$ in $z=-1$, $A$ tale che $H(1)=1$; il modulo digitale è $|H(e^{j\hat\omega})|^2=\frac1{1+(\tan(\hat\omega/2)/\Omega_0)^{2N}}$. Il filtro è fissato da $N$ e $\Omega_0$. Dalle specifiche ($\varepsilon^2=10^{R_p/10}-1$, $A=10^{R_s/20}$, $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$): $N=\Big\lceil\frac{\log_{10}\frac{A^2-1}{\varepsilon^2}}{2\log_{10}(\Omega_s/\Omega_p)}\Big\rceil$ e $\Omega_0$ fra $\Omega_p\varepsilon^{-1/N}$ e $\Omega_s(A^2-1)^{-1/2N}$.Filtri di Butterworth →). Serve una regola per trasformare un sistema analogico, descritto da una funzione di trasferimento $H_a(s)$ (Laplace unilatera, funzione di trasferimento e stabilitàCon la trasformata unilatera $X(s)=\int_{0^-}^\infty x(t)e^{-st}dt$ le derivate diventano $\mathcal L[y']=sY-y(0^-)$: un problema di Cauchy causale si riduce a un'equazione algebrica, $Y=\frac{b}{a}X+\frac{\text{termini iniziali}}{a}$ (forzata + libera). La funzione di trasferimento è $H(s)=b(s)/a(s)$; il sistema causale è BIBO stabile se e solo se $H$ è propria e tutti i poli hanno parte reale negativa.Laplace unilatera, funzione di trasferimento e stabilità →), in uno digitale, descritto da $H(z)$ (Trasformata zeta e sistemi a tempo discretoUn sistema LTI discreto è descritto da un'equazione alle differenze Σ a_i y(n-i) = Σ b_i x(n-i) con condizioni iniziali. L'uscita è evoluzione libera (dalle condizioni iniziali) più risposta forzata g*x. Si calcola in tre modi: soluzione dell'equazione, risposta impulsiva con segnale fittizio, trasformata zeta. FIR: memoria finita, sempre stabile; IIR: stabile se tutti i poli hanno modulo minore di 1. In frequenza G(f) = Σ b_i e^{-i2πfi} / Σ a_i e^{-i2πfi}.Trasformata zeta e sistemi a tempo discreto →). La regola più usata è la trasformazione bilineare. In questa nota $\Omega$ è la pulsazione analogica (rad/s) e $\hat\omega$ quella digitale normalizzata.

Il metodo di trasformazione $s\to z$

Date le specifiche del filtro digitale (modulo desiderato $D(e^{j\hat\omega})$):

1. con un'opportuna mappa $z\to s$, la risposta desiderata digitale diventa una risposta desiderata analogica $D_a(j\Omega)$;
2. si risolve il problema analogico con le tecniche note, $H_a(j\Omega)\approx D_a(j\Omega)$;
3. con la mappa inversa $s\to z$ si ottiene $H(e^{j\hat\omega})\approx D(e^{j\hat\omega})$, che eredita le proprietà della soluzione analogica.

La mappa $s\to z$ deve avere tre proprietà:

1. funzioni razionali in $s$ vanno in funzioni razionali in $z$ (altrimenti non si ottiene un'equazione alle differenze finita);
2. un sistema stabile in $s$ (poli con $\operatorname{Re}s<0$, Causalità, stabilità e connessioni dei sistemi LTIPer un sistema LTI le proprietà si leggono sulla risposta impulsiva $h$: causale se $h$ è nulla per $t<0$; reale se $h$ è reale; statico se $h=K\delta$; BIBO stabile se e solo se $h$ è assolutamente integrabile (sommabile nel discreto). Serie e parallelo corrispondono a $h_1*h_2$ e $h_1+h_2$.Causalità, stabilità e connessioni dei sistemi LTI →) va in un sistema stabile in $z$ (poli con $|z|<1$, Funzione di sistema, poli, zeri e stabilitàLa funzione di sistema H(z) è la trasformata zeta della risposta impulsiva: con ingresso z^n l'uscita è H(z) z^n. Per un FIR H(z) = Σ b_k z^{-k} è un polinomio con M zeri e M poli in z = 0; in generale H = B(z)/A(z) dall'equazione alle differenze. Sulla circonferenza unitaria H(e^{jω̂}) è la risposta in frequenza: |H| = prodotto delle distanze dagli zeri / prodotto delle distanze dai poli, quindi gli zeri bloccano frequenze e i poli le esaltano. Un LTI causale è BIBO stabile se e solo se tutti i poli hanno modulo < 1 (a meno di cancellazioni polo-zero); i FIR sono sempre stabili.Funzione di sistema, poli, zeri e stabilità →);
3. l'asse immaginario $s=j\Omega$ va nella circonferenza unitaria $z=e^{j\hat\omega}$ (così la risposta in frequenza analogica diventa quella digitale).

Una trasformazione che le soddisfa è la bilineare.

La trasformazione bilineare

Definizione. La trasformazione bilineare (o di Möbius) è $s=\dfrac{\alpha z+\beta}{\gamma z+\delta}$, con $\alpha,\beta,\gamma,\delta\in\mathbb C$ e $\alpha\delta-\beta\gamma\ne0$. È una mappa conforme (conserva gli angoli) che porta circonferenze in circonferenze (eventualmente degeneri: rette). Imponendo che il cerchio unitario vada nell'asse immaginario si ottiene $$s=\frac{z-1}{z+1},\qquad z=\frac{1+s}{1-s}.$$

I tre parametri essenziali si fissano con tre corrispondenze:

Esempio. $s=-0{,}5\Rightarrow z=\frac{1-0{,}5}{1+0{,}5}=\frac13$: un polo analogico reale in $-0{,}5$ diventa un polo digitale in $\frac13$ (stabile).

Conservazione della stabilità. Con $s=\sigma+j\Omega$ è $z=\frac{1+\sigma+j\Omega}{1-\sigma-j\Omega}$ e $$|z|^2=\frac{(1+\sigma)^2+\Omega^2}{(1-\sigma)^2+\Omega^2}=\frac{1+2\sigma+\sigma^2+\Omega^2}{1-2\sigma+\sigma^2+\Omega^2}.$$ Il numeratore è minore del denominatore se e solo se $\sigma<0$. Quindi $\operatorname{Re}s<0\iff|z|<1$ (e $\sigma=0\iff|z|=1$). Il semipiano sinistro va nel disco unitario, l'asse immaginario sulla circonferenza, il semipiano destro fuori. Esempio: $s=-0{,}3+0{,}5j$ dà $|z|^2=\frac{0{,}7^2+0{,}5^2}{1{,}3^2+0{,}5^2}=\frac{0{,}74}{1{,}94}=0{,}3814<1$.

Cosa fa al sistema: $H(z)=H_a\big(\frac{z-1}{z+1}\big)$

Dato $H_a(s)$ si definisce $H(z)=H_a(s)\big|_{s=(z-1)/(z+1)}$. La bilineare cambia il dominio della funzione, non i suoi valori: $H(z)$ in $z$ assume gli stessi valori che $H_a$ assume nel punto $s=\frac{z-1}{z+1}$. Di conseguenza se $H_a$ ha uno zero o un polo in $\hat s$, $H$ ha uno zero o un polo in $$\hat z=\frac{1+\hat s}{1-\hat s}.$$

Esempio 1. $H_a(s)=\frac1{s+1}$ (polo in $-1$, zero all'infinito): $H(z)=\frac{1}{\frac{z-1}{z+1}+1}=\frac{z+1}{2z}=\frac{1+z^{-1}}{2}$, con polo in $\hat z=\frac{1-1}{1+1}=0$ e zero in $z=-1$ (immagine dello zero all'infinito).

Esempio 2 (integratore). $H_a(s)=\frac1s$ dà $H(z)=\frac{z+1}{z-1}=\frac{1+z^{-1}}{1-z^{-1}}$, con polo in $z=1$ (immagine di $s=0$) e zero in $z=-1$; l'equazione alle differenze è $y[n]=y[n-1]+x[n]+x[n-1]$, un accumulatore (integrazione col metodo dei trapezi).

Forma generale. Se $H_a(s)=b_0\frac{\prod_{k=1}^M(s-z_{a,k})}{\prod_{k=1}^N(s-p_{a,k})}$ con $M\le N$, sostituendo $s=\frac{z-1}{z+1}$ e usando $s-\hat s=\frac{(1-\hat s)z-(1+\hat s)}{z+1}=\frac{(1-\hat s)\,(z-\hat z)}{z+1}$ si ottiene $$H(z)=A\,\frac{(1+z^{-1})^{N-M}\prod_{k=1}^M(1-z_kz^{-1})}{\prod_{k=1}^N(1-p_kz^{-1})},\quad z_k=\frac{1+z_{a,k}}{1-z_{a,k}},\ p_k=\frac{1+p_{a,k}}{1-p_{a,k}},\ A=b_0\frac{\prod(1-z_{a,k})}{\prod(1-p_{a,k})}.$$ Gli $N-M$ zeri all'infinito di $H_a$ diventano uno zero di molteplicità $N-M$ in $z=-1$. Per un all-pole di ordine $N$ ($M=0$): $H(z)=A(1+z^{-1})^N/\prod(1-p_kz^{-1})$.

Relazione fra le frequenze: il warping

Dove si trasformano le frequenze? Sulla circonferenza $z=e^{j\hat\omega}$: $$s=\frac{e^{j\hat\omega}-1}{e^{j\hat\omega}+1}=\frac{e^{j\hat\omega/2}\big(e^{j\hat\omega/2}-e^{-j\hat\omega/2}\big)}{e^{j\hat\omega/2}\big(e^{j\hat\omega/2}+e^{-j\hat\omega/2}\big)}=\frac{2j\sin(\hat\omega/2)}{2\cos(\hat\omega/2)}=j\tan\frac{\hat\omega}2 .$$ Dunque $s=j\Omega$ con

$$\boxed{\Omega=\tan\frac{\hat\omega}{2}\qquad\Longleftrightarrow\qquad\hat\omega=2\arctan\Omega},$$

una relazione non lineare: $H(e^{j\hat\omega})=H_a(j\Omega)\big|_{\Omega=\tan(\hat\omega/2)}$. È inevitabile: l'asse infinito $\Omega\in(-\infty,+\infty)$ deve entrare nell'intervallo finito $\hat\omega\in(-\pi,\pi)$. Si ha una distorsione dell'asse delle frequenze (frequency warping): vicino a $0$ è quasi lineare ($\tan\frac{\hat\omega}2\simeq\frac{\hat\omega}2$), poi l'asse viene sempre più compresso e tutte le frequenze analogiche molto grandi finiscono vicino a $\pi$.

Esempio. $\hat\omega/\pi=0{,}1,\ 0{,}25,\ 0{,}5,\ 0{,}75,\ 0{,}9$ corrispondono a $\Omega=0{,}158,\ 0{,}414,\ 1{,}000,\ 2{,}414,\ 6{,}314$. La metà dell'asse digitale, $\hat\omega=\pi/2$, è $\Omega=1$.

Per questo la bilineare serve per filtri selettivi in frequenza (a tratti costanti: passa-basso, passa-alto, passa-banda, elimina-banda), dove conta solo dove sono le frequenze di taglio e la distorsione si compensa; non serve per approssimare una risposta con un andamento prescritto continuo (per esempio un differenziatore).

Effetto sulle specifiche di un passa-basso

Le tolleranze in banda passante e in banda oscura sono conservate (stessi valori di $\delta_p,\delta_s$), si spostano solo le frequenze:

• banda passante: $0\le\hat\omega\le\hat\omega_p\iff0\le\Omega\le\Omega_p=\tan\frac{\hat\omega_p}2$, con $1-\delta_p\le|H_a|\le1$ che diventa $1-\delta_p\le|H|\le1$;
• banda oscura: $\hat\omega\ge\hat\omega_s\iff\Omega\ge\Omega_s=\tan\frac{\hat\omega_s}2$, con $|H_a|\le\delta_s$ che diventa $|H|\le\delta_s$.

Il numero di oscillazioni del modulo si conserva: se $|H_a|$ oscilla un certo numero di volte, anche $|H|$ oscilla altrettanto. Il procedimento è quindi: trasforma le frequenze delle specifiche con $\Omega=\tan\frac{\hat\omega}2$ (prewarping), progetta il filtro analogico con queste frequenze, poi applica la bilineare.

Esempio (prewarping). Passa-basso del primo ordine con taglio a $-3$ dB a $1$ kHz e $F_s=8$ kHz. È $\hat\omega_c=2\pi\frac{1000}{8000}=\frac\pi4$ e $\Omega_c=\tan\frac\pi8=0{,}4142$. Il prototipo analogico con taglio $\Omega_c$ è $H_a(s)=\frac{\Omega_c}{s+\Omega_c}$ e $$H(z)=\frac{\Omega_c}{\frac{z-1}{z+1}+\Omega_c}=\frac{\Omega_c}{1+\Omega_c}\cdot\frac{1+z^{-1}}{1-\frac{1-\Omega_c}{1+\Omega_c}z^{-1}}=\frac{0{,}2929\,(1+z^{-1})}{1-0{,}4142\,z^{-1}},$$ con $|H(e^{j\pi/4})|=0{,}7071$: il taglio cade proprio a $1$ kHz. Se invece si fosse usato $\Omega_c=\hat\omega_c=0{,}7854$ senza prewarping, il $-3$ dB cadrebbe a $\hat\omega=0{,}424\pi$, cioè $1695$ Hz: un errore del 70%.

Esempio ($\Omega_0=1$). Con $H_a(s)=\frac1{s+1}$ è $H(z)=\frac{1+z^{-1}}2$ e $|H(e^{j\hat\omega})|=\cos\frac{\hat\omega}2$: vale $0{,}7071$ in $\hat\omega=\pi/2$ (è proprio $\Omega=1$) e $0$ in $\pi$. Stessa forma di $\frac1{\sqrt{1+\Omega^2}}$, ma con l'asse compresso (nel grafico, con l'asse delle frequenze portato alla stessa unità, la curva analogica sarebbe quella tratteggiata).

Effetto sulla fase

Se $H_a(j\Omega)$ ha fase quasi lineare in banda passante, $\angle H_a(j\Omega)\simeq-\tau\Omega$, la fase digitale è $\angle H(e^{j\hat\omega})=\angle H_a\big(j\tan\frac{\hat\omega}2\big)\simeq-\tau\tan\frac{\hat\omega}2$, che non è lineare in $\hat\omega$ se non per $\hat\omega$ piccolo ($\tan\frac{\hat\omega}2\simeq\frac{\hat\omega}2$). Quindi la bilineare conserva il modulo, ma non la fase: un IIR così ottenuto non ha fase lineare (Filtri IIR - definizione e confronto con i FIRUn filtro IIR (infinite impulse response) è un sistema LTI descritto da $y[n]=\sum_{\ell=1}^{N}a_\ell y[n-\ell]+\sum_{k=0}^{M}b_kx[n-k]$: l'uscita usa anche le uscite passate (retroazione), per questo si chiama ricorsivo. Con le condizioni di riposo iniziale è LTI, $H(z)=\frac{\sum b_kz^{-k}}{1-\sum a_\ell z^{-\ell}}$ è un rapporto di polinomi, l'ordine è $N$ (numero di poli) e la risposta impulsiva ha durata infinita. Nel primo ordine $y[n]=a_1y[n-1]+b_0x[n]$ si ha $h[n]=b_0a_1^nu[n]$, ROC $|z|>|a_1|$, stabile se $|a_1|<1$; il gradino dà $b_0\frac{1-a_1^{n+1}}{1-a_1}\to\frac{b_0}{1-a_1}$. Si implementa iterando l'equazione alle differenze, non con la convoluzione. Rispetto ai FIR gli IIR rispettano le stesse specifiche di modulo con ordine molto più basso, ma possono essere instabili e non hanno fase lineare. Progetto: per tentativi (notch), con la trasformazione $s\to z$ dai filtri analogici, o con ottimizzazione numerica.Filtri IIR - definizione e confronto con i FIR →, Sistemi a fase lineare e assenza di distorsioneUn sistema non deforma il segnale se y[n] = K x[n-n0]: modulo costante e fase lineare -n0 w, cioè ritardo di gruppo costante n0. Un sistema reale e causale ha fase lineare (generalizzata) se e solo se è FIR con risposta impulsiva simmetrica h[n] = h[N-n] (ampiezza pari, fase -N w/2) o antisimmetrica h[n] = -h[N-n] (ampiezza dispari, fase -N w/2 + pi/2). Il ritardo di gruppo è N/2: intero se N è pari (vale la condizione di non distorsione), semi-intero se N è dispari (uscita interpolata e ritardata). Un IIR causale non può essere simmetrico, quindi non ha fase lineare esatta.Sistemi a fase lineare e assenza di distorsione →).

Un esempio completo di secondo ordine

Butterworth analogico di ordine $2$ con $\Omega_0=1$: $H_a(s)=\frac{1}{s^2+\sqrt2\,s+1}$, con poli in $s=\frac{-1\pm j}{\sqrt2}$ (Filtri di ButterworthIl filtro di Butterworth analogico di ordine $N$ e pulsazione di taglio a $-3$ dB $\Omega_0$ ha $|H_a(j\Omega)|^2=\frac1{1+(\Omega/\Omega_0)^{2N}}$: massimamente piatto in $\Omega=0$, monotono decrescente, valore $\frac1{\sqrt2}$ in $\Omega_0$. È un filtro tutti-poli: i $2N$ poli di $H_a(s)H_a(-s)$ stanno su una circonferenza di raggio $\Omega_0$, $H_a(s)$ prende i $N$ poli del semipiano sinistro $s_k=\Omega_0e^{j\pi(\frac12+\frac{2k-1}{2N})}$, $H_a(s)=\Omega_0^N/\prod(s-s_k)$. Con la bilineare $H(z)=A\frac{(1+z^{-1})^N}{\prod(1-p_kz^{-1})}$, $p_k=\frac{1+s_k}{1-s_k}$, zero di molteplicità $N$ in $z=-1$, $A$ tale che $H(1)=1$; il modulo digitale è $|H(e^{j\hat\omega})|^2=\frac1{1+(\tan(\hat\omega/2)/\Omega_0)^{2N}}$. Il filtro è fissato da $N$ e $\Omega_0$. Dalle specifiche ($\varepsilon^2=10^{R_p/10}-1$, $A=10^{R_s/20}$, $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$): $N=\Big\lceil\frac{\log_{10}\frac{A^2-1}{\varepsilon^2}}{2\log_{10}(\Omega_s/\Omega_p)}\Big\rceil$ e $\Omega_0$ fra $\Omega_p\varepsilon^{-1/N}$ e $\Omega_s(A^2-1)^{-1/2N}$.Filtri di Butterworth →). Sostituendo $s=\frac{z-1}{z+1}$ e moltiplicando numeratore e denominatore per $(z+1)^2$: $$H(z)=\frac{(z+1)^2}{(z-1)^2+\sqrt2\,(z-1)(z+1)+(z+1)^2}=\frac{(z+1)^2}{(2+\sqrt2)\,z^2+(2-\sqrt2)}=\frac{1}{2+\sqrt2}\cdot\frac{1+2z^{-1}+z^{-2}}{1+\frac{2-\sqrt2}{2+\sqrt2}\,z^{-2}}=0{,}29289\,\frac{1+2z^{-1}+z^{-2}}{1+0{,}17157\,z^{-2}}.$$ Nel denominatore il coefficiente di $z$ è $-2+2=0$. Si verifica che i poli sono $\hat z=\frac{1+s_k}{1-s_k}=\pm0{,}4142\,j$ (dentro il disco, come deve essere), che c'è uno zero doppio in $z=-1$ (immagine dei due zeri all'infinito di $H_a$), che $H(1)=0{,}29289\cdot\frac{4}{1{,}17157}=1$, $|H(e^{j\pi/2})|=0{,}7071$ (il taglio $\Omega_0=1$ cade in $\hat\omega=\pi/2$) e $H(e^{j\pi})=0$ (controllato con scipy.signal.bilinear).

Domande d'esame

1. Quali proprietà deve avere la mappa $s\to z$ per il progetto di filtri IIR e come le soddisfa la trasformazione bilineare? Traccia: razionali in razionali, stabilità conservata ($|z|^2$ in funzione di $\sigma$), asse immaginario nel cerchio unitario; punti notevoli $s=0,\infty,-1$.
2. Ricava la relazione fra la frequenza analogica e quella digitale nella bilineare e discuti le conseguenze. Traccia: $z=e^{j\hat\omega}$, $s=j\tan\frac{\hat\omega}2$, warping; compressione dell'asse; modulo conservato (tolleranze e oscillazioni), fase non lineare; prewarping delle frequenze di taglio; utile solo per filtri selettivi in frequenza.
3. Spiega come si progetta un filtro IIR passa-basso con la bilineare a partire dalle specifiche digitali. Traccia: $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$, progetto analogico, sostituzione $s=\frac{z-1}{z+1}$, poli in $\frac{1+s_k}{1-s_k}$, zeri in $z=-1$ (Filtri di ButterworthIl filtro di Butterworth analogico di ordine $N$ e pulsazione di taglio a $-3$ dB $\Omega_0$ ha $|H_a(j\Omega)|^2=\frac1{1+(\Omega/\Omega_0)^{2N}}$: massimamente piatto in $\Omega=0$, monotono decrescente, valore $\frac1{\sqrt2}$ in $\Omega_0$. È un filtro tutti-poli: i $2N$ poli di $H_a(s)H_a(-s)$ stanno su una circonferenza di raggio $\Omega_0$, $H_a(s)$ prende i $N$ poli del semipiano sinistro $s_k=\Omega_0e^{j\pi(\frac12+\frac{2k-1}{2N})}$, $H_a(s)=\Omega_0^N/\prod(s-s_k)$. Con la bilineare $H(z)=A\frac{(1+z^{-1})^N}{\prod(1-p_kz^{-1})}$, $p_k=\frac{1+s_k}{1-s_k}$, zero di molteplicità $N$ in $z=-1$, $A$ tale che $H(1)=1$; il modulo digitale è $|H(e^{j\hat\omega})|^2=\frac1{1+(\tan(\hat\omega/2)/\Omega_0)^{2N}}$. Il filtro è fissato da $N$ e $\Omega_0$. Dalle specifiche ($\varepsilon^2=10^{R_p/10}-1$, $A=10^{R_s/20}$, $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$): $N=\Big\lceil\frac{\log_{10}\frac{A^2-1}{\varepsilon^2}}{2\log_{10}(\Omega_s/\Omega_p)}\Big\rceil$ e $\Omega_0$ fra $\Omega_p\varepsilon^{-1/N}$ e $\Omega_s(A^2-1)^{-1/2N}$.Filtri di Butterworth →).