Convoluzione a blocchi - overlap-add e overlap-save

Filtrare un segnale $x[n]$ con un FIR di risposta impulsiva $h[n]$ lunga $M$ vuol dire calcolare $y[n]=\sum_{k=0}^{M-1}h[k]x[n-k]$ (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 →). Con la formula diretta ogni campione di uscita costa $M$ moltiplicazioni e $M-1$ somme. Con la FFT si può fare di meglio per filtri lunghi, ma la DFT fornisce una convoluzione circolare (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 →) e lavora su segnali di lunghezza finita: serve un metodo per ottenere quella lineare e per trattare segnali molto lunghi (o infiniti, come un flusso audio). Questa nota riassume il calcolo con la FFT, il confronto di costo e i due metodi a blocchi.

Convoluzione lineare con la DFT

Il procedimento, per $x$ di lunghezza $D_x$ e $h$ di lunghezza $M$ (uscita lunga $D_y=D_x+M-1$):

1. Si sceglie $N\ge D_y$; per sfruttare la FFT, $N$ potenza di due: $N=\min\{2^m:\ 2^m\ge D_y\}$.
2. Si completano $x$ e $h$ con zeri fino a $N$ campioni (zero-padding).
3. Si calcolano le DFT $X[k]$ e $H[k]$ (due FFT), si moltiplica $Y[k]=X[k]H[k]$ ($N$ moltiplicazioni) e si antitrasforma (una FFT inversa).

Con $N\ge D_y$ non c'è aliasing temporale e $y[n]=y_{\rm circ}[n]$ per $0\le n\le N-1$; con $N<D_y$ la coda si ripiega sull'inizio.

