Esercizio - Laboratorio 2, filtri notch e audio

Testo (laboratorio 2 del corso Multimedia Signal Processing, UniPD; lezione 14). Un file audio contiene un messaggio disturbato da toni sinusoidali aggiunti. Si vogliono togliere i toni, in due passi: (a) determinare dallo spettrogramma le frequenze dei disturbi; (b) filtrarle con filtri FIR notch. Più precisamente:

1. leggere il file, plottare lo spettrogramma normalizzato (frequenza in unità di $\pi$ rad/campione) e leggere le frequenze dei toni;
2. ricavare la funzione di sistema di un filtro FIR a tre coefficienti $h=\{1,A,1\}$ che annulla un coseno e normalizzarla a guadagno $1$ in continua ($\hat\omega=0$);
3. progettare un filtro per ogni tono, metterli in cascata (convoluzione), tracciare modulo e fase della risposta totale;
4. filtrare l'audio, plottare lo spettrogramma dell'uscita e riascoltare: l'operazione dà risultati soddisfacenti?

Dati. Il file SunshineSquare.wav sta nella cartella del laboratorio 2 su Moodle. È un file mono a $16$ bit, $F_s=11025$ Hz, $119\,071$ campioni ($10{,}8$ s). I risultati sotto sono ottenuti eseguendo il codice su questo file.

Teoria usata: 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 →, Sistemi LTI e convoluzione discretaUn sistema lineare e tempo-invariante (LTI) è completamente descritto dalla sua risposta impulsiva h[n]: y[n] = Σ x[k] h[n-k] = x[n] * h[n] (somma di convoluzione). Si ricava scomponendo x in impulsi traslati e usando linearità e invarianza. La convoluzione è commutativa, associativa e distributiva; una cascata di LTI equivale a un solo filtro con h = h1 * h2 e l'ordine non conta. Se x ha lunghezza N e h lunghezza L, y ha lunghezza N+L-1. Si calcola col metodo della finestra scorrevole o a tabella. I FIR sono LTI; altri sistemi (n·x[n], x[-n], x²) non lo sono.Sistemi LTI e convoluzione discreta →, Trasformata di Fourier a tempo breve e spettrogrammaPer un segnale lungo il cui contenuto in frequenza cambia nel tempo si calcolano tante DFT brevi: X[k, l] = sum_{n=0}^{L-1} w[n] x[l+n] e^{-j2pi kn/N}, con finestra di analisi w di lunghezza L, istante l = m R (passo R, di solito L/2) e DFT a N >= L punti. Il modulo (o i dB) di X[k, l] come immagine tempo-frequenza è lo spettrogramma. La finestra impone un compromesso: lunga = ottima risoluzione in frequenza (Delta w circa 8pi/L per Hann) ma transizioni temporali sfocate su circa L/2 campioni; corta = buona localizzazione nel tempo ma righe vicine fuse. Si provano più lunghezze.Trasformata di Fourier a tempo breve e spettrogramma →.

Passo 1: le frequenze dei toni

Con le impostazioni del corso (finestra di Hamming di $256$ campioni, $50\%$ di sovrapposizione, nfft=256) lo spettrogramma ha una riga ogni $F_s/256=43{,}07$ Hz, cioè ogni $0{,}0078\pi$. Mediato nel tempo mostra tre righe orizzontali molto marcate (da $40$ a $46$ dB sopra il resto) nei canali di $1593{,}5$, $3143{,}8$ e $4737{,}3$ Hz, cioè a $0{,}289\pi$, $0{,}570\pi$, $0{,}859\pi$. Il risultato è limitato dalla risoluzione del grafico: la FFT dell'intero file ($119\,071$ campioni) dà $$f_1=1575{,}0\ \mathrm{Hz},\quad f_2=3150{,}0\ \mathrm{Hz},\quad f_3=4725{,}0\ \mathrm{Hz},$$ cioè $\hat\omega_1=0{,}2857\pi$, $\hat\omega_2=0{,}5714\pi$, $\hat\omega_3=0{,}8571\pi$: sono esattamente $F_s/7$, $2F_s/7$, $3F_s/7$, cioè $\hat\omega=\frac{2\pi}{7},\frac{4\pi}{7},\frac{6\pi}{7}$. Le frequenze della soluzione del corso ($0{,}2857\pi$, $0{,}5709\pi$, $0{,}8573\pi$) sono quelle lette dal grafico: la prima è esatta, la seconda e la terza sono a $0{,}0005\pi$ e $0{,}0002\pi$ da quelle vere. Il fatto che siano multipli di $\frac{2\pi}7$ non è un caso e si userà nel passo 3.

