If we statescale_by_freq=False, then we are asking for the power spectrum of the signal, instead of the power spectral density; i.e. units of V^2 instead of V^2/Hz, for example.

mlab.psd, which actually ismlab._spectral_helper under the hood, does not return the expected value of the power in a sinusoid. Furthermore, in contrast to the PSD case, it does not correct for the energy of the signal lost to the window function, meaning that the results differ depending on the window used.

Here's an illustration. If we have a simple signal of a cosine oscillation happening at 100Hz with amplitude of 10, we would expect the power spectrum result at 100Hz to be 50, regardless of the window used. (A cosine like this that is identically zero at each integer second means we can use one second FFTs and compare the rectangular window to the hanning window, without worrying about any distortion)

In MATLAB:

fs=2048;tt=[0:1/fs:60-1/fs];x = 10*cos(2*pi*100*tt);[pR, f]=pwelch(x, ones(1,fs), 0, fs, fs, 'power');[pH, f]=pwelch(x, hann(fs), 0, fs, fs, 'power');

givespR(101)==50 andpH(101)==50.000..., as expected.

Inmatplotlib.mlab:

fs = 2048tt = np.arange(fs*64)/float(fs)x = 10*np.cos(2*np.pi*100*tt)# One second FFT segmentspsdHann, ff = mlab.psd(x, Fs=fs, NFFT=fs, window=mlab.window_hanning, scale_by_freq=False)psdRect, ff = mlab.psd(x, Fs=fs, NFFT=fs, window=mlab.window_none, scale_by_freq=False)

givespsdRect[100] == 102400, which too large by a factor ofnfft andpsdHann[100] == 68233.3, which differs further by some window correction factor.

I came across this inscipy/scipy#4682, where I was adapting the code for inclusion in scipy's signal processing repertoire. There, I usedscale = 1.0 / win.sum()**2 for the power spectrum case. I haven't looked for the exact places where_spectral_helper should be fixed, at this time.