Esercizio - Spettrogramma e finestra di analisi

Testo (lezione 18 e laboratorio 2, corso Multimedia Signal Processing, UniPD). Si considera il segnale di prova delle slide, di $10\,000$ campioni: $$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}$$

1. Scrivere una funzione che calcola la STFT con finestra di Hann di lunghezza $L$, passo $R$ e DFT a $N$ punti, e calcolare lo spettrogramma con $L=301$ ($R=30$) e con $L=75$ ($R=8$), $N=1024$.
2. Dire quanti frame si ottengono e quali frequenze sono visibili nei frame centrati in $n=2000$, $6000$, $8000$ nei due casi.
3. Spiegare perché con $L=75$ i due toni vicini non si separano e con $L=301$ sì, e che cosa si perde con $L=301$ nel tempo (transizione in $n=7000$).
4. Scegliere $L$ per separare due toni distanti $0.027\pi$.

Teoria usata: 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 →, 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 →, 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 →.

Il codice

import numpy as np
from scipy.signal import find_peaks

n = np.arange(10000)
x = np.where(n < 5000, np.cos(0.2*np.pi*n) + 1.5*np.cos(0.227*np.pi*n),
    np.where(n < 7000, np.cos(0.2*np.pi*n), np.cos(0.6*np.pi*n) + np.cos(0.7*np.pi*n)))

def stft(x, w, R, N):
    L = len(w)
    nfr = (len(x) - L) // R + 1                                  # numero di frame
    return np.array([np.fft.rfft(x[m*R:m*R + L] * w, N) for m in range(nfr)])   # frame x bin

def picchi(riga, N):                                              # picchi sopra il 10% del massimo, in unita' di pi
    p, _ = find_peaks(riga, height=0.1 * riga.max())
    return np.round(p / N * 2, 3), np.round(riga[p], 1)

N = 1024
for L, R in ((301, 30), (75, 8)):
    X = stft(x, np.hanning(L), R, N)                              # spettrogramma = 20*log10(abs(X).T)
    print("L =", L, " frame:", X.shape[0], " bin:", X.shape[1])
    for centro in (2000, 6000, 8000):
        m = (centro - L // 2) // R
        print("  n =", centro, picchi(np.abs(X[m]), N))

Risultati (eseguendo il codice)

Numero di frame. $\lfloor(10000-301)/30\rfloor+1=324$ per $L=301$, $R=30$; $\lfloor(10000-75)/8\rfloor+1=1241$ per $L=75$, $R=8$. In tutti e due i casi $N/2+1=513$ bin (segnale reale).

Frame centrati negli istanti indicati (frequenze in unità di $\pi$, altezze tra parentesi):

Perché. La finestra di Hann ha lobo principale largo $8\pi/L$ e mezza larghezza $4\pi/L$: $0.0133\pi$ per $L=301$, $0.053\pi$ per $L=75$. I due toni distano $0.027\pi$: $0.027\pi>0.0133\pi$, quindi con $L=301$ i lobi sono separati; $0.027\pi<0.053\pi$, quindi con $L=75$ si fondono in un solo picco a $0.217\pi$ (tra i due toni, più vicino al più forte). I toni a $0.6\pi$ e $0.7\pi$ distano $0.1\pi>0.053\pi$ e si separano anche con $L=75$, con picchi più larghi. Con $L=75$ le altezze sono più basse perché la somma dei pesi della finestra è minore ($\sum w\approx L/2$).

Costo nel tempo. Nell'istante $n=7000$ compare di colpo il tono a $0.6\pi$. L'ampiezza della riga a $0.6\pi$ nei frame sale dal $10\%$ al $90\%$ del massimo in $\approx145$ campioni per $L=301$ e in $\approx38$ per $L=75$: la finestra lunga sfoca la transizione su circa $L/2$ campioni; quella corta la localizza bene ma non separa i toni vicini. Non c'è una scelta di $L$ che faccia bene entrambe le cose: da qui l'abitudine di calcolare spettrogrammi con $L$ diversi.

Scelta di $L$. Si impone che la mezza larghezza del lobo di Hann, $4\pi/L$, sia minore della separazione $\Delta\hat\omega=0.027\pi$: $L>4\pi/(0.027\pi)=148.1$, cioè almeno $L=149$. Nella nostra prova $L=301$ (quasi il doppio) separa bene i toni; $L=75$ no, come previsto.

Nota di laboratorio. Nel laboratorio 2 lo spettrogramma si ottiene con scipy.signal.spectrogram(signal, fs, window=np.hamming(256), nfft=256, noverlap=128): finestra di Hamming di $256$ campioni, passo $R=128$ (sovrapposizione del $50\%$), risoluzione $F_s/256$; l'asse delle frequenze si divide per $F_s/2$ per leggere $\hat\omega/\pi$. Con un segnale di $10\,000$ campioni darebbe $77$ frame e $129$ bin.