Com'è fatto il file:

• i toni sono presenti solo nell'ultimo tratto, da circa $7{,}5$ s alla fine ($10{,}8$ s); nei primi $7{,}5$ s c'è solo il messaggio. Nel tratto $8$-$10{,}5$ s l'ampiezza dei tre toni (stimata con la DFT del tratto) è $1{,}163$, $0{,}287$, $0{,}415$ (scala in cui il fondo scala è $1$): il primo tono da solo supera il fondo scala, ma i tre toni sono sfasati e la loro somma resta sotto (picco del file $0{,}99997$);
• il messaggio ha valore efficace $0{,}0144$ (stimato nei primi $6$ s, senza continua) e il tratto con i toni $0{,}897$: i toni sono circa $36$ dB sopra il messaggio, che per questo non si sente;
• il file ha una componente continua di $0{,}036$ (maggiore del valore efficace del messaggio);
• nel messaggio (primi $6$ s, senza continua) la potenza è concentrata nelle basse frequenze: il $50\%$ sotto $541$ Hz ($0{,}098\pi$) e il $90\%$ sotto $1547$ Hz ($0{,}28\pi$).

Passo 2: il filtro per un tono

Un coseno è somma di due esponenziali a $\pm\hat\omega_0$ (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 →); servono quindi due zeri in $z_{1,2}=e^{\pm j\hat\omega_0}$: $$H(z)=(1-e^{j\hat\omega_0}z^{-1})(1-e^{-j\hat\omega_0}z^{-1})=1-2\cos(\hat\omega_0)\,z^{-1}+z^{-2}.$$ Si sono usate $e^{j\hat\omega_0}+e^{-j\hat\omega_0}=2\cos\hat\omega_0$ e $e^{j\hat\omega_0}e^{-j\hat\omega_0}=1$. I coefficienti sono $b_0=1$, $b_1=A=-2\cos\hat\omega_0$, $b_2=1$, cioè hh = [1, AA, 1]. Il guadagno in continua è $H(1)=2-2\cos\hat\omega_0$, quindi la versione normalizzata è $h=\{1,A,1\}/(2-2\cos\hat\omega_0)$.

Passo 3: la cascata

La risposta impulsiva della cascata è la convoluzione dei tre filtri: ha $3+3-1=5$ coefficienti dopo due filtri e $5+3-1=7$ dopo tre. $$h_{\rm out}=\{0.1430,\ 0.1426,\ 0.1427,\ 0.1432,\ 0.1427,\ 0.1426,\ 0.1430\},\qquad\sum h_{\rm out}[n]=1.$$ Risultato notevole: i sette coefficienti sono praticamente uguali a $1/7=0.1429$ (scarto massimo $0.0004$). I tre toni sono infatti, come si è visto nel passo 1, i multipli $\frac27\pi$, $\frac47\pi$, $\frac67\pi$, cioè gli zeri $e^{j2\pi k/7}$ di una media mobile a 7 punti. Il filtro è simmetrico (tipo I, Filtri FIR a fase lineare - tipi e zeriI FIR a fase lineare di ordine N si dividono in quattro tipi: I (N pari, h simmetrica), II (N dispari, simmetrica), III (N pari, antisimmetrica), IV (N dispari, antisimmetrica). Si scrive H = e^{-jwN/2} e^{jb} A(w), con A somma di coseni (tipi I e II) o di seni (III e IV). Dalla relazione H(z) = ±z^{-N} H(1/z) segue che gli zeri vengono in gruppi (z0, z0*, 1/z0, 1/z0*) e che ci sono zeri forzati: tipo II in z = -1, tipo III in z = 1 e z = -1, tipo IV in z = 1, tipo I nessuno. Quindi II non fa passa-alto, III non fa passa-basso né passa-alto, IV non fa passa-basso. Ogni tipo II, III, IV è un tipo I moltiplicato per (1+z^-1), (1-z^-2), (1-z^-1).Filtri FIR a fase lineare - tipi e zeri →), con fase lineare e ritardo $3$ campioni; il modulo vale $0$ nei tre zeri e $1$ in continua.

