Analisi spettrale con la DFT - finestre e leakage

L'analisi di Fourier (o spettrale) è l'applicazione principale della DFT: dato un segnale $x[n]$ si vuole capire quali frequenze contiene e con quale ampiezza. In teoria servirebbe la DTFT di $x[n]$ su tutto $n\in\mathbb Z$ (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 →); in pratica si dispone solo di un tratto finito di $L$ campioni e si può calcolare solo la DFT (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 →). Le due approssimazioni, finestratura nel tempo e campionamento in frequenza, producono artefatti che bisogna saper riconoscere: lobi laterali, perdita di risoluzione, errore di ampiezza. Il ragionamento è lo stesso per le immagini e per l'audio; per segnali la cui frequenza cambia nel tempo si usa la STFT (Trasformata di Fourier a tempo breve e spettrogrammaPer un segnale lungo il cui contenuto in frequenza cambia nel tempo si calcolano tante DFT brevi: X[k, l] = sum_{n=0}^{L-1} w[n] x[l+n] e^{-j2pi kn/N}, con finestra di analisi w di lunghezza L, istante l = m R (passo R, di solito L/2) e DFT a N >= L punti. Il modulo (o i dB) di X[k, l] come immagine tempo-frequenza è lo spettrogramma. La finestra impone un compromesso: lunga = ottima risoluzione in frequenza (Delta w circa 8pi/L per Hann) ma transizioni temporali sfocate su circa L/2 campioni; corta = buona localizzazione nel tempo ma righe vicine fuse. Si provano più lunghezze.Trasformata di Fourier a tempo breve e spettrogramma →).

L'effetto della finestra

Un esponenziale complesso

Sia $x[n]=Ae^{j\hat\omega_0n}$, $n\in\mathbb Z$, con $A$ costante complessa. La sua DTFT è un treno di impulsi periodico: $X(e^{j\hat\omega})=2\pi A\sum_{r\in\mathbb Z}\delta(\hat\omega-\hat\omega_0+2\pi r)$. Dalla posizione dell'impulso si legge $\hat\omega_0$ e dalla sua area $A$: ogni parametro è determinato. Se però si osservano solo $L$ campioni, il segnale osservato è $$y[n]=x[n]\,w_R[n],\qquad w_R[n]=u[n]-u[n-L]\quad(\text{finestra rettangolare di lunghezza }L).$$ La moltiplicazione nel tempo diventa convoluzione periodica in frequenza (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 l'impulso si converte in una copia della trasformata della finestra: $$Y(e^{j\hat\omega})=A\,W_R\big(e^{j(\hat\omega-\hat\omega_0)}\big),\qquad W_R(e^{j\hat\omega})=\sum_{n=0}^{L-1}e^{-j\hat\omega n}=e^{-j\hat\omega\frac{L-1}2}\,\frac{\sin(L\hat\omega/2)}{\sin(\hat\omega/2)}.$$ L'ultima uguaglianza è la somma di una serie geometrica di ragione $e^{-j\hat\omega}$ (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 →), con lo stesso passaggio di $e^{-j\hat\omega L/2}\big(e^{j\hat\omega L/2}-e^{-j\hat\omega L/2}\big)/\big(e^{-j\hat\omega/2}(e^{j\hat\omega/2}-e^{-j\hat\omega/2})\big)$. Il rapporto di seni è il nucleo di Dirichlet (Dirichlet kernel), una "sinc periodica".

Fatto (spettro di una sinusoide osservata per $L$ campioni). Il picco di $\lvert Y\rvert$ è in $\hat\omega_0$ e vale $L\lvert A\rvert$; $\lvert W_R\rvert$ si annulla in $\hat\omega=2\pi r/L$ ($r\neq0$). La parte tra i due zeri $\pm2\pi/L$ è il lobo principale (main lobe), di larghezza $4\pi/L$; le parti tra zeri consecutivi successivi sono i lobi laterali (side lobes), larghi $2\pi/L$. Il primo lobo laterale vale circa il $22\%$ del principale, cioè $-13.3\ \mathrm{dB}$, qualunque sia $L$.

Esempio. Per $L=32$ e $A=1$: picco $32$, zeri a multipli di $2\pi/32=0.196$ rad, lobo principale largo $4\pi/32=0.393$ rad. Il massimo del primo lobo laterale è $\approx0.217\cdot32\approx6.9$ (calcolo numerico: $-13.3\ \mathrm{dB}$ per $L$ grande).

Aumentando $L$, il lobo principale si restringe (larghezza $\propto1/L$) e si avvicina all'impulso ideale, ma l'ampiezza relativa dei lobi laterali non diminuisce.

Più componenti: risoluzione e leakage

