Esercizio - Laboratorio 1 - FIR e immagini

Testo (laboratorio 1 del corso Multimedia Signal Processing, UniPD; tema "filtri, loro applicazioni e fondamenti di elaborazione delle immagini"; codice del corso in Python con numpy, scipy e matplotlib, anche in MATLAB). Due esperienze: (1) rilevamento dei bordi di un'immagine con la varianza locale calcolata su una finestra scorrevole, con soglia pari alla media della varianza (o scelta dall'utente, per esempio $400$); la domanda del testo: a che cosa assomiglia la procedura per la varianza locale? (2) eliminazione di un'eco da un segnale audio: stima del ritardo $N$ dall'autocorrelazione e ricostruzione con $x[n]=y[n]-a\,x[n-N]$ per $a=0{,}75$; la domanda del testo: è un filtro FIR?

Dati. I file stanno nella cartella del laboratorio 1 su Moodle: eleph2.jpg (una foto di elefanti, $1600\times800$ pixel; esercizio 1) e y_echo.mat (esercizio 2, campionato a $F_s=44100$ Hz). I risultati qui sotto sono ottenuti eseguendo il codice proprio su questi file.

Teoria usata: 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 → (filtraggio 2D, convoluzione, media mobile), Sistemi a tempo discreto e filtri FIRUn sistema a tempo discreto trasforma una sequenza x[n] in una sequenza y[n]. Il filtro FIR causale di ordine M calcola y[n] = Σ b_k x[n-k] (k = 0..M): è una media mobile pesata di L = M+1 campioni, la sua risposta impulsiva h[n] coincide con i coefficienti b_k e l'uscita ha supporto lungo N+M se l'ingresso è lungo N. La media mobile è un passa-basso che ritarda di M/2 campioni; la versione centrata non è causale. Gli schemi a blocchi usano solo moltiplicatori, sommatori e ritardi unitari, senza anelli (feed-forward).Sistemi a tempo discreto e filtri FIR → (media mobile, ritardo), 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à → (FIR e IIR, poli e stabilità), Segnali e sistemi in Python - campioni, convoluzione e filtriAl calcolatore un segnale continuo è un vettore di campioni con un asse dei tempi. Con numpy si definiscono i segnali, si approssimano area ed energia con somme moltiplicate per il passo, si calcola la convoluzione continua con np.convolve(x, y)*dt (asse = somma degli istanti iniziali), si simulano filtri con equazioni alle differenze e si verificano numericamente linearità e tempo-invarianza.Segnali e sistemi in Python - campioni, convoluzione e filtri → (strumenti numpy).

Parte 1 - bordi con la varianza locale

Idea. In una zona dove l'immagine è piatta i pixel vicini hanno quasi lo stesso valore: la varianza calcolata su una piccola finestra è circa zero. Sui bordi la finestra contiene pixel scuri e chiari insieme, e la varianza è grande. Quindi la varianza locale è una mappa dei bordi, che si trasforma in immagine binaria con una soglia. Per ogni pixel $(i,j)$, con le medie fatte su una finestra (per esempio $3\times3$) centrata nel pixel: $$\mathrm{var}[i,j]=\underbrace{\overline{x^2}[i,j]}_{\text{media locale di }x^2}-\Big(\underbrace{\overline{x}[i,j]}_{\text{media locale di }x}\Big)^2.$$ Si usa la formula $\mathrm{Var}=E[x^2]-E[x]^2$, così bastano due medie locali.

A che cosa assomiglia? La media locale su una finestra $3\times3$ è un filtro FIR bidimensionale (media mobile 2D): è la convoluzione dell'immagine con un nucleo $3\times3$ con tutti i coefficienti uguali a $\frac19$, la versione 2D della media mobile di Sistemi a tempo discreto e filtri FIRUn sistema a tempo discreto trasforma una sequenza x[n] in una sequenza y[n]. Il filtro FIR causale di ordine M calcola y[n] = Σ b_k x[n-k] (k = 0..M): è una media mobile pesata di L = M+1 campioni, la sua risposta impulsiva h[n] coincide con i coefficienti b_k e l'uscita ha supporto lungo N+M se l'ingresso è lungo N. La media mobile è un passa-basso che ritarda di M/2 campioni; la versione centrata non è causale. Gli schemi a blocchi usano solo moltiplicatori, sommatori e ritardi unitari, senza anelli (feed-forward).Sistemi a tempo discreto e filtri FIR →. La varianza locale combina due filtraggi lineari (le due medie) con un'operazione non lineare (il quadrato prima della media, e quello della media): complessivamente il sistema non è LTI, ma i suoi mattoni lo sono.

