Trasformata di Fourier a tempo breve e spettrogramma

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 →) dà lo spettro globale di un segnale finito: dice quali frequenze ci sono ma non quando. Un segnale audio, una voce, un fischio che si sposta, hanno un contenuto in frequenza che cambia nel tempo. Calcolare una DFT enorme su tutto il segnale richiederebbe di aspettare tutti i campioni e di fare un calcolo gigantesco, inaccettabile per l'elaborazione in tempo reale. La soluzione è vedere il segnale lungo come la concatenazione di tanti segnali brevi e calcolare molte DFT brevi: la trasformata di Fourier a tempo breve (short-time Fourier transform, STFT, o short-time DFT, STDFT; nelle slide anche time-dependent DFT), la cui rappresentazione grafica è lo spettrogramma (spectrogram).

Definizione

Definizione (STFT, "DFT dipendente dal tempo"). Sia $x[n]$ una sequenza lunga e $w[n]$ una finestra di analisi (analysis window, per esempio di Hann) diversa da zero solo per $n=0,\dots,L-1$, con $L$ molto minore della lunghezza di $x$. Per ogni istante di analisi $\ell$ si definisce $$X[k,\ell]=\sum_{n=0}^{L-1}w[n]\,x[\ell+n]\,e^{-j2\pi kn/N},\qquad k=0,1,\dots,N-1,\quad N\ge L.$$

Contiene due passi:

1. Finestratura: $w[n]x[\ell+n]$ preleva dal segnale il tratto lungo $L$ che comincia in $\ell$ (frame, segmento) e lo pesa. Cambiando $\ell$ si sposta la finestra sul segnale.
2. DFT breve: la DFT a $N$ punti di quel frame ($N>L$ con zero-padding se si vuole una griglia di frequenze più fitta, 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 →) dà lo spettro locale.

Il parametro $k$ è la frequenza: $\hat\omega_k=2\pi k/N$, cioè $f_k=kF_s/N$. Il parametro $\ell$ è il tempo. L'istante $\ell$ si muove a salti di $R$ campioni (hop size), $\ell=mR$ con $m=0,1,2,\dots$, con $1\le R<L$. Se $R<L$ i frame si sovrappongono e ogni campione cade in almeno una finestra; la scelta comune è $R=L/2$ (sovrapposizione del $50\%$). Il numero di frame per un segnale di lunghezza $D$ è $\lfloor(D-L)/R\rfloor+1$.

Esempio. Un segnale di $D=10\,000$ campioni, finestra $L=301$, passo $R=30$ ($\approx0.1L$) e $N=1024$: $\lfloor(10000-301)/30\rfloor+1=324$ frame; per un segnale reale bastano i bin $k=0,\dots,512$ ($N/2+1=513$ valori) per la simmetria coniugata. Con $L=75$ e $R=8$: $1241$ frame. Il numero di valori calcolati è $324\times513$, circa $1.7\cdot10^5$: ogni frame costa una FFT a $N=1024$ punti (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 →).

Lo spettrogramma

$X[k,\ell]$ è una funzione di due variabili: per ogni istante c'è un diverso spettro locale. Il modo di rappresentarla è lo spettrogramma: un'immagine in cui l'asse orizzontale è il tempo ($\ell=mR$), l'asse verticale è la frequenza ($\hat\omega_k$ in rad/campione, $\hat\omega/\pi\in[0,1]$, oppure Hz) e il livello di grigio (o di colore) è proporzionale a $\lvert X[k,\ell]\rvert$ oppure a $20\log_{10}\lvert X[k,\ell]\rvert$ (scala in dB, che mostra anche le componenti deboli) o a $\lvert X\rvert^2$ (densità di potenza). Un'orizzontale chiara indica un tono stabile; una linea inclinata una frequenza che cresce o decresce; un trattino verticale un impulso (che ha spettro piatto).

Esempio (come nel laboratorio 2). scipy.signal.spectrogram con finestra di Hamming di $256$ campioni, nfft=256 e noverlap=128: $R=128$, risoluzione in frequenza $F_s/256$ e per un segnale di $10\,000$ campioni si hanno $\lfloor(10000-256)/128\rfloor+1=77$ frame e $129$ bin; l'asse delle frequenze si normalizza dividendo per $F_s/2$ per leggere $\hat\omega/\pi$. I toni si vedono come righe orizzontali, da cui si ricavano le $\hat\omega_0$ per i filtri notch (Filtri notch e applicazioni dei filtri FIRUno zero di H(z) sulla circonferenza unitaria in z0 = e^{j w0} annulla (a regime) le sinusoidi alla pulsazione w0. Per eliminare un coseno servono i due zeri coniugati: H(z) = (1 - e^{jw0} z^-1)(1 - e^{-jw0} z^-1) = 1 - 2 cos(w0) z^-1 + z^-2, cioè h = {1, -2cos w0, 1}, un FIR simmetrico di tipo I (ritardo 1). Si normalizza dividendo per 2 - 2cos(w0) per avere guadagno 1 in continua; più toni si eliminano mettendo in cascata (convoluzione) un filtro per ogni tono. Segue una panoramica delle applicazioni dei FIR (equalizzazione audio), vantaggi (fase lineare, stabilità) e costo (N moltiplicazioni per campione, molte più di un IIR).Filtri notch e applicazioni dei filtri FIR →).

Il compromesso tra tempo e frequenza

La scelta più importante è la lunghezza $L$ della finestra, perché agisce in due modi opposti.