Se $x[n]=A_1e^{j\hat\omega_1n}+A_2e^{j\hat\omega_2n}$, per linearità $Y=A_1W_R(\hat\omega-\hat\omega_1)+A_2W_R(\hat\omega-\hat\omega_2)$: le due copie del nucleo si sommano e interferiscono. Ne derivano due problemi.

• Risoluzione (resolution): se $\hat\omega_1$ e $\hat\omega_2$ sono troppo vicine, i due lobi principali si fondono in un solo picco e non si riconoscono due righe. Regola pratica: due toni si separano se distano più di circa la mezza larghezza del lobo principale, cioè $2\pi/L$ per la finestra rettangolare. La risoluzione dipende solo dalla durata $L$ dell'osservazione e dalla finestra, non dallo zero-padding.
• Leakage (spectral leakage, "perdita" di energia): l'energia di un tono sfugge nei lobi laterali e può coprire (mask) un tono più debole vicino, rendendo errata anche la stima delle ampiezze.

Per ridurre l'interferenza si prende $L$ più grande possibile (lobi principali stretti) oppure si cambia la finestra, come sotto.

Esempio (lab 3). $x(t)=\cos(2\pi f_1t)+10^{-3}\cos(2\pi f_2t)+\cos(2\pi f_3t)$ con $f_1=1.5$, $f_2=2.5$, $f_3=3.5\ \mathrm{kHz}$ e $F_s=10\ \mathrm{kHz}$: il tono centrale è $20\log_{10}10^{-3}=-60\ \mathrm{dB}$ rispetto agli altri. Per vederlo bisogna che i lobi laterali dei due toni forti stiano sotto $-60\ \mathrm{dB}$ (con un margine, per esempio $-70\ \mathrm{dB}$): non lo fa né la finestra rettangolare ($-13\ \mathrm{dB}$) né quella di Hamming ($-42\ \mathrm{dB}$).

Il campionamento in frequenza

Calcolando la DFT a $N$ punti del segnale osservato si ottengono i valori $Y[k]=Y(e^{j2\pi k/N})$, $k=0,\dots,N-1$: campioni della DTFT a $\hat\omega_k=2\pi k/N$. Un tono a $\hat\omega_0$ qualsiasi non coincide in generale con uno di questi punti, e i campioni possono mancare il vertice del lobo principale: l'ampiezza stimata è sbagliata.

Condizione di picco corretto. Con $N=L$ (DFT di lunghezza pari alla finestra), l'ampiezza è letta bene solo se $\hat\omega_0=\frac{2\pi}{N}k_0$ con $k_0$ intero (il tono cade su un bin). Equivalente: la finestra contiene un numero intero $k_0$ di periodi del tono. Allora il lobo principale è campionato nel vertice e tutti gli altri campioni cadono sugli zeri: la DFT è una riga sola.

La ragione: per $\hat\omega_0=2\pi k_0/N$ i campioni $x[0]=1$ e $x[N]=e^{j2\pi k_0}=1$ coincidono, quindi la ripetizione periodica del tratto osservato non ha salti e il tono è esattamente periodico di periodo $N$.

Esempio (verificato con Python). $L=N=32$, $x[n]=e^{j2\pi k_0n/32}$.

• $k_0=4$ (su un bin): $\lvert X[k]\rvert=32$ per $k=4$ e $0$ per tutti gli altri.
• $k_0=4.5$ (a metà tra due bin): $\lvert X[4]\rvert=\lvert X[5]\rvert=20.38$, $\lvert X[3]\rvert=\lvert X[6]\rvert=6.82$ e tanti altri valori non nulli (leakage). Il picco campionato è $20.38$ invece di $32$: perdita di $-3.92\ \mathrm{dB}$ (scalloping loss).
• $k_0=4.25$: picco $28.81$ ($-0.91\ \mathrm{dB}$).

Zero-padding

Per campionare la DTFT più fitta si completa il segnale con zeri e si calcola una DFT a $N>L$ punti. Il massimo viene campionato più vicino al vertice e l'ampiezza è stimata meglio. Con $L=32$, $k_0=4.5$ (cioè $\hat\omega_0=2\pi\cdot4.5/32$):

Con $N=64$ il punto $\hat\omega_0=2\pi\cdot9/64$ è proprio un campione. Lo zero-padding non aggiunge informazione: i campioni sono punti della stessa curva (il nucleo di Dirichlet), quindi non separa due toni che la finestra non separa, ma permette di localizzare meglio il picco.

Finestre rastremate

La finestra rettangolare ha il lobo principale più stretto (miglior risoluzione) ma i lobi laterali più alti. Una finestra rastremata (tapered window) scende dolcemente a zero ai bordi e ha lobi laterali molto più bassi, a prezzo di un lobo principale più largo: minor leakage, minor risoluzione.

