Esercizio - Convoluzione a blocchi con la DFT

Testo (lezioni 15-17, corso Multimedia Signal Processing, UniPD; esempi svolti a lezione). Si considerano $x[n]=\{1,2,3,4\}$, $h[n]=\{1,1,1,1\}$ e poi $x[n]=\{1,2,3,4,5,6\}$, $h[n]=\{1,1,1\}$.

1. Calcolare la convoluzione lineare $y_{\rm lin}=x*h$ e la convoluzione circolare $y_{\rm circ}$ a $N=4$ di $x=\{1,2,3,4\}$ con $h=\{1,1,1,1\}$; spiegare la differenza e verificare che con $N\ge D_x+M-1$ (zero-padding) le due coincidono.
2. Per $x=\{1,\dots,6\}$ e $h=\{1,1,1\}$, calcolare $x*h$ con overlap-add (blocchi di $L=4$, DFT a $N=6$).
3. Lo stesso con overlap-save ($N=6$, $L=N-M+1=4$).
4. Scrivere in Python i due algoritmi e verificarli su un segnale casuale di $1000$ campioni con un FIR di $33$ coefficienti.

Teoria usata: Convoluzione a blocchi - overlap-add e overlap-saveLa convoluzione lineare y = x * h (h lungo M) si può calcolare con la FFT scegliendo N >= Dx + M - 1 (potenza di 2): costo per campione circa 3 k log2 N contro k M del calcolo diretto, quindi la FFT conviene per M abbastanza grande (con Dx circa M, già da M = 16). Per segnali lunghi o infiniti si divide x in blocchi. Overlap-add: blocchi disgiunti di lunghezza L, ognuno convoluto con h (lunghezza L+M-1 con FFT di N >= L+M-1) e le code di M-1 campioni si sommano. Overlap-save: blocchi di N campioni sovrapposti di M-1, convoluzione circolare, si scartano i primi M-1 campioni di ogni uscita (corrotti dall'aliasing). Stesso costo.Convoluzione a blocchi - overlap-add e overlap-save →, 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 →, 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 →.

Punto 1: lineare e circolare

Lineare: $y_{\rm lin}[n]=\sum_\ell x[\ell]h[n-\ell]$ ha lunghezza $D_x+M-1=4+4-1=7$: $$y_{\rm lin}=\{1,\ 3,\ 6,\ 10,\ 9,\ 7,\ 4\}.$$ Circolare a $N=4$: $y_{\rm circ}[n]=\sum_{\ell=0}^{3}x[\ell]h[(n-\ell)_4]$. Poiché $h$ è costante e vale $1$ su tutti i $4$ posti, ogni uscita è la somma di tutti i campioni, $1+2+3+4=10$: $$y_{\rm circ}=\{10,\ 10,\ 10,\ 10\}.$$ Il confronto: la circolare ha solo $4$ campioni, perché la coda della lineare ($9,7,4$ in $n=4,5,6$) si ripiega sui primi posti: $y_{\rm circ}[0]=y_{\rm lin}[0]+y_{\rm lin}[4]=1+9$, $y_{\rm circ}[1]=3+7$, $y_{\rm circ}[2]=6+4$, $y_{\rm circ}[3]=10$: tutte le somme valgono $10$. Con zero-padding a $N=8\ge7$: $x_p=\{1,2,3,4,0,0,0,0\}$, $h_p=\{1,1,1,1,0,0,0,0\}$ la circolare è $\{1,3,6,10,9,7,4,0\}$, uguale alla lineare (più uno zero). Verificato con np.fft.

Punto 2: overlap-add

Blocchi non sovrapposti di $L=4$: $x_0=\{1,2,3,4\}$, $x_1=\{5,6\}$ (il secondo è più corto: si completa con zeri). Con $N=6\ge L+M-1=4+3-1$: si completano $x_0$, $x_1$ e $h$ a $6$ campioni, si calcola $y_i=\text{IDFT}\big(\text{DFT}(x_i)\cdot\text{DFT}(h)\big)$.

• $y_0=\{1,3,6,9,7,4\}$ (indici $0..5$);
• $y_1=\{5,11,11,6,0,0\}$, da traslare di $L=4$: occupa gli indici $4..9$.

Somma:

L'uscita lineare vera ha $6+3-1=8$ campioni, $\{1,3,6,9,12,15,11,6\}$: gli ultimi due zeri sono padding. La sovrapposizione ($M-1=2$ campioni, in $n=4$ e $5$) è dove si sommano le due code: $7+5=12$ e $4+11=15$.

Punto 3: overlap-save

Blocchi di $N=6$ campioni, sovrapposti di $M-1=2$; il primo parte da due zeri.

• $b_0=\{0,0,1,2,3,4\}$: circolare con $h$: $\{7,4,1,3,6,9\}$. Si scartano i primi $M-1=2$ campioni ($7,4$) e si tengono $\{1,3,6,9\}$.
• $b_1=\{3,4,5,6,0,0\}$ (ultimi due campioni del blocco precedente, poi $5,6$ e padding): circolare $\{3,7,12,15,11,6\}$. Si scartano $3,7$ e si tengono $\{12,15,11,6\}$.

Uscita: $\{1,3,6,9,12,15,11,6\}$, uguale a $x*h$. I campioni scartati sono quelli in cui la convoluzione circolare ha ripiegato la coda: per $b_0$ la convoluzione lineare è $\{0,0,1,3,6,9,7,4\}$ e la coda $7,4$ si somma ai due zeri iniziali.

Punto 4: il codice

import numpy as np

def overlap_add(x, h, L):
    M = len(h)
    N = 1 << int(np.ceil(np.log2(L + M - 1)))          # potenza di 2 >= L+M-1
    H = np.fft.fft(h, N)                               # DFT di h, una volta sola
    y = np.zeros(len(x) + M - 1)
    for s in range(0, len(x), L):
        blk = x[s:s + L]
        yk = np.real(np.fft.ifft(np.fft.fft(blk, N) * H))[:len(blk) + M - 1]
        y[s:s + len(yk)] += yk                         # somma delle code sovrapposte
    return y

def overlap_save(x, h, N):
    M = len(h)
    L = N - M + 1                                       # campioni validi per blocco
    H = np.fft.fft(h, N)
    xp = np.concatenate([np.zeros(M - 1), x, np.zeros(N)])   # M-1 zeri iniziali
    out = []
    for s in range(0, len(x) + M - 1, L):
        blk = xp[s:s + N]
        if len(blk) < N:
            blk = np.concatenate([blk, np.zeros(N - len(blk))])
        yk = np.real(np.fft.ifft(np.fft.fft(blk) * H))
        out.append(yk[M - 1:])                          # si scartano i primi M-1 campioni
    return np.concatenate(out)[:len(x) + M - 1]

rng = np.random.default_rng(3)
x = rng.normal(size=1000); h = rng.normal(size=33)
print(np.allclose(overlap_add(x, h, 100), np.convolve(x, h)))     # True
print(np.allclose(overlap_save(x, h, 128), np.convolve(x, h)))    # True

Entrambi i controlli stampano True (eseguito). Con $x=\{1,\dots,6\}$ e $h=\{1,1,1\}$ le funzioni danno $\{1,3,6,9,12,15,11,6\}$.

Osservazione sul costo. Per ogni blocco servono una FFT diretta, un prodotto di $N$ numeri e una FFT inversa; $H=\text{DFT}(h)$ è calcolata una sola volta. Per $M=33$ e $N=128$ ogni blocco produce $L=96$ campioni validi con circa $2\cdot\frac N2\log_2N+N=1024$ moltiplicazioni complesse, cioè circa $10.7$ per campione (contro $33$ moltiplicazioni reali per campione del calcolo diretto, da confrontare tenendo conto che una moltiplicazione complessa ne vale circa quattro reali).