Esercizio - Laboratorio 3, DFT e analisi spettrale

Testo (laboratorio 3 del corso Multimedia Signal Processing, UniPD; lezione 19). Il laboratorio studia i compromessi tra risoluzione in frequenza e leakage.

Parte 1. Il segnale $x(t)=\cos(2\pi f_1t)+\cos(2\pi f_2t)+\cos(2\pi f_3t)$ con $f_1=4$, $f_2=5$, $f_3=6\ \mathrm{kHz}$ è osservato per $5\ \mathrm{ms}$ con $F_s=20\ \mathrm{kHz}$.

• (a) Plottare $x[n]$ per $0\le n\le L-1$, $L=100$.
• (b) Calcolare la DFT a $256$ punti e plottarne il modulo per $0\le f/F_s\le1$.
• (c) Calcolare l'IDFT a $256$ punti e verificare il recupero del segnale (più gli zeri di padding).
• (d) Finestrare con una finestra di Hamming di lunghezza $L$ e plottare.
• (e) Ripetere (b) per il segnale finestrato; che cosa si osserva?
• (f) Ripetere (b) ed (e) per $L=50,25,10$: effetti della lunghezza.

Parte 2. $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 è $60\ \mathrm{dB}$ più debole degli altri. (g) Confrontare finestra rettangolare, di Hamming e di Kaiser per rivelarlo, con $A_{dB}=60+10=70\ \mathrm{dB}$ (margine di $10\ \mathrm{dB}$), larghezza di transizione $\Delta f=0.2\ \mathrm{kHz}$, $\beta$ e lunghezza della finestra di Kaiser dalle formule $$\beta=0.1102(A_{dB}-8.7)\ (A_{dB}>50),\qquad L\ge\frac{A_{dB}-8}{14.36\,\Delta f}F_s+1,$$ con la DTFT calcolata in $2000$ frequenze equispaziate in $[0,5]\ \mathrm{kHz}$, normalizzata a $0\ \mathrm{dB}$ nel massimo.

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

Il codice

import numpy as np
from scipy.signal import find_peaks

# ---- parti a-f
fs = 20e3
f1, f2, f3 = 4e3, 5e3, 6e3
def segnale(L):
    t = np.arange(L) / fs
    return np.cos(2*np.pi*f1*t) + np.cos(2*np.pi*f2*t) + np.cos(2*np.pi*f3*t)

L = 100
x = segnale(L)                       # a) 100 campioni = 5 ms
X = np.fft.fft(x, 256)               # b) DFT a 256 punti (zero-padding automatico)
f_norm = np.arange(256) / 256        #    asse f/fs; plt.plot(f_norm, np.abs(X))
xr = np.fft.ifft(X)                  # c) x seguito da 156 zeri
xw = x * np.hamming(L)               # d) finestra di Hamming
Xw = np.fft.fft(xw, 256)             # e)

def picchi(L, finestra):             # f) quante righe si distinguono
    x = segnale(L) * finestra(L)
    M = np.abs(np.fft.fft(x, 256))[:129]
    p, _ = find_peaks(M, height=0.2 * M.max())
    return np.round(np.arange(129)[p] * fs / 256 / 1e3, 2)

# ---- parte g
fs = 10e3
f1, f2, f3 = 1.5e3, 2.5e3, 3.5e3
A_dB, df = 70, 0.2e3
beta = 0.1102 * (A_dB - 8.7)
L = int(np.ceil((A_dB - 8) * fs / (14.36 * df) + 1))
n = np.arange(L)
x = np.cos(2*np.pi*f1/fs*n) + 1e-3*np.cos(2*np.pi*f2/fs*n) + np.cos(2*np.pi*f3/fs*n)
f = np.linspace(0, 5e3, 2000)
E = np.exp(-2j*np.pi*np.outer(f/fs, n))      # DTFT in 2000 frequenze
for nome, w in (("rettangolare", np.ones(L)), ("Hamming", np.hamming(L)), ("Kaiser", np.kaiser(L, beta))):
    Xd = np.abs(E @ (x * w))
    dB = 20*np.log10(Xd / Xd.max())          # plt.plot(f, dB)

