Trasformata di Fourier discreta (DFT) e convoluzione circolare

La trasformata di Fourier a tempo discreto (Trasformata di Fourier a tempo discretoLa trasformata di Fourier di un segnale discreto è $S(f)=\sum_nT,s(nT),e^{-i2\pi fnT}$: una funzione continua e periodica in $f$ di periodo $F_p=1/T$. L'antitrasformata è l'integrale su un periodo, $s(nT)=\int_0^{F_p}S(f)e^{i2\pi fnT}df$. Le regole sono quelle del caso continuo (traslazione, convoluzione $\leftrightarrow$ prodotto, Parseval $\sum T|s|^2=\int_0^{F_p}|S|^2$), con incremento e somma corrente al posto di derivata e integrale. Per segnali reali $S(f)=S^*(-f)$, quindi basta $[0,F_p/2]$. Esempi: $\delta\to1$, $1\to\delta_{F_p}$, rect $\to$ sinc periodico, $a^n\mathbf 1_0\to\frac T{1-ae^{-i2\pi fT}}$.Trasformata di Fourier a tempo discreto →) di un segnale è una funzione continua della frequenza: un calcolatore non può calcolarla né memorizzarla. Può invece trattare un numero finito di valori: la DFT (discrete Fourier transform) trasforma $N$ numeri in $N$ numeri e dà campioni della DTFT. In questa nota si ripassano definizione e proprietà (con $T=1$, come nel corso; la trattazione con $T$ esplicito è in Trasformata di Fourier discreta (DFT) e FFTUn segnale discreto periodico di periodo $NT$ ha una trasformata discreta e periodica, la DFT: $S(kF)=\sum_{n=0}^{N-1}T,s(nT)e^{-i2\pi kn/N}$, con $F=1/(NT)$, e $s(nT)=\sum_{k=0}^{N-1}F,S(kF)e^{i2\pi kn/N}$. Dipende da soli $N$ numeri, si calcola senza approssimazioni e con la FFT costa $N\log_2N$ invece di $N^2$. I coefficienti di Fourier del segnale periodico sono $S_k=F,S(kF)$. La DFT dà anche campioni della trasformata di un segnale discreto o continuo a durata limitata (con zero-padding), ma la scalatura $T$ e l'asse delle frequenze vanno gestiti con cura.Trasformata di Fourier discreta (DFT) e FFT → e in FFT e zero-padding - TF, TFtd e asse delle pulsazioniLa FFT dei campioni di un segnale a durata finita, moltiplicata per il passo $T_c$, approssima la trasformata di Fourier: $X(\omega_k)\approx T_c,\mathtt{fft}(x,M)[k]$ con $\omega_k=\frac{2\pi k}{MT_c}$ (e fattore di fase $e^{-j\omega t_0}$ se l'asse parte da $t_0$). Lo zero-padding ($M>N$) infittisce i punti della stessa TFtd senza aggiungere informazione; la risoluzione dipende dalla durata osservata. Per un segnale reale $|X|$ è simmetrico: il picco in $k$ ha un gemello in $M-k$.FFT e zero-padding - TF, TFtd e asse delle pulsazioni →), si collega la DFT alla trasformata zeta e alla DTFT, e si studia la convoluzione circolare, che è il modo in cui la DFT "fa" la convoluzione.

Definizione

Un segnale finito $x[0],\dots,x[N-1]$ si può vedere come un periodo di un segnale periodico $x[n]=x[n+N]$, $n\in\mathbb Z$. Il periodo $T_p=NT$ determina la risoluzione in frequenza $F=1/T_p$ e la trasformata è periodica con periodo $F_s=NF$ ($F_s=1/T$ frequenza di campionamento).

