Filtri di Butterworth

Per usare la Trasformazione bilinearePer progettare un IIR si parte da un filtro analogico noto $H_a(s)$ e si sostituisce $s=\frac{z-1}{z+1}$ (trasformazione bilineare): $H(z)=H_a!\big(\frac{z-1}{z+1}\big)$. Questa mappa manda funzioni razionali in funzioni razionali, l'asse immaginario $s=j\Omega$ nella circonferenza unitaria $z=e^{j\hat\omega}$ e il semipiano sinistro nel disco unitario, quindi conserva la stabilità. Poli e zeri vanno in $\hat z=\frac{1+\hat s}{1-\hat s}$ ($s=0\to z=1$, $s=\infty\to z=-1$). Le frequenze si corrispondono con $\Omega=\tan\frac{\hat\omega}{2}$, cioè $\hat\omega=2\arctan\Omega$: relazione non lineare che comprime l'asse e distorce le frequenze (warping). Il modulo si conserva (oscillazioni, tolleranze $\delta_p,\delta_s$) ma la fase no: si perde la fase lineare. Per i filtri selettivi in frequenza la distorsione si compensa progettando il filtro analogico alle frequenze trasformate $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$ (prewarping).Trasformazione bilineare → serve un prototipo analogico passa-basso. I due grandi gruppi di soluzioni analogiche sono quelli tutti-poli (all-pole) e quelli con zeri finiti. Il filtro di Butterworth è il più semplice del primo gruppo e il principale del corso.

Specifiche analogiche e famiglie di filtri

