matplotlib.mlab.psd() 函数中最大 PSD 频率箱的问题

在使用 matplotlib.mlab.psd() 函数时，开发人员发现返回的频率箱似乎不正确。

以下是测试代码：

import matplotlib.mlab as ml
import numpy as np

sampf = 500.
nfft = 2**4
testdat = np.random.randn(10000,)
p2, f2 = ml.psd(testdat, nfft, sampf, sides='twosided')
p1, f1 = ml.psd(testdat, nfft, sampf, sides='onesided')

print(testdat.shape)
print("Twosided")
print("\tbin1     : {:f} ".format(f2[0]))
print("\tbin2     : {:f} ".format(f2[1]))
print("\tbinlast  : {:f} ".format(f2[-1]))

print("onesided")
print("\tbin1     : {:f} ".format(f1[0]))
print("\tbin2     : {:f} ".format(f1[1]))
print("\tbinlast  : {:f} ".format(f1[-1]))

print("recreate")
f3 = np.arange(nfft) * (sampf / 2.) / nfft
print("\tbin1     : {:f} ".format(f3[0]))
print("\tbin2     : {:f} ".format(f3[1]))
print("\tbinlast  : {:f} ".format(f3[-1]))

输出结果如下：

Twosided
    bin1     : -250.000000
    bin2     : -218.750000
    binlast  : 218.750000
onesided
    bin1     : 0.000000
    bin2     : 31.250000
    binlast  : 250.000000
recreate
    bin1     : 0.000000
    bin2     : 15.625000
    binlast  : 234.375000

在双边模式下，最大频率 (binlast) 应该是采样频率的一半，但从输出结果可以看出，它与采样频率的一半相差一个箱位。
而在单边模式下，最大频率也比采样频率的一半小。

2、解决方案

根据一位答主提供的答案，这个问题与 matplotlib.mlab.psd() 函数的内部实现有关。

在单边模式下，matplotlib.mlab.psd() 函数不会返回负频率的箱位，因此最大频率只能达到采样频率的一半。

而在双边模式下，matplotlib.mlab.psd() 函数会返回负频率和正频率的箱位，但由于采样频率是偶数，因此最大频率无法达到采样频率的一半，而是会比采样频率的一半小一个箱位。

如果需要获得准确的频率箱位，可以使用以下方法：

import numpy as np

def psd(data, nfft=None, sampf=None, sides='default'):
  """
  Calculate the power spectral density (PSD) of a signal.

  Args:
    data: The input signal.
    nfft: The length of the FFT window.
    sampf: The sampling frequency.
    sides: The type of PSD to calculate. 'default' returns a two-sided PSD,
      'onesided' returns a one-sided PSD, and 'onlypositive' returns a PSD
      with only positive frequencies.

  Returns:
    A tuple containing the PSD and the corresponding frequencies.
  """

  # If no window is specified, use the entire signal.
  if nfft is None:
    nfft = len(data)

  # If no sampling frequency is specified, assume 1 Hz.
  if sampf is None:
    sampf = 1

  # Calculate the FFT of the signal.
  fft_data = np.fft.fft(data, nfft)

  # Calculate the power spectrum.
  psd = np.abs(fft_data)**2

  # If a one-sided PSD is requested, only return the positive frequencies.
  if sides == 'onesided':
    psd = psd[:nfft // 2 + 1]

  # If a positive-only PSD is requested, only return the positive frequencies.
  elif sides == 'onlypositive':
    psd = psd[:nfft // 2]

  # Calculate the corresponding frequencies.
  freqs = np.fft.fftfreq(nfft, d=1 / sampf)

  # If a one-sided PSD is requested, only return the positive frequencies.
  if sides == 'onesided':
    freqs = freqs[:nfft // 2 + 1]

  # If a positive-only PSD is requested, only return the positive frequencies.
  elif sides == 'onlypositive':
    freqs = freqs[:nfft // 2]

  return psd, freqs

通过这个函数可以准确地计算出PSD和相应的频率。