Risultati (eseguendo il codice) e commenti

Parte 1.

• (a) Con $F_s=20\ \mathrm{kHz}$, $5\ \mathrm{ms}$ danno $100$ campioni; le pulsazioni sono $0.4\pi$, $0.5\pi$, $0.6\pi$ (somma di tre coseni con ampiezza massima $3$ nell'origine).
• (b) I bin sono distanti $F_s/256=78.125\ \mathrm{Hz}$; i toni cadono alle posizioni $51.2$, $64$, $76.8$ (frazionarie tranne la centrale). I tre toni hanno $20$, $25$ e $30$ periodi interi in $100$ campioni: una DFT a $100$ punti darebbe tre righe pure di altezza $L/2=50$. Con $N=256$ i bin non coincidono con i vertici né con gli zeri della finestra: il modulo mostra tre picchi di altezza $50.6$, $50.0$, $50.6$ (vicini a $L/2$, alterati dalla somma dei lobi dei toni vicini) e, attorno, i lobi laterali della finestra rettangolare, che a $-13\ \mathrm{dB}$ danno picchi spuri a $3.67$ e $6.33\ \mathrm{kHz}$ di altezza $12.1$.
• (c) L'IDFT restituisce $x$ nei primi $100$ campioni con errore $6\cdot10^{-16}$ e zeri (errore $6\cdot10^{-16}$) nei successivi $156$: la DFT con zero-padding è invertibile.
• (d)-(e) Con la finestra di Hamming i lobi laterali spariscono: restano tre picchi puliti a $3.98$, $5.00$, $6.02\ \mathrm{kHz}$ (altezza $\approx26.7$, circa metà del caso rettangolare perché la finestra ha somma $\approx0.54L$); il prezzo è un lobo principale più largo (circa il doppio).
• (f) Effetto della lunghezza (picchi trovati con soglia $20\%$ del massimo):

Più $L$ è corto, più i lobi principali si allargano ($\propto1/L$, risoluzione $F_s/L=200,\ 400,\ 800,\ 2000\ \mathrm{Hz}$ per la rettangolare): la distanza tra i toni è $1\ \mathrm{kHz}$, quindi la rettangolare li separa fino a $L=25$ ($F_s/L=800\ \mathrm{Hz}<1\ \mathrm{kHz}$) ma non per $L=10$; l'Hamming (lobo largo il doppio) li separa bene fino a $L=50$ e a $L=25$ già li sposta e li fonde parzialmente. Con $L=10$ restano solo picchi ampi (rettangolare: due picchi spuri a $4.77$ e $5.23\ \mathrm{kHz}$ per interferenza dei lobi; Hamming: uno solo a $5\ \mathrm{kHz}$). Con la rettangolare, inoltre, i picchi spuri dei lobi laterali compaiono già a $L=100$.

Parte 2 (g). Con $A_{dB}=70$ e $\Delta f=0.2\ \mathrm{kHz}$: $\beta=0.1102\cdot61.3=6.755$ e $L\ge\frac{62}{14.36\cdot0.2}\cdot10+1=216.9$, quindi $L=217$. Livelli normalizzati nell'intorno di $2.5\ \mathrm{kHz}$:

Il tono debole sta a $-60\ \mathrm{dB}$: emerge solo se i lobi laterali dei due toni forti stanno sotto quel livello. Rettangolare ($-13\ \mathrm{dB}$) e Hamming ($-42\ \mathrm{dB}$) lasciano leakage sopra $-60\ \mathrm{dB}$ attorno a $2.5\ \mathrm{kHz}$ e lo coprono; la finestra di Kaiser, progettata per attenuare i lobi laterali di $70\ \mathrm{dB}$, porta il leakage a $-67.4\ \mathrm{dB}$ e il picco del tono, a $-60.9\ \mathrm{dB}$, è chiaramente visibile (valori calcolati con Python).