求支持N维数组的Welch周期图Python实现（类scipy.signal.welch）

现有可用工具与快速实现 #

• 基于scipy和numpy的轻量扩展：
  虽然scipy.signal.welch原生仅支持1D/2D信号，但可以通过遍历目标轴的方式扩展到N维。核心思路是对数组的非目标轴切片逐一调用welch，再拼接结果。代码示例：

  import numpy as np
  from scipy.signal import welch
  
  def n_dim_welch(x, fs=1.0, nperseg=None, noverlap=None, axis=-1):
      if x.ndim == 1:
          return welch(x, fs=fs, nperseg=nperseg, noverlap=noverlap)
      
      freqs = None
      psd_list = []
      # 构建索引器，遍历目标轴外的所有维度
      idx = [slice(None)] * x.ndim
      for i in range(x.shape[axis]):
          idx[axis] = i
          curr_slice = x[tuple(idx)]
          f, p = welch(curr_slice, fs=fs, nperseg=nperseg, noverlap=noverlap)
          if freqs is None:
              freqs = f
          psd_list.append(p)
      
      # 将结果堆叠回原维度结构
      psd = np.stack(psd_list, axis=axis)
      return freqs, psd
  

  这个函数的参数、返回值和scipy.signal.welch完全对齐，可直接替换使用。

• pycwt库的原生支持：
  pycwt（连续小波变换库）内置了支持多维信号的功率谱计算，其中包含与Welch算法逻辑一致的实现，无需手动处理轴迭代，直接传入N维数组即可完成计算。

无额外依赖的复刻方案 #

如果不想引入新库，用numpy手动实现N维Welch也不复杂，核心步骤：

1. 沿目标轴对数组分段，给每个分段加窗（比如np.hanning(nperseg)）
2. 对每个分段执行N维FFT（np.fft.fftn）
3. 计算每个分段的功率谱（取模平方后除以窗的能量和）
4. 对所有分段的功率谱做平均，得到最终周期图
5. 生成对应的频率轴（利用np.fft.fftfreq）

这种方式完全基于numpy，可以自定义窗函数、分段策略，适合需要高度定制的场景。

内容的提问来源于stack exchange，提问作者mocquin