Compromesso (leakage contro risoluzione). Lobi laterali più bassi si pagano con un lobo principale più largo; non esiste una finestra che migliori entrambi. Si prova più finestre e più lunghezze e si confrontano gli spettri.

Con $w[n]$ simmetrica di lunghezza $L$ ($n=0,\dots,L-1$), centro in $(L-1)/2$:

I valori sono stati ricalcolati numericamente su finestre con $L=64$ (FFT a $2^{16}$ punti); i livelli dei lobi laterali dipendono poco da $L$, le larghezze sono $4,\ 8,\ 8,\ 8,\ 12$ volte $\pi/L$. La Hann è una finestra "coseno rialzato" al $100\%$; la Hamming ha un piedistallo ($0.54$) che sposta il primo lobo laterale più in basso ma lascia lobi lontani meno decrescenti; la Blackman aggiunge un termine correttivo. Per i valori di attenuazione in banda proibita e per la finestra regolabile di Kaiser si veda Progetto di filtri FIR con il metodo delle finestreUn filtro FIR a fase lineare di tipo I (ordine $N$ pari, $h[n]=h[N-n]$) ha risposta $H(e^{j\hat\omega})=e^{-j\hat\omega N/2}\bar H(\hat\omega)$ con ampiezza $\bar H(\hat\omega)=\sum_{n=0}^{N/2}p_n\cos(n\hat\omega)$. Per approssimare un'ampiezza desiderata $D(\hat\omega)$ (passa-basso: $1$ in banda passante, $0$ in banda oscura, con tolleranze $\delta_p,\delta_s$ e frequenze $\hat\omega_p,\hat\omega_s$) ci sono tre metodi. Finestre: si tronca la serie di Fourier di $D$, cioè $h[n]=h_d[n-N/2],w[n]$ con $h_d[n]=\frac{\hat\omega_0}{\pi}\operatorname{sinc}\frac{\hat\omega_0 n}{\pi}$; la rettangolare è ottima in errore quadratico ma dà il fenomeno di Gibbs (9%, $21$ dB), le finestre rastremate (Hann, Hamming, Blackman, Kaiser) abbassano i lobi laterali allargando la transizione ($\hat\omega_s-\hat\omega_p\simeq\alpha,2\pi/L$). Campionamento in frequenza: $h=\mathrm{IDFT}$ dei campioni di $D$, esatto solo sui campioni. Minimax (Parks-McClellan): errore pesato minimo nel caso peggiore, soluzione equiripple con almeno $r+2$ alternanze, ordine più basso a parità di specifiche.Progetto di filtri FIR con il metodo delle finestre →.

Le trasformate delle finestre a coseno si ottengono dal nucleo di Dirichlet: poiché $0.5-0.5\cos(\hat\omega_1n)$ con $\hat\omega_1=2\pi/(L-1)$ è la somma di tre termini, nella forma centrata $W(\hat\omega)=0.5\,D(\hat\omega)+0.25\,\big[D(\hat\omega-\hat\omega_1)+D(\hat\omega+\hat\omega_1)\big]$ dove $D(\hat\omega)=\sin(L\hat\omega/2)/\sin(\hat\omega/2)$ è il nucleo senza il fattore di fase $e^{-j\hat\omega(L-1)/2}$ (che i tre termini hanno in comune). I lobi laterali dei tre termini si cancellano in parte, da cui lobi laterali bassi.

Esempio (lab 3, tono debole a $-60\ \mathrm{dB}$). Con $L=217$ campioni (lunghezza ottenuta dalla formula del progetto, $L\ge\frac{A_{dB}-8}{14.36\,\Delta f}F_s+1$ con $A_{dB}=70$, $\Delta f=0.2\ \mathrm{kHz}$, $F_s=10\ \mathrm{kHz}$, che dà $216.9$) e finestra di Kaiser con $\beta=0.1102(70-8.7)=6.755$, il tono a $2.5\ \mathrm{kHz}$ emerge a $-60.9\ \mathrm{dB}$ sopra un leakage di $-67.4\ \mathrm{dB}$. Con la finestra rettangolare e con la Hamming della stessa lunghezza il massimo attorno a $2.5\ \mathrm{kHz}$ è a $-32.5\ \mathrm{dB}$ e a $-48.7\ \mathrm{dB}$, cioè molto sopra il tono da $-60\ \mathrm{dB}$: è leakage dei toni forti (nelle zone vicine il leakage arriva a $-28.5$ e $-45.4\ \mathrm{dB}$) e il tono debole non si distingue (valori calcolati con Python). Il codice è in Esercizio - Laboratorio 3, DFT e analisi spettrale.

Segnali sinusoidali reali

