Perché le stime della densità spettrale di potenza da matplotlib.mlab.psd e scipy.signal.welch differiscono quando il numero di punti per finestra è pari?

Question

matplotlib.mlab.psd(...) e scipy.signal.welch(...) entrambi implementano il metodo di periodogramma medio di Welch per stimare la densità spettrale di potenza (PSD) di un segnale. Supponendo che i parametri equivalenti siano passati a ciascuna funzione, il PSD restituito dovrebbe essere equivalente.

Tuttavia, utilizzando una semplice funzione di test sinusoidale, vi sono differenze sistematiche tra i due metodi quando il numero di punti utilizzati per ogni finestra (la parola nperseg per scipy.signal.welch, la parola NFFT per mlab.psd) è anche, come mostrato di seguito per il caso di 64 punti per finestra

la trama superiore mostra il PSD calcolata mediante entrambi i metodi, mentre la trama inferiore visualizza il loro errore relativo (ossia la differenza assoluta dei due PSD divisa per la loro averag assoluta e). Un errore relativo più ampio è indicativo di un maggiore disaccordo tra i due metodi.

Le due funzioni hanno molto migliore accordo quando il numero di punti utilizzati per ogni finestra è dispari, come mostrato di seguito per lo stesso segnale ma elaborati con 65 punti per finestra

Aggiunta altre caratteristiche del segnale (ad es. rumore) diminuiscono un po 'questo effetto, ma è ancora presente, con errori relativi del ~ 10% tra i due metodi quando si utilizza un numero pari di punti per finestra. (Mi rendo conto che l'errore relativo alla frequenza più alta per i PSD calcolati con 65 punti per finestra sopra è relativamente grande, tuttavia si tratta della frequenza Nyquist del segnale, e non sono troppo preoccupato per le caratteristiche a frequenze così elevate. Sono più preoccupato dell'ampio e sistematico errore relativo nella maggior parte della larghezza di banda del segnale quando il numero di punti per finestra è pari).

C'è una ragione per questa discrepanza? Un metodo è preferibile all'altro in termini di precisione? Sto usando scipy versione 0.16.0 e matplotlib versione 1.4.3.

Il codice utilizzato per generare i grafici sopra è incluso di seguito. Per un segnale sinusoidale puro, impostare A_noise = 0; per un segnale rumoroso, impostare A_noise su un valore limitato.

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import welch 
from matplotlib import mlab 

# Sampling information 
Fs = 200. 
t0 = 0 
tf = 10 
t = np.arange(t0, tf, 1./Fs) 

# Pure sinusoidal signal 
f0 = 27. 
y = np.cos(2 * np.pi * f0 * t) 

# Add in some white noise 
A_noise = 1e-3 
y += (A_noise * np.random.randn(len(y))) 

nperseg = 64 # even 
# nperseg = 65 # odd 

if nperseg > len(y): 
    raise ValueError('nperseg > len(y); preventing unintentional zero padding') 

# Compute PSD with `scipy.signal.welch` 
f_welch, S_welch = welch(
    y, fs=Fs, nperseg=nperseg, noverlap=(nperseg // 2), 
    detrend=None, scaling='density', window='hanning') 

# Compute PSD with `matplotlib.mlab.psd`, using parameters that are 
# *equivalent* to those used in `scipy.signal.welch` above 
S_mlab, f_mlab = mlab.psd(
    y, Fs=Fs, NFFT=nperseg, noverlap=(nperseg // 2), 
    detrend=None, scale_by_freq=True, window=mlab.window_hanning) 

fig, axes = plt.subplots(2, 1, sharex=True) 

# Plot PSD computed via both methods 
axes[0].loglog(f_welch, S_welch, '-s') 
axes[0].loglog(f_mlab, S_mlab, '-^') 
axes[0].set_xlabel('f') 
axes[0].set_ylabel('PSD') 
axes[0].legend(['scipy.signal.welch', 'mlab.psd'], loc='upper left') 

# Plot relative error 
delta = np.abs(S_welch - S_mlab) 
avg = 0.5 * np.abs(S_welch + S_mlab) 
relative_error = delta/avg 
axes[1].loglog(f_welch, relative_error, 'k') 
axes[1].set_xlabel('f') 
axes[1].set_ylabel('relative error') 

plt.suptitle('nperseg = %i' % nperseg) 
plt.show() 

Answer 1

Mentre i parametri possono sembrare equivalenti, il parametro della finestra potrebbe differire leggermente per la finestra di dimensioni uguali. Più specificamente, salva una specifica finestra vettoriale, la finestra utilizzata dalla funzione di SciPy welch viene generato con

win = get_window(window, nperseg) 

che utilizza il parametro predefinito fftbins=True, e secondo scipy documentation:

Se Vero, creare un Finestra "periodica" pronta per l'uso con ifftift e moltiplicata per il risultato di un fft (VEDI ANCHE FFTFREQ).

Questo risultato in una finestra generata diversa per le lunghezze pari.Da this section of the Window function entry on Wikipedia, questo potrebbe offrire un leggero vantaggio in termini di prestazioni rispetto a window_hanning di Matplotlib che restituisce sempre la versione simmetrica.

Per utilizzare la stessa finestra è possibile specificare esplicitamente il vettore della finestra per entrambe le funzioni di stima PSD. Si potrebbe ad esempio calcolare questa finestra con:

win = scipy.signal.get_window('hanning',nperseg) 

Utilizzando questa finestra come parametro (con window=win in entrambe le funzioni) darebbe il seguente grafico in cui si può notare un accordo molto più vicino tra le due funzioni di stima PSD:

In alternativa, a patto che non si desidera che la proprietà della finestra periodica, si potrebbe anche usare:

win = mlab.window_hanning(np.ones(nperseg)) # or 
win = scipy.signal.get_window('hanning',nperseg,fftbins=False)