Bug summary

If my understanding is correct:

• NFFT of mlab._spectral_helper corresponds to nperseg of spectral funtions of scipy.signal.
  It's the actual length of signal (segment) for DFT.
• pad_to of mlab._spectral_helper corresponds to nfft of spectral funtions of scipy.signal.
  It's the length of DFT, i.e. parameter n of fft funtion.
  

  matplotlib/lib/matplotlib/mlab.py

  Line 385 in 3418bad

  result = np.fft.fft(result, n=pad_to, axis=0)[:numFreqs, :]

But the implementation of Matplotlib is really confusion:

• If length of signal is less than NFFT, signal will pad to length of NFFT.
  

  matplotlib/lib/matplotlib/mlab.py

  Lines 342 to 351 in 3418bad

  	# zero pad x and y up to NFFT if they are shorter than NFFT
  	if len(x) < NFFT:
  	n = len(x)
  	x = np.resize(x, NFFT)
  	x[n:] = 0
  	if not same_data and len(y) < NFFT:
  	n = len(y)
  	y = np.resize(y, NFFT)
  	y[n:] = 0

  And the window is determined by NFFT instead of actual length of signal. This would cause wrong window correction.
  

  matplotlib/lib/matplotlib/mlab.py

  Lines 376 to 380 in 3418bad

  	if not np.iterable(window):
  	window = window(np.ones(NFFT, x.dtype))
  	if len(window) != NFFT:
  	raise ValueError(
  	"The window length must match the data's first dimension")

• The Nyquist frequency is NFFT/2 instead of pad_to/2.
  

  matplotlib/lib/matplotlib/mlab.py

  Lines 411 to 416 in 3418bad

  	# if we have a even number of frequencies, don't scale NFFT/2
  	if not NFFT % 2:
  	slc = slice(1, -1, None)
  	# if we have an odd number, just don't scale DC
  	else:
  	slc = slice(1, None, None)

Code for reproduction

import numpy as np
from matplotlib import mlab
fs = 1000
t = np.arange(0, 1, 1/fs)
rng = np.random.default_rng(42)
s = np.sin(2 * np.pi * 100 * t) + 0.5 * rng.normal(size=t.shape)

nfft = 2048
nsignal = len(s) #1000
psd, freqs = mlab.psd(s, NFFT=nfft, Fs=fs, window=mlab.window_none)

mag = np.abs(np.fft.rfft(s, n=nfft))
psd1 = (mag)**2/nsignal/fs
if not nfft % 2:
  psd1[1:-1] *= 2
else: 
  psd1[1:] *= 2
relative_error = np.abs( 2 * (psd-psd1)/(psd + psd1) )
print(relative_error.max())

Actual outcome

0.6876640419947717

Expected outcome

0