Trasformata di Fourier discreta (DFT) e FFT

I segnali discreti periodici sono gli unici rappresentabili esattamente in un calcolatore, perché sono specificati da un numero finito di valori (quelli di un periodo). Anche le loro trasformate sono discrete e periodiche: la trasformata di Fourier discreta (DFT, discrete Fourier transform) è lo strumento con cui si studiano al calcolatore tutti i segnali, previa un'opportuna approssimazione.

Le quattro trasformate: il quadro

Lo studio di un segnale dipende da due caratteristiche: il dominio ($\mathbb R$ continuo oppure $\mathbb Z(T)$ discreto) e la periodicità (aperiodico oppure periodico di periodo $T_p$). Ogni classe ha la propria trasformata, e la regola per trovare il dominio in frequenza è la stessa: se il tempo ha "quanto" $T$ e "periodo" $T_p$ (con le convenzioni $T=0$ per il continuo e $T_p=\infty$ per l'aperiodico), la frequenza ha quanto $F=\frac1{T_p}$ e periodo $F_p=\frac1T$.

In forma unificata, tutte hanno lo stesso aspetto, $S(f)=\int_I dt\,s(t)e^{-i2\pi ft}$ e $s(t)=\int_{\hat I}df\,S(f)e^{i2\pi ft}$, dove $\int_I$ è l'integrale di Haar: integrale su $\mathbb R$, integrale su un periodo, somma $\sum T\,s(nT)$ su tutti gli interi, somma su un periodo, a seconda della classe. Le regole di calcolo (linearità, traslazione, convoluzione $\leftrightarrow$ prodotto, simmetria) e il teorema di Parseval $\int_I|s|^2=\int_{\hat I}|S|^2$ valgono in forma unica; le formule per le singole classi sono casi particolari, e alcune regole (derivazione) esistono solo per i segnali continui.

Definizione della DFT

Sia $s(nT)$ periodico di periodo $T_p=NT$ ($N$ campioni per periodo). La frequenza fondamentale è $F=\frac1{NT}$ e la trasformata vive su $\mathbb Z(F)$ con periodo $F_p=\frac1T=NF$.

Definizione (DFT e DFT inversa). $$S(kF)=\sum_{n=0}^{N-1}T\,s(nT)\,e^{-i2\pi kn/N},\qquad s(nT)=\sum_{k=0}^{N-1}F\,S(kF)\,e^{i2\pi kn/N}.$$ $S(kF)$ è periodica in $k$ di periodo $N$ (perché $e^{-i2\pi(k+N)n/N}=e^{-i2\pi kn/N}$) ed è specificata da $N$ valori, $k=0,\dots,N-1$. I coefficienti di Fourier del segnale periodico sono $S_k=F\,S(kF)=\frac1N\sum_{n=0}^{N-1}s(nT)e^{-i2\pi kn/N}$.

Esempio. Con $N=4$, $T=1$ e $s=(1,2,0,-1)$: $S(kF)=(2,\ 1-3i,\ 0,\ 1+3i)$ e i coefficienti $S_k=\frac14S(kF)=(0{,}5,\ 0{,}25-0{,}75i,\ 0,\ 0{,}25+0{,}75i)$. Il primo è il valor medio $\frac{1+2+0-1}4$.

La seconda formula discende da $\sum_{k=0}^{N-1}e^{i2\pi k(n-m)/N}=N\delta_{nm}$ (ortogonalità su $N$ punti, somma geometrica di ragione $e^{i2\pi(n-m)/N}$). Proprietà:

• Parseval. $\sum_{n=0}^{N-1}T|s(nT)|^2=\sum_{k=0}^{N-1}F|S(kF)|^2$ (per l'esempio sopra $6=\frac14(4+10+0+10)$).
• Segnale reale: $S(-kF)=S^*(kF)$, cioè $S((N-k)F)=S^*(kF)$ (basta metà dei valori).
• Traslazione ciclica di $m$ campioni: moltiplica per $e^{-i2\pi km/N}$; convoluzione ciclica $\leftrightarrow$ prodotto: $\mathrm{DFT}[x\circledast y]=X\cdot Y$ (con le convenzioni del corso, e il fattore $T$ nella convoluzione). Verificato numericamente con vettori casuali.
• Forma matriciale. Con $W$ la matrice $N\times N$ di elementi $W_{kn}=e^{-i2\pi kn/N}$ (matrice di Fourier), il vettore $\mathbf S=T\,W\mathbf s$. Per $N=4$ le righe sono $(1,1,1,1)$, $(1,-i,-1,i)$, $(1,-1,1,-1)$, $(1,i,-1,-i)$. Poiché $WW^H=N\,I$ le righe sono ortogonali e la matrice $\frac1{\sqrt N}W$ è unitaria: la "forma unitaria" della DFT è $\hat{\mathbf s}=\frac1{\sqrt N}W\mathbf s$, che conserva la norma (Basi ortonormali e Gram-SchmidtUna base ortonormale è fatta di vettori di norma 1 a due a due ortogonali: le coordinate si calcolano con prodotti scalari e le proiezioni con una formula diretta. Il procedimento di Gram-Schmidt trasforma una base qualsiasi in una base ortogonale (e poi ortonormale) dello stesso sottospazio.Basi ortonormali e Gram-Schmidt →, Prodotto scalare, norma e angoliIl prodotto scalare aggiunge a uno spazio vettoriale lunghezze e angoli: norma, disuguaglianza di Cauchy-Schwarz, angolo tra vettori in R^n, ortogonalità, proiezione su una retta, aree e volumi con il determinante della matrice dei prodotti scalari.Prodotto scalare, norma e angoli →).

Che cosa rappresenta la DFT

1. Segnale periodico di partenza discreto

È il caso diretto: i $S_k$ sono i coefficienti della serie di Fourier del segnale discreto, le uniche $N$ frequenze distinguibili.

2. Campionamento di un segnale continuo periodico

Se i campioni $s(nT)$, $n=0,\dots,N-1$, vengono da un segnale continuo periodico $s_0$ di periodo $NT$ con coefficienti $S^0_h$, i coefficienti del campionato sono $$S_k=\sum_{m=-\infty}^{+\infty}S^0_{k+mN}:$$ la ripetizione periodica dei coefficienti con periodo $N$ (è il teorema di campionamento in forma di serie, Campionamento e ricostruzioneIl campionamento $s_c(nT)=s(nT)$ trasforma un segnale continuo in uno discreto e, in frequenza, ripete lo spettro con periodo $F_c=1/T$: $S_c(f)=\sum_kS(f-kF_c)$. Se le repliche si sovrappongono si ha aliasing e il segnale non è recuperabile. Un interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g$ produce $\tilde s(t)=\sum_nT g(t-nT)s(nT)$ e $\tilde S=G,S_c$. Teorema del campionamento: se $s$ ha banda $B$ e $F_c\ge2B$, l'interpolatore ideale $G=\operatorname{rect}(f/F_c)$ ricostruisce esattamente $s(t)=\sum s(nT)\operatorname{sinc}(F_c(t-nT))$. Se le ipotesi non valgono c'è un errore, in banda e fuori banda, riducibile con un prefiltro anti-aliasing.Campionamento e ricostruzione →). Se $s_0$ contiene solo armoniche $|h|<N/2$ non c'è sovrapposizione e $S_k=S^0_k$ per $|k|<N/2$ (esatto).

Esempio. L'onda quadra di periodo $1$ (vale $1$ per $|t|<\frac14$, coefficienti $S^0_h=\frac12\operatorname{sinc}\frac h2$) campionata con $N=8$ (valori $1,1,\frac12,0,0,0,\frac12,1$ ai tempi $0,\frac18,\dots$; i campioni nei salti valgono $\frac12$) dà $S_0=0{,}5$, $S_1=0{,}3018$, $S_3=-0{,}0518$, $S_2=S_4=0$: coincidono con $\sum_mS^0_{k+8m}$. Per esempio $S_1$ vale $0{,}3018$ e non $S^0_1=0{,}3183$ perché si sommano i contributi delle armoniche $-7$ e $9$, $-15$ e $17$, ...: $0{,}3183-0{,}0455+0{,}0354-\dots$ (aliasing).

3. Campioni della trasformata di un segnale a durata limitata

Se $s(nT)$ è un segnale discreto aperiodico con estensione contenuta in $N$ campioni, lo si considera il periodo di un segnale periodico (le repliche non si sovrappongono nel tempo). Allora $$S(kF)=\sum_{n=0}^{N-1}T\,s(nT)e^{-i2\pi kn/N}=S(f)\Big|_{f=kF},$$ cioè la DFT è un campionamento della trasformata $S(f)$ del segnale originario (quella di 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 →) nei punti $kF=\frac k{NT}$.

Esempio. Il rect con $4$ campioni unitari ($T=1$) ha $S(f)=\frac{\sin4\pi f}{\sin\pi f}e^{-i3\pi f}$. Con $N=4$ la DFT dà solo $4$ valori ($4,0,0,0$). Aggiungendo zeri fino a $N=16$ (zero-padding) si ottengono $16$ campioni della stessa $S(f)$ in $f=\frac k{16}$: la DFT lunga $16$ coincide con la formula (verificato), perché il segnale è sempre lo stesso con più zeri.

Lo zero-padding non aumenta la risoluzione (due componenti vicine restano indistinguibili se la durata del segnale è breve) ma infittisce i punti in cui si valuta la trasformata, e quindi rende il grafico più liscio. Per spaziare di più i campioni in frequenza ($F=\frac1{NT}$ più piccolo) si deve avere $N$ più grande, cioè estendere il segnale nel tempo.

4. Approssimazione della trasformata di un segnale continuo

Per un segnale continuo $s_0(t)$ aperiodico a durata (circa) limitata in $[-T_1,T_1)$, lo si campiona con passo $T=\frac{2T_1}N$ e si usa $S_0(kF)\approx T\,\mathrm{fft}$ con $F=\frac1{NT}=\frac1{2T_1}$ (come nei laboratori, Esercizio - laboratorio 4, segnali continui aperiodici e risoluzione). Due scelte progettuali:

• il passo in frequenza $F=\frac1{NT}=\frac1{2T_1}$ (risoluzione) dipende solo dalla durata osservata $2T_1$;
• l'estensione dell'asse delle frequenze, $\pm\frac{F_p}2=\pm\frac1{2T}$, dipende dal passo di campionamento: se il segnale non è a banda limitata entro quel valore c'è aliasing.

FFT: come si calcola

Il calcolo diretto di $N$ valori con $N$ moltiplicazioni ciascuno richiede circa $N^2$ operazioni. La FFT (fast Fourier transform) abbassa il costo a circa $N\log_2N$. L'idea (decimazione nel tempo, per $N$ potenza di $2$): si separano i campioni di indice pari e dispari, $$X_k=\underbrace{\sum_m x_{2m}e^{-i2\pi km/(N/2)}}_{E_k}+e^{-i2\pi k/N}\underbrace{\sum_m x_{2m+1}e^{-i2\pi km/(N/2)}}_{O_k},\qquad X_{k+N/2}=E_k-e^{-i2\pi k/N}O_k,$$ dove $E_k$ e $O_k$ sono DFT di lunghezza $N/2$, e si ripete la divisione fino a lunghezza $1$. Per $N=1024$: circa $5120$ moltiplicazioni complesse contro oltre $10^6$. Il risultato è identico alla DFT diretta (verificato: l'implementazione ricorsiva coincide con np.fft.fft per $N=16$).

Convenzioni numeriche (NumPy / MATLAB)

np.fft.fft(x) (come fft di MATLAB) calcola $\sum_{n=0}^{N-1}x_ne^{-i2\pi kn/N}$ senza fattori. Quindi con le convenzioni del corso:

I primi $N/2$ valori sono le frequenze positive $0,\dots,\frac{F_p}2$, gli altri le frequenze negative (per la periodicità $k\equiv k-N$). Se il segnale discreto non parte da $n=0$ ma è centrato in $0$ (segnale non causale), la fase della DFT contiene un termine lineare dovuto allo spostamento: si ottiene la fase attesa riordinando le due metà del vettore (equivale a ifftshift) prima della FFT. Il leakage compare se il segnale periodico non ha un numero intero di periodi nella finestra: la DFT di un coseno con $f_0$ non multiplo di $F$ non è una sola riga ma si distribuisce su molti $k$.

Esercizi collegati

• Esercizio - DFT di un segnale discreto periodico
• Esercizio - laboratorio 3, DFT di segnali periodici
• Esercizio - laboratorio 4, coefficienti di Fourier con la FFT
• Esercizio - laboratorio 4, trasformata di segnali discreti con la FFT

Vedi anche la materia gemella: Trasformata di Fourier discreta (TFD)La TFD è la serie di Fourier dei segnali discreti periodici di periodo $N$: servono solo $N$ armoniche, $X(k)=\frac1N\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}$ e $x(n)=\sum_{k=0}^{N-1}X(k)e^{j2\pi kn/N}$. È un prodotto matrice-vettore in $\mathbb{C}^N$ (la matrice di Fourier); la convoluzione circolare diventa prodotto, $N X(k)Y(k)$. La FFT la calcola in $O(N\log N)$.Trasformata di Fourier discreta (TFD) → e 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 →.