Passi del codice. Si carica l'immagine, la si converte in scala di grigi e in float (per non avere overflow su interi a 8 bit), si calcolano le due medie locali, la varianza, la soglia e si disegna il risultato. Il codice del corso usa due cicli annidati con zero padding (qui con PIL al posto di OpenCV: la conversione in grigi, $0{,}299R+0{,}587G+0{,}114B$, dà gli stessi valori):

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

img_gray = np.asarray(Image.open("eleph2.jpg").convert("L")).astype(np.float64)   # 800 x 1600

def local_mean(image, window_size=3):
    if window_size % 2 == 0:
        window_size += 1                      # finestra dispari: pixel centrale
    mean_matrix = np.zeros_like(image)
    mn = window_size // 2
    padded = np.pad(image, pad_width=mn, mode="constant", constant_values=0)
    rows, cols = image.shape
    for i in range(rows):
        for j in range(cols):
            mean_matrix[i, j] = np.mean(padded[i:i+window_size, j:j+window_size])
    return mean_matrix

sz = 3
local_variance = local_mean(img_gray**2, sz) - local_mean(img_gray, sz)**2
threshold = np.mean(local_variance)           # oppure, ad esempio, 400
edges = np.zeros_like(local_variance)
edges[local_variance < threshold] = 1
edges[local_variance >= threshold] = 0

Risultati sull'immagine del corso (foto in grigi, $800\times1600$ pixel, livello medio $130{,}5$):

