Esercizio - Frequenze bloccate e uscita a un coseno di un IIR

Testo (esercitazione del corso Multimedia Signal Processing, UniPD, lezione 24, esercizio 2.5 sui filtri IIR). Un sistema è definito dall'operatore $$H(z)=(1-z^{-1})\,\frac{1+z^{-4}}{1-0{,}8z^{-2}} .$$

1. Scrivere la descrizione nel tempo come equazione alle differenze.
2. Scrivere la formula della risposta in frequenza.
3. Il sistema può "bloccare" certi ingressi: determinare tutte le frequenze $\hat\omega_0$ per cui la risposta a $x[n]=\cos(\hat\omega_0n)$ è nulla.
4. Quando l'ingresso è $x[n]=\cos(\pi n/3)$ determinare l'uscita nella forma $y[n]=A\cos(\hat\omega_0n+\phi)$, con valori numerici di $A$, $\hat\omega_0$, $\phi$.

Teoria usata: 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à →, 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 →, Risposta in frequenza dei sistemi FIRSe all'ingresso di un FIR c'è un esponenziale complesso A e^{jφ} e^{jω̂n} (per ogni n), l'uscita è lo stesso esponenziale moltiplicato per H(ω̂) = Σ b_k e^{-jω̂k}: la frequenza non cambia, ampiezza e fase sono modificate da |H| (guadagno) e ∠H (sfasamento). Per sovrapposizione si trattano somme di sinusoidi. H è periodica di periodo 2π e, per coefficienti reali, hermitiana (|H| pari, fase dispari). La cascata ha H = H1·H2. Esempi: ritardo (fase lineare), differenza prima (passa-alto), {1,2,1} (passa-basso), media mobile di L punti (Dirichlet: |H| = |sin(Lω̂/2)/(L sin(ω̂/2))|, fase lineare -(L-1)ω̂/2).Risposta in frequenza dei sistemi FIR →, Filtri notch e applicazioni dei filtri FIRUno zero di H(z) sulla circonferenza unitaria in z0 = e^{j w0} annulla (a regime) le sinusoidi alla pulsazione w0. Per eliminare un coseno servono i due zeri coniugati: H(z) = (1 - e^{jw0} z^-1)(1 - e^{-jw0} z^-1) = 1 - 2 cos(w0) z^-1 + z^-2, cioè h = {1, -2cos w0, 1}, un FIR simmetrico di tipo I (ritardo 1). Si normalizza dividendo per 2 - 2cos(w0) per avere guadagno 1 in continua; più toni si eliminano mettendo in cascata (convoluzione) un filtro per ogni tono. Segue una panoramica delle applicazioni dei FIR (equalizzazione audio), vantaggi (fase lineare, stabilità) e costo (N moltiplicazioni per campione, molte più di un IIR).Filtri notch e applicazioni dei filtri FIR →, Numeri complessi, formula di Eulero ed esponenziali complessiUn numero complesso si scrive in forma cartesiana a+jb o polare |x|e^{jφ}; il prodotto moltiplica i moduli e somma le fasi. La formula di Eulero e^{jα}=cos α + j sin α lega esponenziali e sinusoidi e permette di trattare tutti i segnali del corso come somme di esponenziali complessi e^{(σ+jω)t}.Numeri complessi, formula di Eulero ed esponenziali complessi →. I suggerimenti del corso: ricondurre $H$ a un rapporto di polinomi in $z^{-1}$; per le frequenze bloccate cercare gli zeri sulla circonferenza unitaria; per l'uscita a un coseno calcolare modulo e fase di $H$ alla frequenza dell'ingresso.

1. Equazione alle differenze

Si moltiplicano i due fattori del numeratore: $(1-z^{-1})(1+z^{-4})=1-z^{-1}+z^{-4}-z^{-5}$. Da $Y(1-0{,}8z^{-2})=(1-z^{-1}+z^{-4}-z^{-5})X$: $$y[n]=0{,}8\,y[n-2]+x[n]-x[n-1]+x[n-4]-x[n-5].$$ È un IIR di ordine $2$ ($N=2$ coefficienti di retroazione, di cui solo $a_2=0{,}8$ diverso da zero) con $M=5$.

2. Risposta in frequenza

$$H(e^{j\hat\omega})=\frac{(1-e^{-j\hat\omega})(1+e^{-j4\hat\omega})}{1-0{,}8\,e^{-j2\hat\omega}} .$$ I poli sono $z^2=0{,}8$, cioè $z=\pm0{,}894$ (dentro il cerchio unitario: il sistema è stabile e la risposta in frequenza esiste).