Per un passa-basso analogico passivo ($\max|H_a|=1$) le specifiche sono le stesse del caso digitale (Progetto di filtri FIR con il metodo delle finestreUn filtro FIR a fase lineare di tipo I (ordine $N$ pari, $h[n]=h[N-n]$) ha risposta $H(e^{j\hat\omega})=e^{-j\hat\omega N/2}\bar H(\hat\omega)$ con ampiezza $\bar H(\hat\omega)=\sum_{n=0}^{N/2}p_n\cos(n\hat\omega)$. Per approssimare un'ampiezza desiderata $D(\hat\omega)$ (passa-basso: $1$ in banda passante, $0$ in banda oscura, con tolleranze $\delta_p,\delta_s$ e frequenze $\hat\omega_p,\hat\omega_s$) ci sono tre metodi. Finestre: si tronca la serie di Fourier di $D$, cioè $h[n]=h_d[n-N/2],w[n]$ con $h_d[n]=\frac{\hat\omega_0}{\pi}\operatorname{sinc}\frac{\hat\omega_0 n}{\pi}$; la rettangolare è ottima in errore quadratico ma dà il fenomeno di Gibbs (9%, $21$ dB), le finestre rastremate (Hann, Hamming, Blackman, Kaiser) abbassano i lobi laterali allargando la transizione ($\hat\omega_s-\hat\omega_p\simeq\alpha,2\pi/L$). Campionamento in frequenza: $h=\mathrm{IDFT}$ dei campioni di $D$, esatto solo sui campioni. Minimax (Parks-McClellan): errore pesato minimo nel caso peggiore, soluzione equiripple con almeno $r+2$ alternanze, ordine più basso a parità di specifiche.Progetto di filtri FIR con il metodo delle finestre →) con pulsazioni analogiche $\Omega_p<\Omega_s$ ed errori $\delta_p,\delta_s$. In ambito analogico si preferisce chiamare

• $\dfrac1A=\delta_s$, con $A$ detta attenuazione in banda oscura;
• $\dfrac1{\sqrt{1+\varepsilon^2}}=1-\delta_p$, con $\varepsilon$ legato all'ondulazione in banda passante.

Per i due valori in decibel valgono $$R_p=20\log_{10}\frac{1}{1-\delta_p}=10\log_{10}(1+\varepsilon^2),\qquad R_s=-20\log_{10}\delta_s=20\log_{10}A,$$ cioè $\varepsilon^2=10^{R_p/10}-1$ e $A=10^{R_s/20}$. Le due condizioni sul modulo quadro sono $|H_a(j\Omega_p)|^2\ge\frac1{1+\varepsilon^2}$ e $|H_a(j\Omega_s)|^2\le\frac1{A^2}$.

Esempio. $R_p=0{,}1$ dB e $R_s=40$ dB: $\varepsilon^2=10^{0{,}01}-1=0{,}02329$ ($\varepsilon=0{,}1526$) e $A=100$.

Famiglie di filtri analogici passa-basso:

• tutti-poli: $H_a(s)=\frac1{D_a(s)}$, con $D_a$ polinomio; tutti gli zeri sono all'infinito, quindi $H_a(j\Omega)\to0$ per $\Omega\to\infty$. Sono i filtri di Butterworth e di Chebyshev di tipo 1;
• con zeri finiti: $H_a(s)=\frac{N_a(s)}{D_a(s)}$ con $\deg N_a<\deg D_a$: i Chebyshev di tipo 2 e gli ellittici (o di Cauer).

Il modulo quadro di Butterworth

Il filtro di Butterworth è definito dalla massima piattezza del modulo in $\Omega=0$. Si lavora sul modulo quadro, che per un filtro reale è una funzione pari di $\Omega$: $$M_a(\Omega)=|H_a(j\Omega)|^2=\frac{1}{P_N(\Omega^2)},\qquad P_N(\Omega^2)=p_0+p_1\Omega^2+\dots+p_N\Omega^{2N}.$$ Solo potenze pari: se $H_a$ è reale, $M_a=H_a(j\Omega)H_a^*(j\Omega)$ è pari. Tre vincoli, nell'ordine:

1. Guadagno unitario in continua: $M_a(0)=1$, quindi $P_N(0)=p_0=1$.
2. Massima piattezza in $0$: $M_a$ deve valere $1$ il più a lungo possibile vicino all'origine, cioè $P_N(\Omega^2)-1$ deve essere un infinitesimo di ordine massimo in $\Omega=0$. Per $\Omega\ll1$, $3\Omega^2$ è molto più grande di $5\Omega^8$: conviene quindi annullare tutti i coefficienti bassi, $p_1=\dots=p_{N-1}=0$, e lasciare solo l'ultimo: $P_N(\Omega^2)=1+p_N\Omega^{2N}$.
3. Vincolo sul taglio: $|H_a(j\Omega_0)|^2=\frac12$ (cioè $-3$ dB, modulo $\frac1{\sqrt2}\simeq0{,}707$: il segnale esce al $70{,}7\%$ dell'ampiezza). Da $\frac1{1+p_N\Omega_0^{2N}}=\frac12$ segue $p_N=\Omega_0^{-2N}$.

Formula (Butterworth analogico). $$|H_a(j\Omega)|^2=\frac{1}{1+\big(\frac{\Omega}{\Omega_0}\big)^{2N}}$$ con $N$ ordine e $\Omega_0$ pulsazione di taglio a $-3$ dB. Vale $1$ in $0$, $\frac12$ in $\Omega_0$ per ogni $N$, ed è monotona decrescente fino a $0$ per $\Omega\to\infty$; per $N\to\infty$ tende al rettangolo ideale.

Esempio. Con $\Omega_0=1$ e $N=4$: $|H_a(j0{,}5)|^2=\frac1{1+0{,}5^8}=0{,}9961$, $|H_a(j2)|^2=\frac1{1+2^8}=0{,}00389$ (circa $-24$ dB). A $\Omega=2\Omega_0$ (un'ottava sopra il taglio) l'attenuazione è $10\log_{10}(1+2^{2N})\simeq6N$ dB: $24$ dB per $N=4$, $48$ dB per $N=8$.

Dal modulo alla funzione di trasferimento

Il modulo quadro non dà la fase. Si cerca la funzione razionale stabile $H_a(s)$ tale che $H_a(j\Omega)H_a(-j\Omega)=M_a(\Omega)$ (per $H_a$ reale $H_a^*(j\Omega)=H_a(-j\Omega)$). Estendendo a tutto il piano con $j\Omega=s$: $$M_a(s)=H_a(s)\,H_a(-s)=\frac{1}{1+\big(\frac{s}{j\Omega_0}\big)^{2N}} .$$ I poli di $M_a(s)$ sono le soluzioni di $\big(\frac{s}{j\Omega_0}\big)^{2N}=-1=e^{j(\pi+2k\pi)}$, cioè $\frac{s}{j\Omega_0}=e^{j\frac{(2k-1)\pi}{2N}}$ e $$s_k=j\Omega_0\,e^{j\frac{(2k-1)\pi}{2N}}=\Omega_0\,e^{j\pi\left(\frac12+\frac{2k-1}{2N}\right)},\qquad k=1,\dots,2N.$$ Sono $2N$ poli sulla circonferenza di raggio $\Omega_0$, simmetrici rispetto agli assi reale e immaginario (se $s_k$ è un polo, lo sono anche $-s_k$ e $s_k^*$): metà nel semipiano sinistro, metà nel destro. Per ottenere un sistema causale e stabile (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à →) $H_a(s)$ prende i $N$ poli con $\operatorname{Re}s_k<0$ (gli altri appartengono a $H_a(-s)$). Con la formula $\cos(\frac\pi2+x)=-\sin x$ e $\sin(\frac\pi2+x)=\cos x$ si scrivono

$$s_k=\Omega_0\Big[-\sin\frac{(2k-1)\pi}{2N}+j\cos\frac{(2k-1)\pi}{2N}\Big],\quad k=1,\dots,N,\qquad H_a(s)=\frac{K}{\prod_{k=1}^{N}(s-s_k)} .$$

La costante $K$ si fissa con $H_a(0)=1$: $K=\prod_k(-s_k)=\Omega_0^N$ (i poli hanno tutti modulo $\Omega_0$ e vengono in coppie coniugate, quindi il prodotto è reale e positivo). Dunque $H_a(s)=\dfrac{\Omega_0^N}{\prod_{k=1}^N(s-s_k)}$.

Esempio ($\Omega_0=1$). $N=2$: $s_{1,2}=\frac{-1\pm j}{\sqrt2}$, $H_a(s)=\frac1{s^2+\sqrt2\,s+1}$. $N=3$: $s=-1$ e $-\frac12\pm j\frac{\sqrt3}2$, $H_a(s)=\frac1{(s+1)(s^2+s+1)}=\frac1{s^3+2s^2+2s+1}$. $N=4$: gli otto poli sono ai $22{,}5^\circ+k\cdot45^\circ$ e $H_a$ tiene i quattro con parte reale negativa.

Un filtro di Butterworth è completamente specificato da due parametri: $N$ e $\Omega_0$.

Il filtro digitale

Con la bilineare $s=\frac{z-1}{z+1}$ (Trasformazione bilinearePer progettare un IIR si parte da un filtro analogico noto $H_a(s)$ e si sostituisce $s=\frac{z-1}{z+1}$ (trasformazione bilineare): $H(z)=H_a!\big(\frac{z-1}{z+1}\big)$. Questa mappa manda funzioni razionali in funzioni razionali, l'asse immaginario $s=j\Omega$ nella circonferenza unitaria $z=e^{j\hat\omega}$ e il semipiano sinistro nel disco unitario, quindi conserva la stabilità. Poli e zeri vanno in $\hat z=\frac{1+\hat s}{1-\hat s}$ ($s=0\to z=1$, $s=\infty\to z=-1$). Le frequenze si corrispondono con $\Omega=\tan\frac{\hat\omega}{2}$, cioè $\hat\omega=2\arctan\Omega$: relazione non lineare che comprime l'asse e distorce le frequenze (warping). Il modulo si conserva (oscillazioni, tolleranze $\delta_p,\delta_s$) ma la fase no: si perde la fase lineare. Per i filtri selettivi in frequenza la distorsione si compensa progettando il filtro analogico alle frequenze trasformate $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$ (prewarping).Trasformazione bilineare →) si ottiene, con $N$ zeri all'infinito che diventano uno zero di molteplicità $N$ in $z=-1$, $$H(z)=A\,\frac{(1+z^{-1})^N}{\prod_{k=1}^{N}(1-p_kz^{-1})},\qquad p_k=\frac{1+s_k}{1-s_k},\quad A=\frac{\prod_k(1-p_k)}{2^N},$$ dove $A$ è il guadagno scelto in modo che $H(1)=H(e^{j0})=1$ (infatti $H(1)=A\,2^N/\prod(1-p_k)$). I poli $p_k$ giacciono su una circonferenza del piano $z$ che è l'immagine della circonferenza di raggio $\Omega_0$ del piano $s$, e sono dentro il cerchio unitario. Il modulo digitale si ricava sostituendo $\Omega=\tan\frac{\hat\omega}2$:

$$|H(e^{j\hat\omega})|^2=\frac1{1+\Big(\dfrac{\tan(\hat\omega/2)}{\Omega_0}\Big)^{2N}} .$$

Il taglio digitale a $-3$ dB è $\hat\omega_0=2\arctan\Omega_0$: per $\Omega_0=1$ è $\hat\omega_0=\pi/2$ e le curve per ogni $N$ passano tutte per $0{,}707$ in $\pi/2$ (verificato numericamente con scipy.signal.butter: la differenza con la formula è dell'ordine di $10^{-15}$).

Esempio ($\Omega_0=1$, $N=4$). $s_k=e^{j\pi(\frac12+\frac{2k-1}8)}$ con $k=1,2,3,4$ (angoli $112{,}5^\circ$, $157{,}5^\circ$, $202{,}5^\circ$, $247{,}5^\circ$). I poli digitali sono $p_k=\frac{1+s_k}{1-s_k}$, tutti dentro $|z|<1$; lo zero di molteplicità $4$ sta in $z=-1$.

Scelta dell'ordine dalle specifiche

Dato un passa-basso digitale con $\hat\omega_p,\hat\omega_s,R_p,R_s$, si trasformano le frequenze con $\Omega_p=\tan\frac{\hat\omega_p}2$, $\Omega_s=\tan\frac{\hat\omega_s}2$ e si impongono le due condizioni sul modulo quadro con il segno di uguaglianza: $$\frac1{1+\big(\frac{\Omega_p}{\Omega_0}\big)^{2N}}=\frac1{1+\varepsilon^2},\qquad\frac1{1+\big(\frac{\Omega_s}{\Omega_0}\big)^{2N}}=\frac1{A^2},$$ equivalenti a $\big(\frac{\Omega_p}{\Omega_0}\big)^{2N}=\varepsilon^2$ e $\big(\frac{\Omega_s}{\Omega_0}\big)^{2N}=A^2-1$. Dividendo la seconda per la prima si elimina $\Omega_0$: $\big(\frac{\Omega_s}{\Omega_p}\big)^{2N}=\frac{A^2-1}{\varepsilon^2}$, da cui

$$\boxed{N\ \ge\ \frac{\log_{10}\dfrac{A^2-1}{\varepsilon^2}}{2\log_{10}\dfrac{\Omega_s}{\Omega_p}}}$$

e si prende il più piccolo intero maggiore o uguale. L'ordine dipende da $\varepsilon$, da $A$ e dal rapporto $\Omega_s/\Omega_p$ (più la transizione è stretta, più alto è l'ordine). Poi $\Omega_0$ si sceglie con una delle due condizioni: $$\Omega_0=\Omega_s\,(A^2-1)^{-\frac1{2N}}\quad\text{(banda oscura soddisfatta con uguaglianza)},\qquad\Omega_0=\Omega_p\,\varepsilon^{-\frac1N}\quad\text{(banda passante con uguaglianza)} .$$ Poiché $N$ è stato arrotondato per eccesso, una delle due condizioni sarà soddisfatta con margine; ogni valore di $\Omega_0$ fra i due soddisfa entrambe.

Esempio (specifiche delle dispense). $F_s=8$ kHz, $f_p=1$ kHz, $f_s=2$ kHz, $R_p=0{,}1$ dB, $R_s=40$ dB. $\hat\omega_p=\frac\pi4$, $\hat\omega_s=\frac\pi2$, $\Omega_p=\tan\frac\pi8=0{,}4142$, $\Omega_s=\tan\frac\pi4=1$, $\varepsilon=0{,}1526$, $A=100$: $N\ge\frac{\log_{10}(9999/0{,}02329)}{2\log_{10}(1/0{,}4142)}=\frac{5{,}6327}{0{,}7656}=7{,}358$, quindi $N=8$. I due valori di $\Omega_0$ sono $0{,}5623$ e $0{,}5239$. Il progetto completo (poli e coefficienti) è in Esercizio - Progetto di un Butterworth con la bilineare.

Fase e ritardo

Come tutti i filtri IIR il Butterworth non ha fase lineare: il suo ritardo di gruppo cresce vicino al taglio (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 →). Se serve una transizione ripida con ordine basso si usano Chebyshev o ellittici, che ammettono ondulazioni (sono in tabella sopra) e a parità di specifiche hanno ordini molto più bassi: per le specifiche $\hat\omega_p=0{,}25\pi$, $\hat\omega_s=0{,}375\pi$, $R_p=0{,}055$ dB, $A_s=50$ dB si ottengono ordine $17$ (Butterworth), $9$ (Chebyshev 1), $6$ (ellittico).

Domande d'esame

1. Ricava la funzione di trasferimento del filtro di Butterworth analogico di ordine $N$ e discuti la posizione dei poli. Traccia: modulo quadro $\frac1{P_N(\Omega^2)}$ pari, $P_N(0)=1$, massima piattezza $\Rightarrow P_N=1+p_N\Omega^{2N}$, taglio a $-3$ dB $\Rightarrow p_N=\Omega_0^{-2N}$; $M_a(s)=H_a(s)H_a(-s)$, poli $s_k=\Omega_0e^{j\pi(1/2+(2k-1)/(2N))}$ su una circonferenza di raggio $\Omega_0$; si tengono quelli con $\operatorname{Re}s<0$, $K=\Omega_0^N$.
2. Come si determina l'ordine di un filtro di Butterworth digitale dalle specifiche? Traccia: prewarping $\Omega=\tan\frac{\hat\omega}2$, $\varepsilon$ e $A$ da $R_p,R_s$, formula di $N$ con arrotondamento per eccesso, scelta di $\Omega_0$ e perché una delle due condizioni ha margine.
3. Descrivi il progetto di un passa-basso digitale di Butterworth con la bilineare e le sue proprietà. Traccia: poli e zero in $-1$ di molteplicità $N$, $H(1)=1$; $|H|^2=\frac1{1+(\tan(\hat\omega/2)/\Omega_0)^{2N}}$; monotono, $-3$ dB in $2\arctan\Omega_0$; fase non lineare.