python 2.7 - Why do the power spectral density estimates from matplotlib.mlab.psd and scipy.signal.welch differ when the number of points per window is even? -

matplotlib.mlab.psd(...) , scipy.signal.welch(...) both implement welch's average periodogram method estimate power spectral density (psd) of signal. assuming equivalent parameters passed each function, returned psd should equivalent.

however, using simple sinusoidal test function, there systematic differences between 2 methods when number of points used per window (the nperseg keyword scipy.signal.welch; nfft keyword mlab.psd) even, shown below case of 64 points per window

the top plot shows psd computed via both methods, while bottom plot displays relative error (i.e. absolute difference of 2 psds divided absolute average). larger relative error indicative of larger disagreement between 2 methods.

the 2 functions have better agreement when number of points used per window odd, shown below same signal processed 65 points per window

adding other features signal (e.g. noise) diminishes effect, still present, relative errors of ~10% between 2 methods when number of points used per window. (i realize relative error @ highest frequency psds computed 65 points per window above relatively large. however, @ signal's nyquist frequency, , i'm not concerned features @ such high frequencies. i'm more concerned large , systematic relative error across majority of signal's bandwidth when number of points per window even).

is there reason discrepancy? 1 method preferable other in terms of accuracy? using scipy version 0.16.0 , matplotlib version 1.4.3.

the code used generate above plots included below. pure sinusoidal signal, set a_noise = 0; noisy signal, set a_noise finite value.

import numpy np import matplotlib.pyplot plt scipy.signal import welch 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 white noise a_noise = 1e-3 y += (a_noise * np.random.randn(len(y)))  nperseg = 64    # # nperseg = 65  # odd  if nperseg > len(y):     raise valueerror('nperseg > len(y); preventing unintentional 0 padding')  # compute psd `scipy.signal.welch` f_welch, s_welch = welch(     y, fs=fs, nperseg=nperseg, noverlap=(nperseg // 2),     detrend=none, scaling='density', window='hanning')  # compute psd `matplotlib.mlab.psd`, using parameters # *equivalent* 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() 

while parameters may appear equivalent, window parameter may differ sized window. more specifically, unless provided specific window vector, window used scipy's welch function generated with

win = get_window(window, nperseg) 

which uses default parameter fftbins=true, , according scipy documentation:

if true, create “periodic” window ready use ifftshift , multiplied result of fft (see fftfreq).

this result in different generated window lengths. this section of window function entry on wikipedia, give slight performance advantage on matplotlib's window_hanning returns symmetric version.

to use same window can explicitly specify window vector both psd estimation functions. example compute window with:

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

using window parameter (with window=win in both functions) give following plot may notice closer agreement between 2 psd estimation functions:

alternatively, assuming not want periodic window property, use:

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