Un segnale sinusoidale reale $x[n]=A\cos(\hat\omega_0n)=\frac A2e^{j\hat\omega_0n}+\frac A2e^{-j\hat\omega_0n}$ è la somma di due esponenziali a pulsazioni opposte ($\hat\omega_0$ e $-\hat\omega_0\equiv2\pi-\hat\omega_0$). Tutto quanto detto vale per ciascuno: servono due coppie di lobi e basta stimare quello in $0\le\hat\omega_0\le\pi$ perché l'altro è il simmetrico (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 →). Con finestra rettangolare e tono su un bin: $\lvert X[k_0]\rvert=\lvert X[N-k_0]\rvert=AL/2$, quindi $A=2\lvert X[k_0]\rvert/L$. Con una finestra generica il picco vale $\frac A2\sum_nw[n]$ e l'ampiezza è $A=2\lvert X\rvert/\sum_nw[n]$ (guadagno coerente).

Esempio. $A=3$, $L=N=32$, $\hat\omega_0=2\pi\cdot5/32$: $\lvert X[5]\rvert=\lvert X[27]\rvert=48$ e $A=2\cdot48/32=3$ ✓. Con una finestra di Hann ($\sum w\approx16$ per $L=32$) si leggerebbe circa $\frac A2\cdot16=24$.

Esempio: due toni vicini

Si osservano $L$ campioni di $x[n]=\cos(0.2\pi n)+1.5\cos(0.227\pi n)$ con una finestra di Hann centrata nell'istante $n_s=2000$ (come nelle slide). La separazione è $0.027\pi$. La finestra di Hann ha lobo principale largo $8\pi/L$, mezza larghezza $4\pi/L$:

• $L=75$: $4\pi/75=0.053\pi>0.027\pi$, i lobi si fondono: un solo picco, a $0.217\pi$, di altezza $39.6$.
• $L=301$: $4\pi/301=0.0133\pi<0.027\pi$: due picchi separati, a $0.200\pi$ (altezza $74.9$) e $0.227\pi$ (altezza $112.4$).

Nota. Se due toni sono più vicini della risoluzione, il risultato dipende anche dalla fase relativa con cui i due lobi si sommano: lo stesso $L=75$ con un'altra posizione della finestra può dare un picco unico oppure due gobbe false a $0.183\pi$ e $0.236\pi$ (interferenza distruttiva a metà), come si verifica numericamente. Per questo si dice che sotto la risoluzione non si può dire se ci siano uno o due toni.

Come procedere in pratica

1. Fissare la risoluzione voluta $\Delta\hat\omega$ e scegliere $L\ge c\cdot2\pi/\Delta\hat\omega$ con $c$ della finestra scelta (circa $1$ per la rettangolare, $2$ per Hann/Hamming, $3$ per Blackman, in modo che il lobo principale sia più stretto della separazione).
2. Scegliere la finestra in base al dinamismo necessario: se un tono forte può nascondere uno debole di $-60\ \mathrm{dB}$ serve una finestra con lobi laterali sotto $-60\ \mathrm{dB}$ (Blackman, Kaiser).
3. Scegliere $N\ge L$ (potenza di due, zero-padding) per leggere meglio i picchi.
4. Correggere l'ampiezza con il guadagno coerente $\sum w[n]$ e ricordare la perdita di $\le3.9\ \mathrm{dB}$ se il tono non cade su un bin (tolta con lo zero-padding o con finestre a picco piatto).

Domande d'esame

1. Spiegare l'effetto della finestratura sullo spettro di una sinusoide e definire risoluzione e leakage. Traccia: $y=xw_R$, $Y=AW_R(\hat\omega-\hat\omega_0)$; nucleo di Dirichlet, zeri in $2\pi k/L$, lobo principale $4\pi/L$, lobi laterali con il primo a $-13.3\ \mathrm{dB}$ indipendente da $L$; due toni: interferenza, risoluzione $\sim2\pi/L$ legata a $L$, leakage che maschera toni deboli.
2. Perché la DFT di una sinusoide può dare un'ampiezza sbagliata? Come si corregge? Traccia: DFT = campioni della DTFT a $2\pi k/N$; solo se $\hat\omega_0=2\pi k_0/N$ (periodi interi) il vertice è campionato; altrimenti perdita fino a $3.92\ \mathrm{dB}$; zero-padding infittisce i campioni ma non migliora la risoluzione; esempio numerico $k_0=4.5$, $L=32$ (picco $20.38$ invece di $32$).
3. Confrontare finestra rettangolare e finestre rastremate. Traccia: compromesso lobo principale/lobi laterali; tabella con larghezze $4\pi/L,8\pi/L,12\pi/L$ e livelli $-13,-26,-31,-42,-58\ \mathrm{dB}$; esempio del tono debole a $-60\ \mathrm{dB}$ (Kaiser sì, Hamming no).