Esercizio - Laboratorio 5 - notch IIR e vuvuzela

Testo (laboratorio 5 del corso Multimedia Signal Processing, UniPD, lezione 23; con codice Python numpy/scipy). Una registrazione di una partita dei Mondiali del 2010 è disturbata dal suono della vuvuzela, un corno di plastica che emette una nota monotona a circa $235$ Hz (prima armonica) con le armoniche successive a circa $465$, $705$, ... Hz, e impedisce di sentire il commentatore. Con un filtro IIR a tacca (notch) si deve eliminare il disturbo:

1. caricare il segnale, tracciarlo nel tempo e in frequenza (modulo in dB, asse in Hz) e verificare le armoniche;
2. progettare un notch IIR del secondo ordine per la prima armonica, per $r\in[0{,}95,\,0{,}99]$: tracciare modulo e fase e il diagramma poli-zeri;
3. filtrare il segnale e confrontare gli spettri; ascoltare;
4. mettere in cascata più notch per togliere anche le armoniche successive;
5. (facoltativo) confrontare il ritardo di gruppo del notch IIR con quello del notch FIR del laboratorio 2.

Dati. Il file vuvuzela.wav sta nella cartella del laboratorio 5 su Moodle. È mono, a $16$ bit, con $F_s=48000$ Hz, $500\,000$ campioni ($10{,}4$ s); in scala $[-1,1)$ ha picco $0{,}373$ e valore efficace $0{,}0611$. I risultati sotto sono ottenuti eseguendo il codice su questo file. Teoria usata: 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à →, 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 →, 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 →, Campionamento e ricostruzioneIl campionamento $s_c(nT)=s(nT)$ trasforma un segnale continuo in uno discreto e, in frequenza, ripete lo spettro con periodo $F_c=1/T$: $S_c(f)=\sum_kS(f-kF_c)$. Se le repliche si sovrappongono si ha aliasing e il segnale non è recuperabile. Un interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g$ produce $\tilde s(t)=\sum_nT g(t-nT)s(nT)$ e $\tilde S=G,S_c$. Teorema del campionamento: se $s$ ha banda $B$ e $F_c\ge2B$, l'interpolatore ideale $G=\operatorname{rect}(f/F_c)$ ricostruisce esattamente $s(t)=\sum s(nT)\operatorname{sinc}(F_c(t-nT))$. Se le ipotesi non valgono c'è un errore, in banda e fuori banda, riducibile con un prefiltro anti-aliasing.Campionamento e ricostruzione →.

Il filtro notch del secondo ordine

Per cancellare la frequenza $f_0$ si pone $\hat\omega_0=2\pi f_0/F_s$ (con $f_0=235$ Hz e $F_s=48$ kHz: $\hat\omega_0=0{,}03076=0{,}00979\pi$). Il notch IIR ha due zeri sulla circonferenza unitaria in $e^{\pm j\hat\omega_0}$ (annullano esattamente $\hat\omega_0$) e due poli coniugati vicini, in $re^{\pm j\hat\omega_0}$ con $r<1$ ($r$ vicino a $1$), che fanno tornare il modulo a $1$ appena ci si allontana da $\hat\omega_0$: $$H(z)=b_0\,\frac{(1-e^{j\hat\omega_0}z^{-1})(1-e^{-j\hat\omega_0}z^{-1})}{(1-re^{j\hat\omega_0}z^{-1})(1-re^{-j\hat\omega_0}z^{-1})}=b_0\,\frac{1-2\cos\hat\omega_0\,z^{-1}+z^{-2}}{1-2r\cos\hat\omega_0\,z^{-1}+r^2z^{-2}}.$$ Con $b_0=\dfrac{1-2r\cos\hat\omega_0+r^2}{2-2\cos\hat\omega_0}$ si ottiene $H(1)=1$ (guadagno $0$ dB in continua). A differenza del notch FIR con soli zeri (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 →), che ha una tacca larga e altera il modulo anche lontano da $\hat\omega_0$, qui i poli vicini agli zeri "compensano" lo zero appena ci si allontana dalla frequenza eliminata, e la tacca resta stretta. La larghezza della tacca a $-3$ dB è circa $2(1-r)$ rad, cioè $2(1-r)F_s/2\pi$ Hz, purché sia molto minore di $\hat\omega_0$: se $1-r$ è confrontabile con $\hat\omega_0$ la stima non vale più (vedi sotto, $r=0{,}95$).

Le equazioni alle differenze sono $y[n]=2r\cos\hat\omega_0\,y[n-1]-r^2y[n-2]+b_0\big(x[n]-2\cos\hat\omega_0\,x[n-1]+x[n-2]\big)$. Il filtro è stabile perché i poli hanno modulo $r<1$.

