Esercizio - Progetto di un Butterworth con la bilineare

Testo (esempio 7.2 delle dispense del corso Multimedia Signal Processing, UniPD, e lezioni 22-23). Progettare un filtro IIR passa-basso di Butterworth con le specifiche: frequenza di campionamento $F_s=8$ kHz, limite della banda passante $f_p=1$ kHz, limite della banda oscura $f_s=2$ kHz, ondulazione in banda passante $R_p=0{,}1$ dB, attenuazione in banda oscura $R_s=40$ dB. Determinare l'ordine $N$, la pulsazione di taglio, i poli e i coefficienti di $H(z)$ e verificare le specifiche.

Teoria usata: 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 →, 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 →, 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 →, 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à →.

1. Dalle frequenze analogiche alle digitali, e viceversa (prewarping)

Le pulsazioni digitali normalizzate sono $\hat\omega=2\pi f/F_s$: $$\hat\omega_p=2\pi\frac{1000}{8000}=\frac\pi4,\qquad\hat\omega_s=2\pi\frac{2000}{8000}=\frac\pi2 .$$ La bilineare le mette in corrispondenza con $\Omega=\tan\frac{\hat\omega}2$, quindi il filtro analogico va progettato per $$\Omega_p=\tan\frac\pi8=0{,}4142,\qquad\Omega_s=\tan\frac\pi4=1 .$$ Gli errori sono conservati dalla trasformazione, quindi si usano direttamente: $$\varepsilon^2=10^{R_p/10}-1=10^{0{,}01}-1=0{,}02329\ \ (\varepsilon=0{,}1526),\qquad A=10^{R_s/20}=100 .$$

2. Ordine

$$N\ge\frac{\log_{10}\dfrac{A^2-1}{\varepsilon^2}}{2\log_{10}\dfrac{\Omega_s}{\Omega_p}}=\frac{\log_{10}\dfrac{9999}{0{,}02329}}{2\log_{10}\dfrac{1}{0{,}4142}}=\frac{5{,}6327}{0{,}7656}=7{,}358\ \Rightarrow\ N=8 .$$ Si arrotonda per eccesso: con $N=7$ non si rispetterebbero le specifiche.

3. Pulsazione di taglio $\Omega_0$

Due scelte possibili: $$\Omega_0=\Omega_s\,(A^2-1)^{-\frac1{2N}}=1\cdot9999^{-1/16}=0{,}5623\quad\text{(banda oscura esatta)},\qquad\Omega_0=\Omega_p\,\varepsilon^{-\frac1N}=0{,}4142\cdot0{,}1526^{-1/8}=0{,}5239\quad\text{(banda passante esatta)}.$$ Poiché $N=8>7{,}358$, ogni valore fra $0{,}5239$ e $0{,}5623$ è accettabile. Si sceglie $\Omega_0=0{,}5623$; il taglio digitale a $-3$ dB è $\hat\omega_0=2\arctan0{,}5623=0{,}326\pi$ (cioè $1{,}30$ kHz).

4. Poli analogici e digitali

I poli di $H_a(s)$ stanno sul semipiano sinistro, sulla circonferenza di raggio $\Omega_0$: $s_k=\Omega_0\big[-\sin\frac{(2k-1)\pi}{16}+j\cos\frac{(2k-1)\pi}{16}\big]$, $k=1,\dots,8$:

Tutti i poli digitali sono dentro il cerchio unitario (come garantito dalla bilineare). Esempio di calcolo per $k=1$: $s_1=-0{,}1097+0{,}5515j$, $1+s_1=0{,}8903+0{,}5515j$, $1-s_1=1{,}1097-0{,}5515j$, il rapporto dà $0{,}4453+0{,}7183j$.

5. Funzione di sistema

$$H(z)=A\,\frac{(1+z^{-1})^8}{\prod_{k=1}^8(1-p_kz^{-1})},\qquad A=\frac{\prod_k(1-p_k)}{2^8}=6{,}1595\cdot10^{-4},$$ dove $A$ si ricava da $H(1)=1$. Svolgendo i prodotti: $$H(z)=\frac{6{,}1595\cdot10^{-4}\,\big(1+8z^{-1}+28z^{-2}+56z^{-3}+70z^{-4}+56z^{-5}+28z^{-6}+8z^{-7}+z^{-8}\big)}{1-2{,}76773z^{-1}+4{,}16903z^{-2}-3{,}91878z^{-3}+2{,}48927z^{-4}-1{,}06826z^{-5}+0{,}30055z^{-6}-0{,}05019z^{-7}+0{,}00379z^{-8}}.$$ I coefficienti binomiali $1,8,28,56,70,\dots$ sono quelli di $(1+z^{-1})^8$: otto zeri in $z=-1$. Numeratore: $0{,}00062;\,0{,}00493;\,0{,}01725;\,0{,}03449;\,0{,}04312;\dots$ (simmetrico). Gli stessi coefficienti escono da scipy.signal.butter(8, 0.32612) ($0{,}32612=\frac{2}{\pi}\arctan0{,}5623$).

6. Verifica delle specifiche

import numpy as np
from scipy import signal
O0 = 0.5623448400103591
b, a = signal.butter(8, 2*np.arctan(O0)/np.pi)           # Butterworth digitale con taglio 2 arctan(O0)
w, H = signal.freqz(b, a, worN=[np.pi/4, np.pi/2])
print(20*np.log10(np.abs(H)))                            # -0.0325 dB in wp, -40.00 dB in ws
print(signal.buttord(0.25, 0.5, 0.1, 40))                # (8, 0.30724...): N = 8

Risultati: con $\Omega_0=0{,}5623$ il modulo vale $0{,}99627$ ($-0{,}0325$ dB) in $\hat\omega_p$, quindi la banda passante è rispettata con margine (serve almeno $-0{,}1$ dB), e $0{,}0100$ ($-40{,}00$ dB) in $\hat\omega_s$: la banda oscura è soddisfatta con uguaglianza. Con $\Omega_0=0{,}5239$ (taglio digitale $0{,}3072\pi$, valore restituito da buttord) si ha esattamente $-0{,}1000$ dB in $\hat\omega_p$ e $-44{,}92$ dB in $\hat\omega_s$. Con il valore intermedio $\Omega_0=0{,}5431$: $-0{,}0565$ dB e $-42{,}41$ dB. In tutti e tre i casi le specifiche sono soddisfatte, come previsto.

Osservazione. Il filtro ha $8$ poli e $8$ zeri in $z=-1$: costa $17$ moltiplicazioni per campione (o meno sfruttando gli zeri). Un FIR a fase lineare che realizzi le stesse specifiche richiederebbe un ordine molto più alto. Nelle dispense (esempio 7.2) il parametro è scritto $\varepsilon=\sqrt{10^{R_p/10}}\simeq0{,}1526$: manca il "$-1$" sotto la radice, perché $\sqrt{10^{0{,}01}}=1{,}0116$. Il valore $0{,}1526$ è quello di $\varepsilon=\sqrt{10^{R_p/10}-1}$, cioè $\varepsilon^2=10^{0{,}01}-1=0{,}02329$, ed è l'unico coerente con $|H_a(j\Omega_p)|^2=\frac1{1+\varepsilon^2}$. Il risultato $N\simeq7{,}3576$ delle dispense coincide con il nostro $7{,}358$ a meno dell'arrotondamento.