Risoluzione in frequenza (frequency resolution). Ogni riga spettrale compare nello spettrogramma come il lobo principale della trasformata della finestra, la cui larghezza è inversamente proporzionale a $L$: $\Delta\hat\omega\simeq c\,\dfrac{2\pi}{L}$, con $c$ costante che dipende dalla finestra (lobo principale $4\pi/L$ rettangolare, $8\pi/L$ per Hann e Hamming, $12\pi/L$ Blackman). Due frequenze che distano meno di circa la mezza larghezza del lobo si fondono. Per separare righe vicine serve una finestra lunga: serve tempo per vedere che le due sinusoidi hanno periodi diversi. Una finestra più corta di un periodo del segnale non dà alcuna informazione sulla periodicità.

Risoluzione nel tempo (time resolution). Il frame $[\ell,\ell+L-1]$ vede un tratto lungo $L$: un cambiamento improvviso del segnale (per esempio una frequenza che parte di colpo) comparirà sfocato su un intervallo di durata confrontabile con $L$. Per localizzare con precisione i cambiamenti serve una finestra corta.

Compromesso (time-frequency trade-off). La finestra deve essere corta per seguire i cambiamenti nel tempo e lunga per risolvere frequenze vicine: i due requisiti sono opposti e non si possono soddisfare insieme. Per un segnale sconosciuto si calcolano spettrogrammi con lunghezze diverse: ciò che è evidente in uno può essere nascosto in un altro.

Esempio: due toni vicini e una transizione

Il segnale di prova delle slide (lunghezza $10\,000$) è $$x[n]=\begin{cases}\cos(0.2\pi n)+1.5\cos(0.227\pi n)&0\le n<5000\\ \cos(0.2\pi n)&5000\le n<7000\\ \cos(0.6\pi n)+\cos(0.7\pi n)&7000\le n<10000.\end{cases}$$ Si calcola la STFT con finestra di Hann, $N=1024$, in tre istanti (finestra centrata nell'istante indicato; valori calcolati con Python).

I toni a $0.2\pi$ e $0.227\pi$ distano $0.027\pi$. Con $L=301$ la mezza larghezza del lobo di Hann è $4\pi/301=0.0133\pi$, quindi li risolve; con $L=75$ vale $4\pi/75=0.053\pi$, quindi li fonde. Quelli a $0.6\pi$ e $0.7\pi$ distano $0.1\pi$ e si separano anche con $L=75$.

Il costo nel tempo. Nell'istante $n=7000$ il tono a $0.6\pi$ compare di colpo. Misurando l'ampiezza della riga a $0.6\pi$ nei frame, la salita dal $10\%$ al $90\%$ del massimo avviene in $\approx145$ campioni per $L=301$ e in $\approx38$ per $L=75$: la finestra lunga sfoca la transizione di circa $L/2$ campioni, quella corta la localizza bene. Un modello: per una finestra di Hann centrata in $t$, con $u=(t-7000)/L$ limitato in $[-\tfrac12,\tfrac12]$, la riga ha ampiezza relativa $\tfrac12+u+\tfrac{\sin2\pi u}{2\pi}$ (è l'area della finestra che cade dopo l'istante $7000$); in $t=7000$ vale $0.5$ (con Python $0.51$).

Questa relazione è un esempio del principio generale: finestra lunga $\Rightarrow$ risoluzione in frequenza $\propto1/L$ e sfocatura temporale $\propto L$; il prodotto delle due incertezze resta costante.

Come scegliere i parametri

• Finestra $L$: dalla risoluzione voluta, $L\ge c\cdot2\pi/\Delta\hat\omega$ ($c\simeq2$ per Hann/Hamming) e dalla durata dei fenomeni da seguire ($L$ non più lunga dell'evento).
• Forma di $w$: Hann o Hamming per un buon equilibrio; rettangolare solo se serve il lobo più stretto e il leakage non conta (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 →).
• Passo $R$: $R=L/2$ come predefinito; $R\approx0.1L$ per un'immagine più liscia (costa più frame).
• $N$: potenza di due, $N\ge L$; zero-padding per tracciare meglio i picchi (non cambia la risoluzione).

Il codice Python è in Esercizio - Spettrogramma e finestra di analisi.

Domande d'esame

1. Definire la STFT e lo spettrogramma. Perché si usano tante DFT brevi invece di una sola lunga? Traccia: $X[k,\ell]=\sum w[n]x[\ell+n]e^{-j2\pi kn/N}$, finestra, passo $R<L$ (tipicamente $L/2$), frame; segnale non stazionario; calcolo di una DFT enorme e impossibilità del tempo reale; rappresentazione in dB come immagine tempo-frequenza; numero di frame.
2. Spiegare il compromesso tra risoluzione in tempo e in frequenza nello spettrogramma. Traccia: $\Delta\hat\omega\simeq c\,2\pi/L$ (lobo della finestra); due toni più vicini della mezza larghezza si fondono; finestra corta per transizioni brevi (sfocatura $\sim L$); esempio con $0.2\pi$ e $0.227\pi$ ($L=301$ risolve, $L=75$ no); strategia di più lunghezze.
3. Come si sceglie la lunghezza della finestra per risolvere due toni a distanza $\Delta\hat\omega$? Traccia: mezza larghezza del lobo principale $<\Delta\hat\omega$: $L>8\pi/(2\Delta\hat\omega)$ per Hann (con $\Delta\hat\omega=0.027\pi$ si ricava $L>148$).