Codice

import numpy as np
from scipy import signal
from scipy.io import wavfile
import matplotlib.pyplot as plt

Fs, y = wavfile.read('vuvuzela.wav')                    # int16, mono, Fs = 48000
x = y / 32768.0
f_h = np.array([235, 465, 694, 932, 1160, 1389, 1632.5, 1864, 2115, 2325, 2568, 2784])  # armoniche del corso (Hz)

def notch(w0, r):
    """Notch del secondo ordine: zeri in e^{+-j w0}, poli in r e^{+-j w0}, guadagno 1 in continua."""
    b0 = (1 - 2 * r * np.cos(w0) + r ** 2) / (2 - 2 * np.cos(w0))
    return b0 * np.array([1, -2 * np.cos(w0), 1]), np.array([1, -2 * r * np.cos(w0), r ** 2])

def potenza(sig, f1, f2):                               # potenza spettrale tra f1 e f2 Hz (si scarta il transitorio, 0,5 s)
    s = sig[Fs // 2:]
    X = np.abs(np.fft.rfft(s * np.hanning(len(s)))) ** 2
    fr = np.fft.rfftfreq(len(s), 1 / Fs)
    return X[(fr >= f1) & (fr < f2)].sum()

# ---- 0) spettro: dove stanno le armoniche? (Welch: 65536 campioni per segmento, risoluzione 0,73 Hz)
f, P = signal.welch(x, Fs, nperseg=2 ** 16, noverlap=2 ** 15)
for fh in f_h:
    m = (f > fh - 25) & (f < fh + 25)
    print(f"armonica del corso {fh:7.1f} Hz: picco a {f[m][np.argmax(P[m])]:7.1f} Hz, livello {10*np.log10(P[m].max()):.1f} dB")

# ---- 1) un solo notch sulla prima armonica, per due valori di r ----
w0 = 2 * np.pi * f_h[0] / Fs
for r in (0.95, 0.99):
    b, a = notch(w0, r)
    w, H = signal.freqz(b, a, worN=2000000); mag = np.abs(H)
    sotto = np.where(mag < mag.max() / np.sqrt(2))[0]
    larg = (w[sotto.max()] - w[sotto.min()]) * Fs / (2 * np.pi)
    y1 = signal.lfilter(b, a, x)
    print(f"r={r}: b={np.round(b, 4)} a={np.round(a, 4)}")
    print(f"   |H(0)|={mag[0]:.3f}, max|H|={mag.max():.4f}, larghezza a -3 dB = {larg:.1f} Hz (stima 2(1-r)Fs/2pi = {2*(1-r)*Fs/(2*np.pi):.1f} Hz)")
    print("   attenuazione della prima armonica (220-250 Hz):", round(10 * np.log10(potenza(y1, 220, 250) / potenza(x, 220, 250)), 1), "dB")
    plt.plot(w * Fs / (2 * np.pi), 20 * np.log10(mag + 1e-12), label=f"r={r}")
plt.xlim(0, 1000); plt.ylim(-40, 14); plt.xlabel("f (Hz)"); plt.ylabel("dB"); plt.legend(); plt.grid(); plt.savefig("lab5_notch.png")

# ---- 2) cascata di 12 notch ----
fuori = (lambda fr: (fr >= 100) & (fr < 8000) & np.all(np.abs(fr[:, None] - f_h) > 60, axis=1))
for r in (0.95, 0.99, 0.998):
    yy = x.copy(); Htot = np.ones(400000, complex); ww = np.linspace(1e-4, np.pi, 400000)
    for f in f_h:
        b, a = notch(2 * np.pi * f / Fs, r)
        yy = signal.lfilter(b, a, yy)
        Htot *= signal.freqz(b, a, worN=ww)[1]
    arm = sum(potenza(yy, f - 20, f + 20) for f in f_h) / sum(potenza(x, f - 20, f + 20) for f in f_h)
    s = yy[Fs // 2:]; sx = x[Fs // 2:]
    fr = np.fft.rfftfreq(len(s), 1 / Fs); k = fuori(fr)
    Y = np.abs(np.fft.rfft(s * np.hanning(len(s)))) ** 2; Xs = np.abs(np.fft.rfft(sx * np.hanning(len(sx)))) ** 2
    print(f"r={r}: armoniche (+-20 Hz) {10*np.log10(arm):.1f} dB; potenza totale {10*np.log10(Y.sum()/Xs.sum()):.1f} dB; "
          f"fuori dalle armoniche {10*np.log10(Y[k].sum()/Xs[k].sum()):.2f} dB; guadagno massimo {20*np.log10(np.abs(Htot).max()):.2f} dB; "
          f"valore efficace {yy.std():.4f} (ingresso {x.std():.4f}), picco {np.abs(yy).max():.3f}")

# ---- 3) ritardo di gruppo: notch IIR contro notch FIR (zeri sul cerchio, b = [1, -2cos w0, 1]) ----
wg = np.linspace(0.0005, np.pi, 400000)
for r in (0.95, 0.99):
    b, a = notch(w0, r)
    _, gd = signal.group_delay((b, a), w=wg)
    i = np.argmax(gd)
    print(f"r={r}: ritardo di gruppo IIR: massimo {gd[i]:.1f} campioni a {wg[i]*Fs/(2*np.pi):.0f} Hz;"
          f" a 100 Hz {np.interp(2*np.pi*100/Fs, wg, gd):.1f}, a 1000 Hz {np.interp(2*np.pi*1000/Fs, wg, gd):.2f}, a 3000 Hz {np.interp(2*np.pi*3000/Fs, wg, gd):.2f}")
_, gdf = signal.group_delay((np.array([1, -2 * np.cos(w0), 1]), [1]), w=wg[wg > 0.3])
print("ritardo di gruppo del notch FIR (b=[1,-2cos w0,1]) lontano da w0:", np.round(gdf.min(), 3), "-", np.round(gdf.max(), 3), "campioni")

Risultati e spiegazione

Lo spettro della registrazione. Quasi tutta la potenza è tra $100$ e $1000$ Hz ($72{,}6\%$) e sopra $4$ kHz resta lo $0{,}7\%$. Le dodici armoniche della lista del corso, prese in una banda di $\pm20$ Hz, contengono il $53{,}5\%$ della potenza totale ($69{,}4\%$ con $\pm40$ Hz), la sola prima armonica ($220$-$250$ Hz) il $13{,}4\%$: il disturbo domina, e la voce sta sotto di esso. I picchi dello spettro (Welch, risoluzione $0{,}73$ Hz; livelli in dB rispetto al fondo scala) sono:

Nove armoniche su dodici coincidono con la lista del corso entro $0{,}5$ Hz; la seconda se ne scosta di $3{,}6$ Hz ($461{,}4$ contro $465$), l'ottava e la nona di $17{,}6$ e $18{,}8$ Hz ($1881{,}6$ e $2096{,}2$ contro $1864$ e $2115$ Hz): la lista è una lettura da grafico. Le armoniche si indeboliscono di circa $20$ dB dalla prima alla dodicesima, e la serie continua oltre i $3$ kHz (picchi a $3036$, $3286$, $3711$, $3943$ Hz, tra $-66$ e $-70$ dB). Ci sono anche picchi non armonici a $84$ e $348$ Hz.

La nota non è stabile. La vuvuzela non emette una frequenza esatta: la prima armonica misurata su finestre di $1$ s si sposta tra $226$ e $240$ Hz durante la registrazione, e la riga nello spettro dell'intero file ha larghezza di $9{,}5$ Hz a $-3$ dB e di $36$ Hz a $-10$ dB (da $213$ a $249$ Hz). Una tacca più stretta della riga non la copre tutta: la scelta di $r$ è un compromesso, e con un tono sinusoidale ideale lo zero lo cancellerebbe del tutto, mentre qui l'attenuazione resta di una ventina di dB.

Progetto. Per la prima armonica ($\hat\omega_0=0{,}03076$) con i due valori estremi (con $F_s=48$ kHz):

I poli sono in $0{,}94955\pm0{,}02922j$ ($r=0{,}95$) e $0{,}98953\pm0{,}03045j$ ($r=0{,}99$), gli zeri in $0{,}99953\pm0{,}03076j$ (sul cerchio unitario, a $\pm\hat\omega_0$): poli e zeri sono vicinissimi, quindi l'effetto del filtro è locale. La costante di tempo del transitorio è $\frac1{1-r}$ campioni: $20$ campioni ($0{,}4$ ms) contro $100$ ($2{,}1$ ms).

Nota sul guadagno. Il testo del laboratorio dice che $b_0$ è scelto per avere "guadagno massimo uno". In realtà la normalizzazione impone $|H(1)|=1$ in continua, e a $48$ kHz, con la tacca vicinissima alla continua ($0{,}0098\pi$), questo ha un effetto grande: il modulo sale molto sopra $1$ appena si esce dalla tacca. Il massimo vale $3{,}78$ ($+11{,}5$ dB, in $\hat\omega=\pi$) per $r=0{,}95$ e $1{,}107$ per $r=0{,}99$; nella cascata di $12$ notch i massimi diventano $+23{,}4$ dB e $+1{,}4$ dB. Si vede nel grafico: con $r=0{,}95$ il filtro amplifica di $8{,}5$ dB a $700$ Hz e di $10{,}5$ dB a $1150$ Hz.

Filtraggio con un solo notch. La potenza nella banda $220$-$250$ Hz (la prima armonica) scende di $25{,}8$ dB con $r=0{,}95$ e di $20{,}2$ dB con $r=0{,}99$; il picco dello spettro passa da $-44{,}5$ dB a $-69{,}8$ dB e a $-64{,}1$ dB. La prima armonica diventa molto più debole ma non sparisce (la riga è larga e oscilla), e le altre undici armoniche restano: il disturbo si sente ancora. Con $r=0{,}95$ la potenza totale del segnale sale di $8{,}3$ dB (amplificazione fuori dalla tacca), con $r=0{,}99$ è invariata ($-0{,}1$ dB).

Cascata. Si ripete con i dodici notch (uno per armonica) uno dopo l'altro: la funzione di sistema complessiva è il prodotto delle $12$ funzioni di sistema (ordine $24$), e le tacche si sommano in dB. Per misurare l'effetto si confronta la potenza (spettro del segnale senza il primo mezzo secondo) dentro $\pm20$ Hz dalle armoniche e nel resto tra $100$ e $8000$ Hz (a più di $60$ Hz dalle armoniche, dove ci sono voce e folla):

(Ingresso: valore efficace $0{,}0611$, picco $0{,}373$.) Con $r=0{,}95$, il risultato è peggiore dell'ingresso: le armoniche scendono di $29$ dB, ma il resto sale di $10$ dB e il segnale in uscita ha valore efficace $1{,}75$ volte quello d'ingresso, con picco $2{,}03$, cioè oltre il fondo scala (salvando il file a $16$ bit andrebbe in saturazione). Con $r=0{,}99$ il segnale diventa $7{,}3$ dB più debole, le armoniche scendono di $18$ dB e il resto perde solo $2{,}2$ dB (il guadagno cade a $-2{,}75$ dB a $600$ Hz e $-4{,}4$ dB a $1100$ Hz, dove la voce ha energia, perché la tacca è larga $158$ Hz e le $12$ tacche sono distanti circa $230$ Hz). Con $r=0{,}998$ il segnale utile è praticamente intatto ($-0{,}14$ dB) ma le armoniche calano di soli $6{,}5$ dB: a $r=0{,}998$ la tacca è di $31$ Hz, più stretta della oscillazione della nota. Scegliere $r$ vuol dire scegliere tra voce intatta e disturbo ridotto; dalla tabella un buon compromesso è $r\approx0{,}99$ (e $r=0{,}995$ dà $-12{,}5$ dB sulle armoniche perdendo $0{,}77$ dB sul resto).

Un'avvertenza sulla normalizzazione. La scelta del corso (guadagno $1$ in continua) non è l'unica: se si normalizza a $1$ in $\hat\omega=\pi$, cioè $b_0=\frac{1+2r\cos\hat\omega_0+r^2}{2+2\cos\hat\omega_0}$, il guadagno massimo è $1$ e la cascata con $r=0{,}95$ riduce le armoniche di $52{,}7$ dB, ma il resto del segnale cala di $13{,}4$ dB (il guadagno a $600$ Hz è $-36{,}9$ dB): una tacca da $700$ Hz, ripetuta ogni $230$ Hz, equivale a un filtro che attenua di oltre $13$ dB tutto fino a $3$ kHz, voce compresa. Per $r=0{,}99$ le cifre sono $-19{,}4$ dB sulle armoniche e $-3{,}6$ dB sul resto.

Ritardo di gruppo $\tau(\hat\omega)=-\frac{d\angle H}{d\hat\omega}$ (opzionale). Il notch FIR con soli zeri sulla circonferenza ($b=[1,-2\cos\hat\omega_0,1]$, fase lineare) ha ritardo costante di $1$ campione a tutte le frequenze (lontano da $\hat\omega_0$). Il notch IIR no: con $r=0{,}99$ il ritardo di gruppo vale $0{,}13$ campioni a $3000$ Hz e $1{,}38$ a $1000$ Hz, ma cresce molto vicino alla tacca, fino a $104{,}3$ campioni ($2{,}2$ ms) a $235$ Hz ($29{,}3$ a $100$ Hz). Con $r=0{,}95$ la tacca è così larga che il massimo, $28{,}7$ campioni, si sposta a $77$ Hz, e anche a $1000$ Hz il ritardo è $5{,}8$ campioni. Le componenti vicine alla frequenza eliminata vengono ritardate molto più delle altre. È il costo dell'IIR rispetto ai filtri a fase lineare (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 →); nel caso di una tacca stretta il ritardo è grande ma riguarda solo una banda molto stretta.