• La soglia alla media taglia molto in alto. La varianza è molto asimmetrica: la mediana è $21{,}2$, la media $166{,}9$, il massimo $14\,318$. Quindi la soglia pari alla media vale circa $8$ volte la mediana e lascia come "bordi" solo i pixel con varianza alta ($21{,}8\%$ dell'immagine). Il massimo possibile per pixel a $8$ bit in una finestra $3\times3$ è $16\,056$ (quattro o cinque pixel a $255$ e gli altri a $0$): il massimo misurato ne è l'$89\%$, cioè nelle zone più contrastate (per esempio le zanne e i contorni degli elefanti contro il cielo) la finestra vede quasi il contrasto massimo.
• Finestra più grande, bordi più spessi. Con $5\times5$ e $7\times7$ cresce la varianza media (ogni finestra raccoglie più contrasto) e crescono i pixel sopra la soglia fissa $400$ (da $9{,}2\%$ a $19{,}0\%$): i bordi diventano linee più larghe.
• Sulla seconda immagine della cartella (test.jpg, $225\times225$, un ritratto) la varianza media è $299{,}0$ e il $20{,}5\%$ dei pixel supera la soglia: una foto con più dettagli fini ha varianza media più alta.

Osservazioni sulla soluzione del corso

1. Colori invertiti. La soluzione mette $1$ (bianco) dove la varianza è sotto soglia, cioè sulle zone piatte, e $0$ (nero) sui bordi: i bordi risultano neri su fondo bianco. Per avere i bordi bianchi si scrive edges = local_variance >= threshold.
2. Bordi finti sul perimetro dell'immagine. Lo zero padding aggiunge una cornice di zeri: le finestre che escono dall'immagine contengono pixel veri e zeri, e la varianza è alta anche senza nessun bordo. Sull'immagine del corso, con lo zero padding $4793$ dei $4796$ pixel del perimetro risultano "bordo" (quasi tutti, e sono falsi); con il riflesso (mode="reflect") i falsi bordi scendono a $1274$, che sono per lo più bordi veri vicini ai lati dell'immagine. Cambia anche la soglia: $166{,}9$ con gli zeri, $145{,}9$ con il riflesso.
3. Velocità. I due cicli annidati su un'immagine $H\times W$ fanno $HW$ medie, ciascuna con una finestra: su questa immagine ($1{,}28$ milioni di pixel) il codice del corso impiega circa $20$ secondi per le due medie. La stessa media mobile 2D è scipy.ndimage.uniform_filter, che impiega circa $0{,}03$ secondi per ciascuna. Con mode="constant" dà lo stesso risultato dei cicli (verificato con np.allclose su una porzione).

from scipy.ndimage import uniform_filter

def local_variance(img, size=3, mode="reflect"):
    m  = uniform_filter(img,    size=size, mode=mode)      # E[x]
    m2 = uniform_filter(img**2, size=size, mode=mode)      # E[x^2]
    return m2 - m**2

v = local_variance(img_gray, 3, "reflect")
edges = v >= v.mean()                                     # True = bordo
plt.imshow(edges, cmap="gray"); plt.show()

Controllo a mano su un bordo ideale. Se un gradino separa una zona a $50$ da una a $200$, per un pixel fuori dal gradino la finestra $3\times3$ contiene $3$ pixel a $200$ e $6$ a $50$: media $100$, $\overline{x^2}=15000$, varianza $15000-100^2=5000$; per un pixel dentro, $6$ a $200$ e $3$ a $50$: media $150$, $\overline{x^2}=27500$, varianza $27500-150^2=5000$. La varianza locale vale $0$ nelle zone piatte e $5000$ sui due pixel che affiancano il gradino: un bordo diventa una striscia spessa $2$ pixel, e le zone di contorno risultano strisce di questo tipo.

Parte 2 - autocorrelazione ed eco

Modello. Il segnale registrato è il segnale originale più una copia ritardata e attenuata: $$y[n]=x[n]+a\,x[n-N],$$ con $a$ l'ampiezza dell'eco e $N$ il ritardo in campioni. Si assume $x[n]=0$ per $n<0$. È un FIR con $h=\delta[n]+a\,\delta[n-N]$ (due soli coefficienti non nulli, di ordine $N$).

Correlazione. La correlazione di $y$ e $z$ è $R_{yz}[n]=y[n]*z[-n]$ (convoluzione con la versione ribaltata); è una misura di somiglianza tra le sequenze. Con $z=y$ è l'autocorrelazione $R_y[n]=y[n]*y[-n]$. Sostituendo $y=x+a\,x[n-N]$, per bilinearità $$R_y[n]=(1+a^2)\,R_x[n]+a\,R_x[n-N]+a\,R_x[n+N].$$ Se $x$ è "rumoroso" (per esempio un segnale audio, la cui autocorrelazione è concentrata attorno a $0$), $R_x$ ha un picco in $0$ e quindi $R_y$ ha picchi in $n=0$ (il più alto), $n=+N$ e $n=-N$ (da $R_x[n\mp N]$): dal grafico di $R_y$ si legge $N$ come la posizione del picco a destra dello zero.

Il file del corso. loadmat('y_echo.mat') restituisce un dizionario con molte chiavi (Exercise1, Exercise2, Fs, N, data, x, y, y_echo, ...): il segnale con eco è solo y_echo, di $166\,255$ campioni, cioè $3{,}77$ s a $44100$ Hz, con valore massimo $1{,}0$ (è normalizzato) e valore efficace $0{,}1465$. Le altre variabili sono avanzi del file di lavoro di chi l'ha preparato: in particolare la chiave N vale $8820$ ma non è il ritardo dell'eco di y_echo, che va misurato con l'autocorrelazione.

Lettura di $N$. Con $R_y$ normalizzata ($R_y[0]=1$):

• il massimo sui ritardi positivi, preso senza cautele, non è l'eco: è nel lobo centrale e cade a $1$ campione ($R_y[1]=0{,}947$), perché un segnale audio è molto correlato con se stesso a pochi campioni di distanza;
• l'eco è il picco secondario lontano dallo zero. Cercando il massimo per ritardi maggiori di $2000$ campioni si trova $N=22050$ campioni, cioè $0{,}5$ s ($22050/44100$), con altezza $0{,}4916$ (i campioni accanto valgono $0{,}467$ e $0{,}465$: il picco è netto); il picco successivo più alto vale solo $0{,}10$ (a $21399$, ancora sulle code del lobo di $R_x$). Il valore coincide con quello della soluzione del corso;
• l'altezza è vicina a $\frac a{1+a^2}=\frac{0{,}75}{1{,}5625}=0{,}48$, il valore che si avrebbe con un segnale "bianco" (senza correlazione tra campioni lontani): da $R_y[N]=a\,(R_x[0]+R_x[2N])+(1+a^2)R_x[N]$ con $R_x[N]\approx0$ si ottiene $R_y[N]/R_y[0]\approx a/(1+a^2)$. Qui $0{,}4916$ è appena sopra, perché $R_x[N]$ non è esattamente zero;
• a $2N$ non c'è nessun picco ($R_y[2N]=0{,}018$): c'è una sola eco, come nel modello.

Rimozione dell'eco. Dato che $x[n]=0$ per $n<0$, i primi $N$ campioni di $y$ coincidono con quelli di $x$ ($x[n]=y[n]$ per $n<N$). Per $n\ge N$, da $y[n]=x[n]+a\,x[n-N]$ si ricava $x[n]=y[n]-a\,x[n-N]$, dove $x[n-N]$ è già noto: si procede campione per campione.

from scipy.io import loadmat
from scipy.signal import correlate, lfilter

y = loadmat('y_echo.mat')['y_echo'].flatten()        # Fs = 44100
Fs = 44100

r = correlate(y, y, mode="full", method="fft")      # veloce: FFT
r = r / r.max()
lag = np.arange(-len(y) + 1, len(y))
lontano = lag > 2000                                 # si salta il lobo centrale
N = lag[lontano][np.argmax(r[lontano])]              # picco a destra di 0: 22050

a = 0.75
x = np.zeros_like(y); x[:N] = y[:N]                 # i primi N campioni coincidono
for n in range(N, len(y)):
    x[n] = y[n] - a * x[n - N]                      # filtraggio inverso

np.correlate in modo full è lento: solo con $40\,000$ campioni impiega circa $14$ s e il tempo cresce col quadrato della lunghezza (su $166\,255$ campioni, circa $4$ minuti, per stima); scipy.signal.correlate con FFT impiega $0{,}014$ s (Trasformata di Fourier discreta (DFT) e convoluzione circolareLa DFT di N campioni x[0..N-1] è X[k] = sum x[n] e^{-j2pi kn/N}, k = 0..N-1 (IDFT: x[n] = (1/N) sum X[k] e^{j2pi kn/N}); vede il segnale come periodico di periodo N e dà N campioni equispaziati della DTFT, cioè X(z) valutata sulle radici N-esime dell'unità. Se x ha durata M <= N non c'è aliasing temporale e lo zero-padding (N > M) infittisce solo i campioni della stessa DTFT. Proprietà: traslazione circolare, simmetria coniugata X[k] = X*[N-k], prodotto = convoluzione circolare. La convoluzione circolare coincide con la lineare solo se N >= Dx + M - 1; altrimenti la coda si ripiega sull'inizio. La FFT calcola la DFT con (N/2) log2 N moltiplicazioni invece di N^2.Trasformata di Fourier discreta (DFT) e convoluzione circolare →). Il ciclo di ricostruzione impiega $0{,}025$ s.

Risultato della rimozione. L'autocorrelazione del segnale ricostruito $x$ a $N$ scende da $0{,}4916$ a $0{,}0215$ (l'eco è praticamente sparita), il valore efficace passa da $0{,}1465$ a $0{,}1245$ (sui primi $144\,205$ campioni) e il massimo resta $1{,}0$.

• Verifica che $a=0{,}75$ sia giusto. Il segnale originale dura $L=166\,255-22\,050=144\,205$ campioni: gli ultimi $N$ campioni di $y$ contengono solo l'eco. Se $a$ è corretto, la ricostruzione deve dare zero dopo $L$. Il valore efficace di quella coda è $0{,}0102$ per $a=0{,}75$, $0{,}0215$ per $a=0{,}5$, $0{,}0778$ per $a=1$: il minimo è proprio a $a=0{,}75$ (con una scansione fine). In alternativa, $R_x[N]$ si annulla per $a\approx0{,}77$.
• File y_echo_1.mat (stessa cartella, $80\,391$ campioni, $1{,}82$ s). Qui il picco dell'autocorrelazione è a $N=8820$ ($0{,}2$ s, altezza $0{,}4463$); prima ci sono picchi ogni $112$ campioni (altezze $0{,}88;0{,}77;0{,}66;\dots$) dovuti alla periodicità del suono, non all'eco. Con $a=0{,}75$ la cancellazione è eccessiva ($R_x[N]=-0{,}23$); la coda e la correlazione indicano $a\approx0{,}57$-$0{,}61$. Il valore $a=0{,}75$ della consegna vale solo per y_echo.mat.

È un filtro FIR? No: $x[n]=y[n]-a\,x[n-N]$ usa l'uscita passata $x[n-N]$ (retroazione): è un filtro ricorsivo, quindi IIR. La sua funzione di sistema è $H(z)=\dfrac1{1+a\,z^{-N}}$, l'inverso del FIR dell'eco $1+a\,z^{-N}$. I poli sono le $N$ radici di $z^N=-a$, tutte di modulo $a^{1/N}$; con $a=0{,}75<1$ e $N=22050$ il raggio vale $0{,}999987<1$: i poli sono dentro il cerchio unitario e il filtro è stabile (per $\lvert a\rvert\ge1$ non lo sarebbe: l'eco compensata crescerebbe; con $a=1$ la coda sopra vale $0{,}078$, molto più alta). Equivalentemente il sistema inverso si realizza in forma diretta con lfilter(b=[1], a=[1, 0, ..., 0, a]) (coefficiente $a$ in posizione $N$): sui dati veri il risultato coincide con il ciclo (differenza massima $0$), ma impiega circa $2$ s perché il denominatore ha $22\,051$ coefficienti. Il raggio dei poli molto vicino a $1$ spiega perché l'eco si spegne lentamente: ogni ripetizione è attenuata di $a=0{,}75$ ogni $N$ campioni.