Poiché le frequenze del corso sono arrotondate, gli zeri della cascata cadono a $0{,}2857\pi$, $0{,}5709\pi$, $0{,}8573\pi$ (il modulo vale lì $10^{-16}$) e non esattamente sulle frequenze vere dei toni. Calcolato in $\frac{2\pi}7,\frac{4\pi}7,\frac{6\pi}7$ il modulo vale $5{,}2\cdot10^{-5}$ ($-85{,}7$ dB), $1{,}06\cdot10^{-3}$ ($-59{,}5$ dB), $2{,}5\cdot10^{-4}$ ($-71{,}9$ dB): il tono a frequenza intermedia è quello attenuato meno, perché è quello letto peggio dal grafico.

La curva continua è il modulo della cascata del corso: tre zeri, e fra uno zero e l'altro i lobi laterali arrivano a $0{,}23$ (cioè è una media mobile a $7$ punti, il cui modulo è $\left|\frac{\sin(7\hat\omega/2)}{7\sin(\hat\omega/2)}\right|$). Il guadagno a $\hat\omega=\pi$ è $0{,}143$; il guadagno è $0{,}81$ a $0{,}1\pi$ e scende a $0{,}37$ a $0{,}2\pi$. La curva tratteggiata è la cascata di tre notch IIR del laboratorio 5 (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 →) sulle stesse frequenze: resta vicina a $1$ quasi dappertutto.

Passo 4: il codice

import numpy as np
from scipy.io import wavfile
from scipy.signal import spectrogram, lfilter, convolve, freqz

fs, xx = wavfile.read('SunshineSquare.wav')      # int16, fs = 11025
if xx.ndim > 1:
    xx = xx[:, 0]                                  # primo canale (qui il file e' mono)
xx = xx / 32768.0                                  # float in [-1, 1) (come soundfile.read)

# spettrogramma normalizzato (Hamming 256, 50 % di sovrapposizione, come nel corso)
f, t, Sxx = spectrogram(xx, fs, window=np.hamming(256), nfft=256, noverlap=128)
w_norm = f / (fs / 2)                              # in unita' di pi rad/campione

# frequenze lette dal corso (unita' di pi) e filtri notch normalizzati
w_toni = np.array([0.2857, 0.5709, 0.8573]) * np.pi

def notch(w0):
    A = -2 * np.cos(w0)
    return np.array([1, A, 1]) / (2 - 2 * np.cos(w0))     # guadagno 1 in continua

hh = [notch(w) for w in w_toni]
h_out = convolve(convolve(hh[0], hh[1]), hh[2])           # cascata = convoluzione

ww = np.linspace(-np.pi, np.pi, 201)
_, HH = freqz(h_out, 1, worN=ww)                          # modulo: np.abs(HH), fase: np.unwrap(np.angle(HH))

yy = lfilter(h_out, 1, xx)                                # filtraggio

def ampiezza_tono(s, w0):                                 # ampiezza stimata con la DFT del tratto
    k = np.arange(s.size)
    return 2 * abs(np.sum(s * np.exp(-1j * w0 * k))) / s.size