Definizione (DFT e IDFT, con $T=1$). $$X[k]=\sum_{n=0}^{N-1}x[n]\,e^{-j2\pi kn/N},\qquad x[n]=\frac1N\sum_{k=0}^{N-1}X[k]\,e^{j2\pi kn/N},\qquad k,n=0,\dots,N-1.$$ Con la notazione $W_N=e^{-j2\pi/N}$ (radice $N$-esima primitiva dell'unità): $X[k]=\sum_nx[n]W_N^{kn}$ e $x[n]=\frac1N\sum_kX[k]W_N^{-kn}$. Gli estremi delle somme possono essere qualsiasi $N$ interi consecutivi, perché i termini sono periodici.

Esempio. $x=\{1,2,3,4\}$, $N=4$: $X[0]=1+2+3+4=10$; $X[1]=1+2(-j)+3(-1)+4(j)=-2+2j$; $X[2]=1-2+3-4=-2$; $X[3]=-2-2j$. Si ottiene $X=\{10,\ -2+2j,\ -2,\ -2-2j\}$ (verificato con Python, con l'IDFT che restituisce $x$).

Asse delle frequenze. L'indice $k$ corrisponde alla pulsazione $\hat\omega_k=2\pi k/N$ e alla frequenza $f_k=kF_s/N$. La risoluzione (frequency spacing) è $F=F_s/N$. Gli indici $k>N/2$ rappresentano frequenze negative: $k$ e $k-N$ coincidono per periodicità.

Esempio (lezione). $T=0.001\ \mathrm s$ ($F_s=1000\ \mathrm{Hz}$) e $N=4$: $F=1000/4=250\ \mathrm{Hz}$ e i bin $k=0,1,2,3$ stanno a $0,\ 250,\ 500,\ 750\ \mathrm{Hz}$; $k=4$ ricade a $1000\ \mathrm{Hz}$, cioè di nuovo a $0$ per periodicità.

L'ortogonalità $\sum_{k=0}^{N-1}e^{j2\pi k(n-m)/N}=N\delta_{n-m}$ (somma geometrica, Serie notevoli - geometrica, telescopica, armonicaLe serie di cui si conosce il carattere e da usare come termine di paragone: geometrica (converge a 1/(1-q) se |q|<1), telescopiche (somma b_1 - lim b_n, come Mengoli), armonica generalizzata (1/n^alpha converge se e solo se alpha>1).Serie notevoli - geometrica, telescopica, armonica →) garantisce che l'IDFT inverta la DFT. $X[k]$ è un numero complesso: il modulo dice quanta potenza c'è a quella frequenza, la fase di quanto è sfasata nel tempo.

DFT, trasformata zeta e DTFT

Sia $x[n]$ un segnale aperiodico a supporto finito $\{0,\dots,M-1\}$. La sua trasformata zeta è $X(z)=\sum_{n=0}^{M-1}x[n]z^{-n}$ (Trasformata zeta - definizione e regione di convergenzaLa trasformata zeta bilatera X(z) = Σ x[n] z^{-n} associa a una sequenza una funzione della variabile complessa z, definita nella regione di convergenza (ROC), sempre una corona circolare |z| in (R1, R2). Segnale a durata finita: ROC tutto il piano (tranne eventualmente 0 e ∞); causale: |z| > R1 (teorema di Abel); anticausale: |z| < R2; bilatero: intersezione, se non vuota. La stessa espressione algebrica con ROC diverse è la trasformata di segnali diversi: la ROC fa parte della trasformata. Sulla circonferenza unitaria, se è nella ROC, X(e^{jθ}) è la trasformata di Fourier. Le ROC non contengono poli.Trasformata zeta - definizione e regione di convergenza →) e la DTFT è $X(e^{j\hat\omega})$.

Teorema (DFT = campionamento della zeta sulla circonferenza). Per $N\ge M$, la DFT a $N$ punti di $x[n]$ è $$X[k]=X(z)\Big|_{z=e^{j2\pi k/N}}=X\big(e^{j\hat\omega}\big)\Big|_{\hat\omega=2\pi k/N},\qquad k=0,\dots,N-1,$$ cioè la zeta valutata nelle $N$ radici $N$-esime dell'unità: $N$ punti equispaziati sulla circonferenza unitaria a partire da $z=1$.

Esempio. $x=\{1,1,1,1\}$ e $N=8$: $X[k]=\frac{\sin(\pi k/2)}{\sin(\pi k/8)}e^{-j3\pi k/8}$, con moduli $\{4,\ 2.613,\ 0,\ 1.082,\ 0,\ 1.082,\ 0,\ 2.613\}$.

Il collegamento si capisce passando dal periodo: la DFT è la trasformata di Fourier del segnale periodico $\tilde x[n]=\sum_{\ell\in\mathbb Z}x[n-\ell N]$ ripetizione periodica di $x$ (periodic repetition) con periodo $N$. Infatti, con $m=n-\ell N$ e $e^{-j2\pi k(m+\ell N)/N}=e^{-j2\pi km/N}$, $$\tilde X[k]=\sum_{n=0}^{N-1}\sum_{\ell}x[n-\ell N]\,e^{-j2\pi kn/N}=\sum_{m\in\mathbb Z}x[m]\,e^{-j2\pi km/N}=X\big(e^{j2\pi k/N}\big),$$ perché al variare di $n=0..N-1$ e di $\ell\in\mathbb Z$ l'indice $m$ percorre una volta sola tutti gli interi.

Aliasing temporale. Campionare la DTFT in frequenza equivale a ripetere periodicamente il segnale nel tempo (dualità del Teorema del campionamento, interpolazione e aliasingTeorema di Shannon: un segnale a banda limitata $\omega_M$ si ricostruisce esattamente dai campioni se $T_c<\pi/\omega_M$ (frequenza di campionamento maggiore di quella di Nyquist $2f_{\max}$), con la formula di interpolazione ideale $x(t)=\sum_nx(nT_c)\operatorname{sinc}\left(\frac{t-nT_c}{T_c}\right)$. Sotto Nyquist c'è aliasing: le frequenze alte si confondono con quelle basse e l'informazione è persa.Teorema del campionamento, interpolazione e aliasing →).

• Se $M\le N$ le repliche non si sovrappongono: $\tilde x[n]=x[n]$ per $0\le n\le N-1$ e $x$ si recupera esattamente dall'IDFT; $N$ campioni della DTFT bastano a descrivere $x$.
• Se $M>N$ (o $x$ ha durata infinita) le repliche si sovrappongono (temporal aliasing): l'IDFT dà $\tilde x[n]\neq x[n]$.

Esempio (aliasing). $x=\{1,2,3,4,5,6\}$ ($M=6$) e $N=4$: $\tilde x[0..3]=\{1+5,\ 2+6,\ 3,\ 4\}=\{6,8,3,4\}$ (gli ultimi due campioni $5,6$ si ripiegano sull'inizio).

Zero-padding

Se $N>M$ si aggiungono $N-M$ zeri in coda. Non si aggiunge informazione ma si infittisce il campionamento della stessa DTFT: più punti sulla curva continua, per vedere meglio dove sono i picchi. Non si migliora la risoluzione (capacità di separare due righe vicine), che dipende dalla durata $M$ del segnale osservato (Analisi spettrale con la DFT - finestre e leakageAnalizzare lo spettro di un segnale con la DFT vuol dire osservarne solo L campioni (finestra rettangolare) e campionare in frequenza la DTFT. Una sinusoide non dà una riga ma il nucleo di Dirichlet della finestra centrato in w0: lobo principale largo 4pi/L e lobi laterali (il primo a -13.3 dB). Due toni più vicini della mezza larghezza del lobo non si separano (risoluzione, che dipende solo da L), e i lobi laterali di un tono mascherano i toni deboli (leakage). La DFT campiona la DTFT a 2pi k/N: l'ampiezza è corretta solo se il tono cade su un bin (periodi interi nella finestra), altrimenti si perde fino a 3.9 dB; lo zero-padding infittisce i campioni ma non migliora la risoluzione. Le finestre rastremate (Hann, Hamming, Blackman) abbassano i lobi laterali (-31, -42, -58 dB) a prezzo di un lobo principale più largo.Analisi spettrale con la DFT - finestre e leakage →).

Esempio. $x=\{1,1\}$ (due campioni). Con $N=4$: $X[k]=1+e^{-j\pi k/2}=\{2,\ 1-j,\ 0,\ 1+j\}$. Con $N=8$ (zero-padding): $Y[k]=1+e^{-j\pi k/4}$, otto punti della stessa curva $1+e^{-j\hat\omega}$, e $Y[0]=X[0]$, $Y[2]=X[1]$, $Y[4]=X[2]$, $Y[6]=X[3]$. In generale raddoppiando $N$ i valori di indice pari coincidono con i vecchi e quelli dispari sono nuovi.

Proprietà

Le proprietà sono quelle della trasformata di Fourier, con la differenza che gli indici sono modulo $N$: $(n)_N=n\bmod N$ è il resto della divisione per $N$ (per esempio $(-1)_4=3$, $(5)_4=1$).

Proprietà della DFT ($x[n]\leftrightarrow X[k]$, $N$ punti).

• Linearità: $\alpha x+\beta y\leftrightarrow\alpha X+\beta Y$.
• Traslazione circolare (circular shift): $x[(n-n_0)_N]\leftrightarrow e^{-j2\pi kn_0/N}X[k]$: il modulo non cambia, la fase sì. Traslare di $N$ campioni (o multipli) non cambia il segnale.
• Simmetria coniugata: se $x$ è reale, $X[k]=X^*[N-k]$ ($=X^*[(-k)_N]$): parte reale pari, parte immaginaria dispari, modulo pari, fase dispari. Allora $X[0]=\sum x[n]$ è reale e, se $N$ è pari, anche $X[N/2]=\sum x[n](-1)^n$ è reale. Bastano $\lfloor N/2\rfloor+1$ valori complessi per descrivere la DFT di un segnale reale.
• Convoluzione circolare: $X[k]H[k]\leftrightarrow y[n]=\sum_{\ell=0}^{N-1}x[\ell]\,h[(n-\ell)_N]=(x\circledast h)[n]$.
• Parseval: $\sum_n\lvert x[n]\rvert^2=\frac1N\sum_k\lvert X[k]\rvert^2$.

Esempio (traslazione). $x=\{1,2,3,4\}$, $X=\{10,-2+2j,-2,-2-2j\}$. Il segnale traslato di un campione, $x[(n-1)_4]=\{4,1,2,3\}$, ha DFT $\{10,\ 2+2j,\ 2,\ 2-2j\}$, uguale a $X[k]e^{-j2\pi k/4}$: i moduli sono invariati ($10,\ 2\sqrt2,\ 2,\ 2\sqrt2$).

Esempio (simmetria e Parseval). Per $x=\{1,2,3,4\}$: $X[3]=X^*[1]=-2-2j$ ✓ e $\sum x^2=30=\frac14(100+8+4+8)$ ✓.

Esempio (sinusoide su un bin). $x[n]=\cos(2\pi\cdot2n/8)$ con $N=8$: $X[2]=X[6]=N/2=4$ e tutti gli altri zero. In generale $\cos(2\pi Pn/N)\leftrightarrow\frac N2\{\delta[k-P]+\delta[k-(N-P)]\}$, e se $P=0$ o $P=N/2$ le due righe coincidono e danno $X=N$.

Per la dimostrazione della convoluzione circolare: partendo da $Y[k]=X[k]H[k]$ si antitrasforma $y[n]=\frac1N\sum_kX[k]H[k]e^{j2\pi kn/N}$, si sostituisce $H[k]=\sum_mh[m]e^{-j2\pi km/N}$ e si scambiano le somme: $y[n]=\sum_mh[m]\cdot\frac1N\sum_kX[k]e^{j2\pi k(n-m)/N}=\sum_mh[m]\,x[(n-m)_N]$, dove l'ultima somma è l'IDFT di $X$ calcolata in $n-m$, periodica di periodo $N$.

Convoluzione circolare e convoluzione lineare

La convoluzione circolare (circular convolution) è una convoluzione "su una circonferenza": per ogni $n$ si prende $h$ riflessa e traslata, ma gli indici fuori da $0..N-1$ rientrano dall'altra parte. Nell'esempio con $N=4$: $$y[0]=x[0]h[0]+x[1]h[(-1)_4]+x[2]h[(-2)_4]+x[3]h[(-3)_4]=x[0]h[0]+x[1]h[3]+x[2]h[2]+x[3]h[1],$$ perché $(-1)_4=3$, $(-2)_4=2$, $(-3)_4=1$. È commutativa e si calcola come IDFT del prodotto delle DFT.

Per un filtro FIR serve invece la convoluzione lineare, di lunghezza $D_y=D_x+M-1$ (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 →). Le due coincidono solo se non c'è aliasing temporale.

Teorema (convoluzione lineare tramite DFT). Se $x$ ha lunghezza $D_x$ e $h$ lunghezza $M$, e si sceglie $$N\ \ge\ D_y=D_x+M-1,$$ allora, dopo aver completato $x$ e $h$ con zeri fino a $N$ campioni, la convoluzione circolare a $N$ punti coincide con quella lineare: $y_{\rm lin}[n]=y_{\rm circ}[n]$ per $0\le n\le N-1$. Con $N<D_y$ la coda di $y_{\rm lin}$ (gli ultimi $D_y-N$ campioni) si ripiega sull'inizio.

Esempio (lezione). $x=\{1,2,3,4\}$ e $h=\{1,1,1,1\}$. Convoluzione lineare: $y_{\rm lin}=\{1,3,6,10,9,7,4\}$, lunghezza $4+4-1=7$. Convoluzione circolare con $N=4$: $y_{\rm circ}=\{10,10,10,10\}$ (lunghezza $4$): infatti $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$ e $y_{\rm circ}[3]=10$: le somme valgono tutte $10$. Per evitare l'aliasing si completa con zeri a $N=8\ge7$ ($x=\{1,2,3,4,0,0,0,0\}$, $h=\{1,1,1,1,0,0,0,0\}$) e si ottiene $y_{\rm circ}=\{1,3,6,10,9,7,4,0\}=y_{\rm lin}$ più uno zero.

Esempio (filtro corto). $x=\{1,2,3,4\}$, $h=\{1,1\}$, $N=4$: $y_{\rm circ}=\{5,3,5,7\}$, mentre $y_{\rm lin}=\{1,3,5,7,4\}$: l'ultimo campione $4$ si somma al primo ($1+4=5$).

Quando conviene la DFT. Con la FFT il costo della convoluzione cala molto per filtri lunghi; il confronto dei costi e i metodi per segnali lunghi sono nella nota 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 →.

La FFT

Calcolare la DFT con la definizione richiede, per ogni $k$, $N$ moltiplicazioni complesse e $N-1$ somme: in tutto $\mathcal O(N^2)$ operazioni. La FFT (fast Fourier transform) è una famiglia di algoritmi con costo $\mathcal O(N\log_2N)$, basati su "divide et impera" (divide and conquer): si spezza il problema in due di dimensione metà, si risolvono e si ricombinano le soluzioni con costo lineare; iterando fino a dimensione $1$ ci sono $\log_2N$ livelli.

Decimazione in tempo (decimation in time, $N=2^m$). Si dividono i campioni in quelli di indice pari $x[2r]$ e dispari $x[2r+1]$. Con $W_N^2=W_{N/2}$: $$X[k]=\underbrace{\sum_{r}x[2r]W_{N/2}^{rk}}_{E[k]}+W_N^k\underbrace{\sum_rx[2r+1]W_{N/2}^{rk}}_{O[k]},\qquad k=0,\dots,N-1,$$ dove $E$ e $O$ sono due DFT a $N/2$ punti, periodiche di periodo $N/2$. Poiché $W_N^{k+N/2}=-W_N^k$, bastano $E[k]$ e $O[k]$ per $k=0,\dots,N/2-1$: $$X[k]=E[k]+W_N^kO[k],\qquad X\!\left[k+\tfrac N2\right]=E[k]-W_N^kO[k].$$ Questa coppia di operazioni (una moltiplicazione, una somma e una differenza) è la farfalla (butterfly). La ricombinazione costa $N/2$ moltiplicazioni per livello, in tutto $\frac N2\log_2N$ moltiplicazioni invece di $N^2$.

Esempio. $x=\{1,2,3,4\}$: pari $\{1,3\}$ danno $E=\{4,-2\}$, dispari $\{2,4\}$ danno $O=\{6,-2\}$, con $W_4^0=1$ e $W_4^1=-j$. Allora $X[0]=4+6=10$, $X[1]=-2+(-j)(-2)=-2+2j$, $X[2]=4-6=-2$, $X[3]=-2-2j$: coincide con la DFT calcolata sopra. Per $N=1024$: $\frac N2\log_2N=5120$ moltiplicazioni contro $1\,048\,576$, un fattore $\approx205$.

Domande d'esame

1. Definire la DFT e spiegarne il legame con la DTFT e con la trasformata zeta. Che cosa succede se il segnale è più lungo di $N$? Traccia: $X[k]=X(e^{j2\pi k/N})$, radici $N$-esime dell'unità; ripetizione periodica $\tilde x[n]=\sum_\ell x[n-\ell N]$ e dimostrazione del cambio di indice; $M\le N$ nessun aliasing, $M>N$ aliasing temporale (esempio numerico); zero-padding infittisce ma non migliora la risoluzione.
2. Enunciare la proprietà della convoluzione circolare e le condizioni sotto cui coincide con la lineare. Traccia: $Y=XH\leftrightarrow\sum x[\ell]h[(n-\ell)_N]$ con derivazione; lunghezza $D_x+M-1$; $N\ge D_x+M-1$ con zero-padding; esempio $\{1,2,3,4\}$ e $\{1,1,1,1\}$ ($\{10,10,10,10\}$ contro $\{1,3,6,10,9,7,4\}$).
3. Ricavare la FFT a decimazione in tempo e il suo costo. Traccia: separazione pari/dispari, $W_N^2=W_{N/2}$, $X[k]=E[k]+W_N^kO[k]$ e $X[k+N/2]=E[k]-W_N^kO[k]$, $\frac N2\log_2N$ moltiplicazioni, confronto con $N^2$.