Esempio. Se $D_y=9$ serve $N=16$; se $D_y=27$, $N=32$ (come nell'appunto del corso).

Costo: quando conviene la FFT

Il costo del metodo diretto per campione è $T_D=k_DM$ ($k_D$ dipende dal processore). Il metodo FFT ha costo per campione $$T_F=\frac1N\big(3\,T_{\rm FFT}+N\,k_{\rm mult}\big)\simeq3\,k_F\log_2N,$$ dove $T_{\rm FFT}=k_FN\log_2N$ è il costo di una FFT, il fattore $3$ conta le due FFT e la FFT inversa, e il termine $N\,k_{\rm mult}$ delle $N$ moltiplicazioni si trascura. Assumendo $k_F\simeq k_D$: $$\frac{T_F}{T_D}\simeq\frac{3\log_2N}{M}.$$ La FFT conviene quando questo rapporto è minore di $1$. Due casi tipici.

Caso 1: $D_x\simeq M$. Allora $N\simeq D_x+M-1\simeq2M$ e $\frac{T_F}{T_D}\simeq\frac{3\log_2(2M)}{M}=\frac{3(1+\log_2M)}{M}$.

Per $M$ piccolo (8) è più veloce il calcolo diretto; per $M=16$ si è in pareggio; per $M$ grande la FFT è molto più veloce (verificato con Python: $3(k+1)/2^k$ per $M=2^k$).

Caso 2: $D_x\gg M$. Allora $N\simeq D_x$ e $\frac{T_F}{T_D}\simeq\frac{3\log_2D_x}{M}<1$ se $D_x<2^{M/3}$.

Quindi, se il filtro è corto, il calcolo diretto è meglio; se è lungo, la FFT è molto più veloce, anche per segnali lunghissimi. Resta però il problema pratico: una FFT enorme richiede di aspettare tutto il segnale (non va per l'elaborazione in tempo reale) e occupa molta memoria. Si risolve a blocchi.

Convoluzione a blocchi

Si suddivide $x[n]$ in blocchi, si convolve ogni blocco con $h[n]$ (con la FFT, su $N$ piccolo) e si ricombinano i risultati in modo opportuno. I due metodi, con la stessa complessità, sono overlap-add e overlap-save.

Overlap-add

Metodo (overlap-add, "sovrapponi e somma"). Sia $h$ di lunghezza $M$ e si divida $x$ in blocchi non sovrapposti di lunghezza $L$ ($L>M$ per semplicità): $$x_i[n]=x[n]\,w_R[n-iL],\qquad w_R[n]=u[n]-u[n-L],\qquad x[n]=\sum_{i=0}^{\infty}x_i[n]\ \ (x\text{ causale}).$$ Per la linearità $y[n]=\sum_ix_i[n]*h[n]=\sum_iy_i[n]$ con $y_i=x_i*h$ di lunghezza $L+M-1$ e supporto $\{iL,\dots,(i+1)L+M-2\}$.

I blocchi $y_i$ e $y_{i+1}$ hanno supporti $[iL,(i+1)L+M-2]$ e $[(i+1)L,(i+2)L+M-2]$: si sovrappongono per $M-1$ campioni, quelli di indice $(i+1)L\le n\le(i+1)L+M-2$. Quindi in ogni istante $n$ la somma ha uno o due termini non nulli: $$y[n]=\begin{cases}y_{i-1}[n]+y_i[n]&iL\le n\le iL+M-2\quad(\text{sovrapposizione con il blocco precedente})\\ y_i[n]&iL+M-2<n<(i+1)L\quad(\text{solo il blocco }i)\end{cases}$$

Per calcolare ogni $y_i=x_i*h$ con la FFT, per evitare l'aliasing la DFT deve avere lunghezza $N\ge L+M-1$. La DFT di $h$ ($H[k]$, a $N$ punti) si calcola una volta sola. Per ogni blocco: FFT di $x_i$ completato con zeri a $N$, prodotto con $H[k]$, FFT inversa; poi si traslano le uscite di $iL$ e si sommano le code sovrapposte.

Esempio (lezione). $x=\{1,2,3,4,5,6\}$, $h=\{1,1,1\}$ ($M=3$), $L=4$ e $N=6\ge L+M-1=6$. Blocchi: $x_0=\{1,2,3,4\}$ e $x_1=\{5,6\}$. Con zero-padding a $6$ campioni: $y_0=x_0*h=\{1,3,6,9,7,4\}$ (indici $0$-$5$) e $y_1=x_1*h=\{5,11,11,6\}$, che parte da $n=4$. Somma:

Il risultato coincide con la convoluzione lineare diretta $\{1,3,6,9,12,15,11,6\}$ (verificato con Python). Nell'appunto del corso compaiono due zeri finali: sono il padding della DFT a $N=6$ sul secondo blocco, non fanno parte del risultato.

Overlap-save

L'idea è opposta: invece di sommare le code, si evita di avere code. I blocchi di ingresso sono sovrapposti di $M-1$ campioni, si calcola la convoluzione circolare a $N$ punti, e si scartano i campioni che non coincidono con la convoluzione lineare.

Metodo (overlap-save, "sovrapponi e conserva"). Sia $N$ la lunghezza della DFT e $L=N-M+1$ il numero di campioni validi per blocco. Ogni blocco di ingresso ha lunghezza $N$ ed è formato da $[\,M-1$ campioni sovrapposti (gli ultimi del blocco precedente) $\mid$ $L$ campioni nuovi$\,]$. Il primo blocco ha $M-1$ zeri iniziali. Per ogni blocco: $Y_i[k]=X_i[k]H[k]$, $y_i=$ IDFT, e si scartano i primi $M-1$ campioni di $y_i$; i restanti $L$ vanno in uscita.

Perché funziona. La convoluzione circolare a $N$ punti ripiega la coda dell'uscita lineare, lunga $N+M-1$, sull'inizio: sono corrotti i primi $M-1$ campioni (esattamente quelli in cui il filtro "vede" l'ingresso oltre il bordo del blocco), mentre gli altri $N-M+1=L$ coincidono con la convoluzione lineare. Nei primi $M-1$ posti del blocco ci sono i campioni del blocco precedente, che servono proprio a rendere corretti i campioni successivi.

Esempio (lezione). Stessi dati, $x=\{1,\dots,6\}$, $h=\{1,1,1\}$ ($M=3$), $N=6$, $L=N-M+1=4$. Blocchi sovrapposti di $M-1=2$:

• $b_0=\{0,0,1,2,3,4\}$ (due zeri iniziali): convoluzione circolare $\{7,4,1,3,6,9\}$; si scartano i primi due ($7,4$) e si tengono $\{1,3,6,9\}$.
• $b_1=\{3,4,5,6,0,0\}$ (gli ultimi due campioni del blocco, poi i nuovi $5,6$ e zeri): circolare $\{3,7,12,15,11,6\}$; si scartano $3,7$ e si tengono $\{12,15,11,6\}$.

Uscita: $\{1,3,6,9\}\cup\{12,15,11,6\}=\{1,3,6,9,12,15,11,6\}$, uguale alla convoluzione lineare (verificato con Python). I campioni scartati ($7,4$ e $3,7$) sono proprio quelli corrotti dall'aliasing: per $b_0$ la convoluzione lineare è $\{0,0,1,3,6,9,7,4\}$ e la coda $7,4$ si ripiega sui primi due posti, dove c'erano gli zeri.

Confronto e scelta

I due metodi hanno la stessa complessità: per ogni blocco due FFT (blocco e inversa; $H[k]$ è calcolata una volta sola) e $N$ moltiplicazioni. Si sceglie $N$ potenza di due, di solito alcune volte $M$ (come regola pratica $N\approx4M$-$8M$): con $N$ troppo vicino a $M$ i campioni validi $L=N-M+1$ sono pochi e si butta via molto lavoro; con $N$ troppo grande cresce il fattore $\log_2N$ e la memoria. In tempo reale il filtro introduce una latenza di almeno $L$ campioni (il blocco va raccolto prima di essere elaborato) più il ritardo del filtro stesso, $(M-1)/2$ se è 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 →).

Il codice Python dei due metodi, confrontati con np.convolve, è in Esercizio - Convoluzione a blocchi con la DFT.

Domande d'esame

1. Spiegare come si calcola la convoluzione lineare di due sequenze finite con la DFT e quando conviene rispetto al calcolo diretto. Traccia: $N\ge D_x+M-1$ (potenza di 2), zero-padding, $Y=XH$; costo $T_F\simeq3k\log_2N$ per campione contro $T_D=kM$; rapporto $3\log_2N/M$; caso $D_x\simeq M$ (tabella: pareggio a $M=16$) e $D_x\gg M$ ($D_x<2^{M/3}$).
2. Descrivere il metodo overlap-add e dimostrare la formula di ricombinazione. Traccia: $x=\sum x_i$, $y=\sum y_i$, supporto $[iL,(i+1)L+M-2]$, sovrapposizione di $M-1$ campioni, tre intervalli della somma, $N\ge L+M-1$, esempio numerico $\{1..6\}*\{1,1,1\}$.
3. Descrivere overlap-save e confrontarlo con overlap-add. Traccia: blocchi di $N$ con $M-1$ campioni sovrapposti, $L=N-M+1$, convoluzione circolare, primi $M-1$ campioni scartati perché corrotti dal ripiegamento; stessa complessità; scelta di $N$; latenza.