Signal.welch 函数：f 和 pxx 结果为 0

Question

我有一个长度为 240 的 BOLD 活动数组。每 2 秒进行一次记录（您可以下载部分数据 这里（该文件可以用.txt编辑器正常打开）。我想要用 signal.welch 函数分析这个时间序列，

f, pxx = signal.welch(data, fs=0.5, window='hanning', nperseg=50, noverlap=25, scaling='density', average='mean')

它给了我以下错误

。

ValueError：noverlap 必须小于 nperseg。

当我设置 noverlap = None 时，没有显示错误，但 f 等于 0 并且 pxx 是一个 0 数组。

非常感谢您的建议！！

Answer 1

例如，如果我读取数据文件，

data = np.loadtxt('rest1_LeftInsula.1D')

data 将是一个形状为 (240,) 的一维 numpy 数组，并且调用

f, pxx = signal.welch(data, fs=0.5, window='hanning', nperseg=50, noverlap=25, scaling='density', average='mean')

不会出现错误。

如果我将 data 重塑为形状为 (240, 1) 的数组，并将其传递给 welch （这可以通过将其索引为 data[: , None] 或使用 reshape 方法 data.reshape((240, 1)))，我得到与您报告的相同的错误：

In [11]: f, pxx = signal.welch(data[:,None], fs=0.5, window='hanning', nperseg=5, 0, noverlap=25, scaling='density', average='mean')
[...]/scipy/signal/spectral.py:1964: UserWarning: nperseg = 50 is greater than input length  = 1, using nperseg = 1
  warnings.warn('nperseg = {0:d} is greater than input length '
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-11-0e5235350cfd> in <module>
----> 1 f, pxx = signal.welch(data[:,None], fs=0.5, window='hanning', nperseg=50, noverlap=25, scaling='density', average='mean')

[...]/scipy/signal/spectral.py in welch(x, fs, window, nperseg, noverlap, nfft, detrend, return_onesided, scaling, axis, average)
    446 
    447     """
--> 448     freqs, Pxx = csd(x, x, fs=fs, window=window, nperseg=nperseg,
    449                      noverlap=noverlap, nfft=nfft, detrend=detrend,
    450                      return_onesided=return_onesided, scaling=scaling,

[...]/scipy/signal/spectral.py in csd(x, y, fs, window, nperseg, noverlap, nfft, detrend, return_onesided, scaling, axis, average)
    580 
    581     """
--> 582     freqs, _, Pxy = _spectral_helper(x, y, fs, window, nperseg, noverlap, nfft,
    583                                      detrend, return_onesided, scaling, axis,
    584                                      mode='psd')

[...]/scipy/signal/spectral.py in _spectral_helper(x, y, fs, window, nperseg, noverlap, nfft, detrend, return_onesided, scaling, axis, mode, boundary, padded)
   1756         noverlap = int(noverlap)
   1757     if noverlap >= nperseg:
-> 1758         raise ValueError('noverlap must be less than nperseg.')
   1759     nstep = nperseg - noverlap
   1760 

ValueError: noverlap must be less than nperseg.

问题是，默认情况下，welch 沿着输入数组的最后轴。如果该数组的形状为 (240, 1)，welch 会尝试将计算应用于二维数组的每个“行”。但每行的长度为 1，对于给定的 nperseg 和 noverlap 值来说太小，这会导致（有点神秘）错误。

您还没有展示如何读取文件并创建 data，但我怀疑它正在创建形状为 (240, 1) 的数组（或其他一些数据结构，例如 Pandas DataFrame） 。

要解决此问题，您可以将数据“展平”为一维数组（例如将 data.ravel() 传递给 welch），或传递参数 axis=0 到 welch 告诉它沿着第一个维度而不是最后一个维度行动。如果您选择后者，请注意 pxx 的形状将为 (26, 1)，而不是当 data 具有形状时获得的形状 (26,) （240，）。