3. Frequenze bloccate

Il coseno $\cos(\hat\omega_0n)$ viene annullato se $H(e^{j\hat\omega_0})=0$. Il denominatore non si annulla mai sulla circonferenza (poli a $\pm0{,}894$, e infatti $|1-0{,}8e^{-j2\hat\omega}|^2=1{,}64-1{,}6\cos2\hat\omega\ge0{,}04$), quindi bastano gli zeri del numeratore. Per il modulo quadro si usa la regola $|1+ae^{-jN\hat\omega}|^2=1+a^2+2a\cos(N\hat\omega)$: $$|H(e^{j\hat\omega})|^2=\frac{(2-2\cos\hat\omega)\,(2+2\cos4\hat\omega)}{1{,}64-1{,}6\cos2\hat\omega}.$$ Gli zeri sono le radici di $(1-z^{-1})(1+z^{-4})$:

• $1-z^{-1}=0\Rightarrow z=1$, cioè $\hat\omega=0$;
• $1+z^{-4}=0\Rightarrow z^4=-1=e^{j\pi}=e^{j(\pi+2\pi l)}\Rightarrow z=e^{j(\frac\pi4+\frac{\pi}2l)}$, $l=0,1,2,3$, cioè $\hat\omega=\frac\pi4,\ \frac{3\pi}4,\ \frac{5\pi}4,\ \frac{7\pi}4$.

I cinque zeri sono tutti sulla circonferenza unitaria. In $[-\pi,\pi]$ le frequenze bloccate sono $$\hat\omega_0=0,\ \pm\frac\pi4,\ \pm\frac{3\pi}4$$ (le frequenze negative sono equivalenti perché $\cos(-x)=\cos x$, e $\frac{5\pi}4,\frac{7\pi}4$ coincidono con $-\frac{3\pi}4,-\frac\pi4$ a meno di $2\pi$). Nota: la soluzione del corso scrive "$0,\frac\pi4,\pm\frac{3\pi}4$", ma il $\pm$ vale anche per $\frac\pi4$.

4. Uscita per $x[n]=\cos(\pi n/3)$

Per un sistema LTI stabile, $y[n]=|H(e^{j\hat\omega_0})|\cos\big(\hat\omega_0n+\angle H(e^{j\hat\omega_0})\big)$ con $\hat\omega_0=\frac\pi3$. Si calcolano i tre fattori di $H(e^{j\pi/3})$:

• $1-e^{-j\pi/3}=1-\big(\frac12-j\frac{\sqrt3}2\big)=\frac12+j\frac{\sqrt3}2=e^{j\pi/3}$ (modulo $1$, fase $\frac\pi3$);
• $1+e^{-j4\pi/3}=1+\big(-\frac12+j\frac{\sqrt3}2\big)=\frac12+j\frac{\sqrt3}2=e^{j\pi/3}$ (modulo $1$, fase $\frac\pi3$), perché $e^{-j4\pi/3}=\cos\frac{4\pi}3-j\sin\frac{4\pi}3=-\frac12+j\frac{\sqrt3}2$;
• $1-0{,}8e^{-j2\pi/3}=1-0{,}8\big(-\frac12-j\frac{\sqrt3}2\big)=1{,}4+j\,0{,}4\sqrt3$, di modulo $\sqrt{1{,}4^2+(0{,}4\sqrt3)^2}=\sqrt{2{,}44}=1{,}562$ e fase $\arctan\frac{0{,}4\sqrt3}{1{,}4}=0{,}4595=0{,}1463\pi$.

Il modulo è il rapporto dei moduli, la fase è la somma delle fasi del numeratore meno quella del denominatore: $$|H|=\frac{1\cdot1}{1{,}562}=0{,}640,\qquad\angle H=\frac\pi3+\frac\pi3-0{,}4595=2{,}0944-0{,}4595=1{,}635\ \text{rad}=0{,}5204\,\pi .$$ Quindi $$y[n]=0{,}640\,\cos\!\Big(\frac\pi3n+0{,}52\pi\Big),\qquad A=0{,}64,\ \ \hat\omega_0=\frac\pi3,\ \ \phi=0{,}52\pi\simeq1{,}635\ \text{rad}.$$ Verifica con freqz ($|H|=0{,}64018$, $\angle H=1{,}63486$ rad) e con una simulazione: filtrando $\cos(\pi n/3)$ con lfilter per $400$ campioni, gli ultimi valori ($0{,}5738;\,0{,}5328;\,-0{,}0410;\,-0{,}5738;\dots$) coincidono con la formula.