FFT e zero-padding - TF, TFtd e asse delle pulsazioni

Una domanda di calcolo numerico compare in ogni tema d'esame di Segnali e Sistemi: scrivere uno script che calcoli la trasformata di Fourier di un segnale campionato con la FFT, con l'asse delle pulsazioni corretto e un "opportuno zero-padding". Questa nota spiega perché si fa così, ed è corredata da codice Python eseguito (equivalente al Matlab degli esami: vedi la tabella in 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 →).

Dal segnale continuo alla FFT

Sia $x(t)$ un segnale continuo, nullo fuori da un intervallo $[t_0,t_0+T_{obs}]$ (o trascurabile fuori), campionato con passo $T_c$ negli istanti $t_n=t_0+nT_c$, $n=0,\ldots,N-1$, con campioni $x_n$. Si approssima l'integrale della TF con la somma rettangolare: $$X(\omega)=\int x(t)e^{-j\omega t}dt\approx T_c\sum_{n=0}^{N-1}x_n\,e^{-j\omega t_n}=T_c\,e^{-j\omega t_0}\sum_{n=0}^{N-1}x_n\,e^{-j\omega nT_c}.$$ L'ultima somma è la TFtd della sequenza dei campioni calcolata in $\theta=\omega T_c$ (Trasformata di Fourier a tempo discreto (TFtd)La TFtd di una sequenza è $X(\omega)=\sum_nx(n)e^{-j\omega n}$, funzione continua e periodica di periodo $2\pi$; si inverte con $x(n)=\frac1{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{j\omega n}d\omega$. Ha le stesse proprietà della TF continua (la convoluzione diventa prodotto, $n,x(n)\leftrightarrow jX'$). È la risposta in frequenza dei sistemi discreti; con la TFD e lo zero-padding se ne ottengono campioni arbitrariamente fitti.Trasformata di Fourier a tempo discreto (TFtd) →): $\hat X(\theta)=\sum x_ne^{-j\theta n}$. Quindi $$\boxed{X(\omega)\approx T_c\,e^{-j\omega t_0}\,\hat X(\omega T_c)}$$ La somma è anche il campionamento di $X$ visto in Campionamento e formula di PoissonCampionare un segnale continuo $x(t)$ con passo $T_c$ dà la sequenza $x(nT_c)$. La formula di Poisson lega gli spettri: $\hat W(\omega)=\frac1{T_c}\sum_kX\left(\frac{\omega+2\pi k}{T_c}\right)$, cioè lo spettro del segnale campionato è la ripetizione periodica (di periodo $2\pi/T_c$ in pulsazione analogica) dello spettro originale, riscalata. Se le repliche si sovrappongono si ha aliasing.Campionamento e formula di Poisson → (la TFtd dei campioni è la ripetizione periodica di $X$ riscalata: $\hat X(\theta)=\frac1{T_c}\sum_kX\left(\frac{\theta+2\pi k}{T_c}\right)$), quindi l'approssimazione è buona se l'aliasing è trascurabile (lo spettro è piccolo oltre $\frac\pi{T_c}$) e se il segnale è davvero nullo fuori dall'intervallo osservato.

Cosa calcola la FFT

np.fft.fft(x, M) ($=$ fft(x,M) in Matlab) calcola la TFD (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) →) di $M$ punti del vettore completato con $M-N$ zeri. Con la convenzione della libreria, senza il fattore $\frac1M$: $$\mathtt{fft}(x,M)[k]=\sum_{n=0}^{N-1}x_n\,e^{-j\frac{2\pi}{M}kn}=\hat X\left(\frac{2\pi k}M\right),\qquad k=0,\ldots,M-1.$$ Sono $M$ campioni equispaziati della TFtd $\hat X(\theta)$ nel periodo $[0,2\pi)$, ai valori $\theta_k=\frac{2\pi k}M$. Con $\theta=\omega T_c$ ($\omega_k=\frac{\theta_k}{T_c}$) si ottengono le pulsazioni:

$$\omega_k=\frac{2\pi k}{M\,T_c},\qquad f_k=\frac{k}{M\,T_c},\qquad X(\omega_k)\approx T_c\,e^{-j\omega_kt_0}\,\mathtt{fft}(x,M)[k].$$

Il passo in pulsazione è $\Delta\omega=\frac{2\pi}{MT_c}$ e il massimo valore rappresentabile è $\frac{\omega_c}2=\frac\pi{T_c}$ (metà della pulsazione di campionamento $\omega_c=\frac{2\pi}{T_c}$): oltre c'è aliasing.

Frequenze negative. Gli indici $k>M/2$ corrispondono, per periodicità, a pulsazioni negative: $\omega_k\equiv\omega_k-\frac{2\pi}{T_c}$. Per avere l'asse centrato in zero si riordina con fftshift e si costruisce l'asse da $-\frac\pi{T_c}$ a $\frac\pi{T_c}-\Delta\omega$:

M = 2**(nextpow2(N) + 3)                 # zero-padding
X = np.fft.fftshift(Tc*np.fft.fft(x, M))
w = 2*np.pi*np.arange(-M/2, M/2)/(M*Tc)  # asse delle pulsazioni [rad/s]

(in Matlab: w = -pi/Tc : 2*pi/(M*Tc) : pi/Tc - 2*pi/(M*Tc)). Con nextpow2 che è $\lceil\log_2N\rceil$: nextpow2 = lambda N: int(np.ceil(np.log2(N))).

Fattore di fase. Se l'asse dei tempi parte da $t_0\ne0$ lo spettro calcolato dalla FFT è quello del segnale traslato a partire da $0$: il modulo $|X|$ è identico (un ritardo non lo cambia) e la fase differisce per $e^{-j\omega t_0}$. Nei temi d'esame si traccia quasi sempre il solo $|X|$, quindi non serve.

Se serve una frequenza in Hz (e non una pulsazione): $f_k=\frac{k}{MT_c}$ con $f\in\left[-\frac{F_c}2,\frac{F_c}2\right)$, $F_c=\frac1{T_c}$ (tema d'esame febbraio 2024: "attenzione, le ascisse sono frequenze, non pulsazioni").

Esempio: impulso modulato con la FFT

Segnale $x(t)=A\cos(\omega_0t)\operatorname{rect}\left(\frac{t-T}{2T}\right)$ con $A=2$, $T=10$, $\omega_0=3$ (cioè un coseno di durata $20$ s tra $0$ e $20$), calcolato in $(-2,30)$ con $T_c=10^{-3}$ (tema d'esame gennaio 2025). Teoricamente $x=A\cos(\omega_0t)w(t)$ con $w=\operatorname{rect}\left(\frac{t-T}{2T}\right)$, $W(\omega)=2T\operatorname{sinc}\left(\frac{\omega T}\pi\right)e^{-j\omega T}$ e $$X(\omega)=\frac A2\left[W(\omega-\omega_0)+W(\omega+\omega_0)\right],$$ con un picco di modulo $\approx\frac A2\cdot2T=AT=20$ in $\omega=\pm3$ (le due copie sono ben separate).

A, T, w0, Tc = 2.0, 10.0, 3.0, 1e-3
t = np.arange(-2, 30 + Tc/2, Tc)
x = A*np.cos(w0*t) * (np.abs((t - T)/(2*T)) < 0.5)
N = len(t); M = 2**(3 + nextpow2(N))
X = Tc * np.fft.fft(x, M)
w = 2*np.pi*np.arange(M)/(M*Tc)

Risultati ($N=32001$, $M=262144$): il picco di $|X|$ vale $20{,}08$ in $\omega=2{,}996$, e in tre punti prossimi al picco la FFT e la formula esatta danno $|X(2{,}5)|=3{,}345$ contro $3{,}345$, $|X(3)|=20{,}080$ contro $20{,}080$, $|X(3{,}5)|=4{,}091$ contro $4{,}091$ (ottimo accordo).

Esempio: confronto con lo spettro teorico

$x(t)=2e^{-3t}u(t)$ osservato in $[0,20]$ con $T_c=10^{-3}$ (tema d'esame febbraio 2026). Teoria: $X(\omega)=\frac{2}{3+j\omega}$. Il segnale non è a banda limitata, ma decade come $\frac1{|\omega|}$ e $\frac{\pi}{T_c}\approx3142\gg3$:

Tc = 1e-3; t = np.arange(0, 20 + Tc/2, Tc); x = 2*np.exp(-3*t)
N = len(t); M = 2**(nextpow2(N) + 3)
X = np.fft.fftshift(Tc*np.fft.fft(x, M))
w = 2*np.pi*np.arange(-M/2, M/2)/(M*Tc)
Xe = 2/(3 + 1j*w)

Per $|\omega|<20$ l'errore massimo $|X_{\text{FFT}}-X_{\text{esatta}}|$ è $1{,}0\cdot10^{-3}$, contro un picco $|X(0)|=\frac23$: errore relativo $0{,}15\%$. La parte principale è dovuta al valore $x(0)=2$ del primo campione: con la convenzione dell'emivalore il primo campione dovrebbe pesare $\frac12$ (errore $\approx\frac{T_c}{2}\cdot2=10^{-3}$, esattamente lo scarto trovato). L'aliasing contribuisce molto meno ($\approx\frac{2}{6283}$ per la prima replica a $\omega=0$, Conversione A-D e campionamento di segnali non limitati in bandaUn convertitore A/D è la catena filtro anti-aliasing passa-basso, campionatore, quantizzatore a $B$ bit. Nessun segnale reale è a banda esattamente limitata: si sceglie la banda essenziale (che contiene quasi tutta l'energia) e un passo $T_c$ tale che le repliche aliasing siano trascurabili. La quantizzazione introduce un errore con potenza $\Delta^2/12$ e un rapporto segnale-rumore di circa $6{,}02B+1{,}76$ dB.Conversione A-D e campionamento di segnali non limitati in banda →).

Zero-padding: cosa fa e cosa non fa

Zero-padding significa passare da $N$ a $M>N$ punti aggiungendo $M-N$ zeri in coda. Poiché il segnale è nullo comunque fuori da $[t_0,t_0+T_{obs}]$, la TFtd $\hat X(\theta)$ non cambia: cambia solo il numero di punti $M$ in cui viene campionata, con passo $\frac{2\pi}M$. Dimostrazione (Trasformata di Fourier a tempo discreto (TFtd)La TFtd di una sequenza è $X(\omega)=\sum_nx(n)e^{-j\omega n}$, funzione continua e periodica di periodo $2\pi$; si inverte con $x(n)=\frac1{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{j\omega n}d\omega$. Ha le stesse proprietà della TF continua (la convoluzione diventa prodotto, $n,x(n)\leftrightarrow jX'$). È la risposta in frequenza dei sistemi discreti; con la TFD e lo zero-padding se ne ottengono campioni arbitrariamente fitti.Trasformata di Fourier a tempo discreto (TFtd) →): per $x$ a supporto in $\{0,\ldots,N-1\}$, la TFD della ripetizione periodica di periodo $M\ge N$ vale $\frac1M\hat X\left(\frac{2\pi k}M\right)$, quindi scegliendo $M$ grande a piacere la risoluzione con cui si "disegna" $\hat X$ è arbitrariamente fine (tema d'esame ricorrente: mostrare che con la TFD e lo zero-padding si ottiene una rappresentazione di $\hat X$ con risoluzione arbitrariamente fine).

• Serve a: disegnare una curva liscia (senza lo zero-padding $M=N$ pochi punti nascondono i lobi del sinc), trovare con precisione la posizione di un picco, avere una potenza di 2 per la FFT.
• Non serve a: distinguere due frequenze vicine. La risoluzione (capacità di separare due righe) dipende dalla durata osservata $T_{obs}=NT_c$: due sinusoidi a distanza $\Delta f$ si separano solo se $\Delta f\gtrsim\frac1{T_{obs}}$, perché il segnale è moltiplicato per una finestra rettangolare lunga $T_{obs}$, che nello spettro è un sinc con lobo principale largo $\sim\frac{2}{T_{obs}}$ Hz.

Verifica numerica. Due coseni a $10$ Hz e $10{,}5$ Hz ($\Delta f=0{,}5$ Hz), $T_c=10^{-3}$:

Con $T_{obs}=1$ s nessuno zero-padding separa le righe; aumentando il tempo osservato appaiono. (Le posizioni dei massimi sono leggermente alterate dalla sovrapposizione dei due lobi.)

Simmetria dello spettro di un segnale reale

Per un segnale reale $|X(\omega)|$ è pari (Proprietà della trasformata di FourierCon le proprietà (linearità, simmetrie, ritardo $\leftrightarrow e^{-j\omega t_0}$, modulazione $\leftrightarrow$ traslazione in frequenza, scala, dualità, convoluzione $\leftrightarrow$ prodotto, Parseval $E=\frac1{2\pi}\int|X|^2$, derivata $\leftrightarrow j\omega$, moltiplicazione per $t\leftrightarrow j,d/d\omega$, integrazione) quasi tutte le trasformate si ottengono da poche coppie base senza integrare.Proprietà della trasformata di Fourier →) e quindi, per i campioni, $|\mathtt{fft}(x,M)[k]|=|\mathtt{fft}(x,M)[M-k]|$ per $k=1,\ldots,M-1$ (indici da $0$). In Matlab (indici da $1$): $|X(k)|=|X(M-k+2)|$. Esempio d'esame (febbraio 2023): "se il massimo di $|X|$ è nel 350-esimo campione, dove ci si aspetta un altro picco?". Risposta: nel campione $M-350+2=M-348$ (accettata anche $M-350$), perché l'indice $k$ (da 1) corrisponde a $\omega=(k-1)\frac{2\pi}{MT_c}$ e $-\omega\equiv\omega_c-\omega$ è l'indice $M-k+2$. Con indici da $0$ (Python) la corrispondenza è $p\leftrightarrow M-p$ (verificata numericamente per un segnale reale casuale con $p=350$, $M=4096$).

Dall'antitrasformata: da $X$ a $x$

Il procedimento inverso (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) →): se X contiene i campioni di $X(\omega)$ sulle pulsazioni $\omega_0\cdot(-\frac N2:\frac N2-1)$ (con $N$ pari), il passo temporale è $T_c=\frac{2\pi}{N\omega_0}$ e il segnale è x = ifft(ifftshift(X))/Tc con tempi t = (0:N-1)*Tc (compito del canale Ing. Biomedica). Controllato: partendo da un esponenziale campionato, fft $\to$ fftshift $\to$ ifftshift $\to$ ifft restituisce esattamente il segnale (True, con $N=2000$).

Prove al calcolatore: tre compiti tipici

Oltre alla domanda dello scritto, il corso prevede una prova al calcolatore facoltativa (punti bonus). I compiti dei due ultimi anni sono di tre tipi.

1. Somma parziale della serie di Fourier di un'onda quadra (prova del 2023). Si traccia l'onda quadra $x(t)=\sum_k\operatorname{rect}\left(\frac{t-kT}{T/2}\right)$ (con $T=1$) in $\left(-\frac12,\frac12\right)$, dove coincide con $\operatorname{rect}(2t)$, e la somma parziale $x_M(t)=\frac12+\frac2\pi\sum_{m=0}^{M-1}\frac{(-1)^m}{2m+1}\cos\left[(2m+1)\omega_0t\right]$ con un ciclo for che aggiunge una armonica per volta, inizializzando $x_M$ alla costante $\frac12$. Il codice Python è nella nota 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 → (somme parziali e fenomeno di Gibbs): le armoniche sono quelle dispari di coefficiente $a_k=\frac{\sin(k\pi/2)}{k\pi}$ (Serie di Fourier - analisi e sintesiUn segnale periodico di periodo $T$ si scrive come somma di esponenziali in relazione armonica, $x(t)=\sum_ka_ke^{jk\omega_0t}$ con $\omega_0=2\pi/T$, e i coefficienti si ottengono per proiezione, $a_k=\frac1T\int_Tx(t)e^{-jk\omega_0t}dt$. L'ortogonalità degli esponenziali dà la formula; la convergenza è in media quadratica (Riesz-Fischer), con il fenomeno di Gibbs nei salti.Serie di Fourier - analisi e sintesi →).

2. Spettro di un segnale audio e frequenza principale (prova del 2024). Dati i campioni y e la frequenza di campionamento $F_s$ ($T_s=\frac1{F_s}$): zero-padding $M=2^{\text{nextpow2}(N)+2}$ (cioè $M\ge4N$), Y = Ts*fft(y, M), asse delle frequenze positive F = np.arange(M//2)*Fs/M (passo $\frac{F_s}{M}$, fino a $\frac{F_s}2$), grafico di $|Y(\omega)|$ in funzione di $\omega=2\pi f$, di $|Y(f)|$ e in scala logaritmica ($10\log_{10}|Y|$ o semilogy), e frequenza principale F[np.argmax(np.abs(Y[:M//2]))] in Hz. Prova su un tono sintetico ($440{,}3$ Hz con smorzamento e rumore, $F_s=8192$ Hz, 2 s, $M=65536$): il massimo cade a $440{,}25$ Hz, entro il passo di griglia $\frac{F_s}M=0{,}125$ Hz.

3. Stima della frequenza di un esponenziale complesso e risoluzione (prova del 2026). Si genera $x(n)=e^{j2\pi\nu n}$ con $\nu=0{,}1234567$, $N=64$ campioni, si calcola fft(x, M) con zero-padding e si stima $\nu$ con l'indice del massimo, $\hat\nu=\frac{k_{\max}}M$. Per un segnale complesso lo spettro non è pari e va guardato l'intero intervallo $\left(-\frac12,\frac12\right)$ (con fftshift). L'errore è al più mezzo passo di griglia, $\frac1{2M}$: con $M=2048$ si trova $\hat\nu=0{,}123535$ (errore $7{,}8\cdot10^{-5}\le2{,}4\cdot10^{-4}$); con $M=N=64$ solo $\hat\nu=0{,}125$ (errore $1{,}5\cdot10^{-3}\le7{,}8\cdot10^{-3}$); con $M=16384$ l'errore scende a $1{,}7\cdot10^{-5}$: lo zero-padding affina la stima del picco di un tono isolato, perché la forma del lobo è nota e liscia. Poi si aggiunge una seconda componente a $\nu+\Delta\nu$ con $\Delta\nu=0{,}015$ e ampiezza $\frac12$ e ci si chiede con quanti campioni $N_1$ si distinguono i due picchi: la finestra rettangolare ha lobo principale con zeri a $\pm\frac1{N}$, quindi serve circa $\Delta\nu\gtrsim\frac1{N_1}$ (criterio di Rayleigh, $N_1\gtrsim67$); la soluzione del corso è più prudente ($\frac2{N_1}<\Delta\nu$, $N_1>134$). Numericamente: con $N_1=32$ si vede un picco, con $N_1\ge48$ due ($M=2^{14}$, soglia $0{,}35$ del massimo). Lo zero-padding, anche enorme, non separa i picchi se $N_1$ è troppo piccolo: conta la durata osservata.

Procedura per l'esame

1. Definisci Tc, l'asse t, il segnale x (con rect = lambda t: abs(t)<0.5).
2. N = len(t); M = 2**(nextpow2(N)+3) (potenza di due, almeno $8N$; l'esame accetta $2N$–$4N$ o "opportuno").
3. X = Tc*fft(x, M) (il fattore $T_c$ è essenziale per un segnale continuo).
4. Asse: w = 2*pi*(0:M-1)/(M*Tc) oppure, con fftshift, da $-\pi/T_c$ a $\pi/T_c$ con passo $\frac{2\pi}{MT_c}$. Etichetta in rad/s (o in Hz dividendo per $2\pi$).
5. Grafico del modulo abs(X) e limiti d'asse con xlim.

Errori comuni

• Dimenticare il fattore $T_c$ (lo spettro risulta $\frac1{T_c}$ volte più grande).
• Dimenticare che fft non ha il fattore $\frac1M$ (e sbagliare il confronto con la TFD "con $\frac1N$" di 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) →).
• Costruire l'asse con $\frac{2\pi k}{N}$ invece di $\frac{2\pi k}{MT_c}$ (si dimentica il passo di campionamento e lo zero-padding).
• Non usare fftshift e disegnare lo spettro "con le frequenze negative a destra".
• Credere che lo zero-padding migliori la risoluzione: aumenta solo il numero di punti.