tratto = slice(int(8 * fs), int(10.5 * fs))               # tratto in cui i toni sono presenti
for w in 2 * np.pi * np.array([1, 2, 3]) / 7:
    print(w / np.pi, ampiezza_tono(xx[tratto], w), ampiezza_tono(yy[tratto], w))

Uscita ottenuta eseguendo il codice sul file del corso.

• Ampiezze dei tre toni nel tratto $8$-$10{,}5$ s: da $1{,}163$, $0{,}287$, $0{,}415$ nell'ingresso a $8{,}2\cdot10^{-5}$, $3{,}0\cdot10^{-4}$, $9{,}8\cdot10^{-5}$ nell'uscita, cioè attenuazioni di $83{,}0$, $59{,}6$ e $72{,}6$ dB (coerenti con il modulo della cascata nelle tre frequenze vere, $85{,}7$, $59{,}5$, $71{,}9$ dB).
• Valore efficace nel tratto $8$-$10{,}5$ s: da $0{,}897$ a $0{,}0100$, cioè $-39$ dB. Il valore efficace del messaggio nei primi $6$ s è $0{,}0144$ in ingresso e $0{,}0109$ in uscita: l'uscita nel tratto con i toni ($0{,}0100$) ha lo stesso ordine di grandezza dell'uscita senza toni, quindi i toni sono davvero spariti e resta il messaggio.
• Sul file intero il valore efficace passa da $0{,}476$ a $0{,}069$ ($-16{,}8$ dB): quasi tutta l'energia del file ($97{,}7\%$) era nei tre toni.

Passo 5: il filtraggio è soddisfacente?

• Sì per i toni: gli zeri sono praticamente esatti, quindi dopo un transitorio di $6$ campioni (ordine $6$, meno di $1$ ms a $11025$ Hz) i toni spariscono, con attenuazioni da $60$ a $83$ dB.
• Ma si perde parte del messaggio: la cascata è una media mobile a $7$ punti, cioè un passabasso con frequenza di taglio a $-3$ dB di $0{,}128\pi$ ($705$ Hz). Sul messaggio dei primi $6$ s la potenza in uscita rispetto all'ingresso, per fasce, è: $-0{,}2$ dB sotto $0{,}05\pi$, $-1{,}1$ dB tra $0{,}05\pi$ e $0{,}1\pi$, $-3{,}3$ dB tra $0{,}1\pi$ e $0{,}2\pi$, poi $-13{,}5$ dB ($0{,}2\pi$-$0{,}3\pi$), $-15{,}4$ dB ($0{,}3\pi$-$0{,}5\pi$) e $-19{,}0$ dB sopra $0{,}5\pi$. In totale il valore efficace del messaggio cala di $2{,}5$ dB. Poiché il $90\%$ della potenza del messaggio è sotto $1547$ Hz il suono resta riconoscibile, ma perde le componenti sopra circa $1$ kHz: risulta più cupo.
• Il guadagno non è $1$ ovunque: il singolo notch ha guadagno $\cot^2(\hat\omega_0/2)$ in $\hat\omega=\pi$ (per esempio $5.83$ a $\pi/4$); nella cascata il massimo è $1$ in continua e il guadagno in $\hat\omega=\pi$ è $0.143$.
• Alternativa con notch IIR (laboratorio 5): tre notch del secondo ordine sulle frequenze esatte ($\frac{2\pi}7,\frac{4\pi}7,\frac{6\pi}7$) con $r=0{,}9$ attenuano i toni di $95{,}4$, $93{,}5$ e $94{,}3$ dB e riducono il valore efficace del messaggio di soli $0{,}39$ dB (con $r=0{,}95$: $0{,}18$ dB; con $r=0{,}99$: $0{,}03$ dB). Il costo è un transitorio più lungo e la fase non 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 →).