先交代一个背景我在做某嵌入式平台的频谱监测模块时目标环境里装不了 scipy只有 numpy 和一套精简的数学库。可业务那边验收脚本偏偏是用 scipy.signal.welch 的标准输出校的PSD 曲线差一点都对不上。没办法我只好把 welch 从里到外捋了一遍用 python 手写了一个纯 numpy 版本然后逐点对拍最终把最大误差压到了 1e-12 这个量级。整个过程走下来最磨人的反而不是 FFT而是分段、加窗、归一化、单边谱补偿这些“一眼看上去很简单”的细节。这篇文章就把整个复现过程、参数换算逻辑和踩过的坑完整写出来。你如果也想自己实现功率谱密度估计或者需要把 scipy 的行为搬到受限环境里可以直接参考这个版本省去翻源码和试错的功夫。1. 为什么要复现一个已经封装好的 welch1.1 一行函数背后的复杂度很多人第一次用 scipy.signal.welch 的时候会觉得这东西太“黑盒”了给一段信号传个 fs出来频率轴和 PSD完事。但真正到了需要复现的时候才发现这一行调用背后至少藏了六件事检查输入信号维度和长度决定是否允许二维数组按轴计算根据 nperseg 和 noverlap 计算分段索引把长信号切成一堆重叠片段对每一段做去趋势detrend默认是减掉均值生成窗函数并让每一段乘上窗对加窗后的每一段做 FFT然后按 density 或 spectrum 模式归一化处理单边谱的 2 倍补偿最后把各段结果平均起来。任何一个环节的顺序错了或者归一化因子差个 fs输出就会整体偏移。更隐蔽的是窗函数的对称模式、奇数段长度、noverlap 等于 0 这些边界情况scipy 都做了特殊处理手动实现时一个没注意对拍误差就会从 1e-15 直接跳到 1e-2。1.2 哪些场景必须自己写我这次是因为目标环境没有 scipy属于典型的平台限制。但除了这个还有几类场景也值得自己写一版嵌入式或实时系统需要在流式数据上边采边算每次只处理固定长度的最新窗口算法移植要把 python 原型转成 C 或固件先把 welch 的逻辑用纯数组操作拆开才好逐行翻译教学和论文复现需要明确每一步的数学含义不能直接调库当黑盒自定义改造比如想用非均匀频率轴、自定义重叠策略、或者对每一段做不同的预whiteningscipy 不给你这个自由度。说白了复现不是为了“重复造轮子”而是为了把一个 API 变成一套自己能完全控制的算法。理解到这一层后面所有参数细节才有意义。2. welch 的算法骨架分段、加窗、平均到底在干什么2.1 从周期图到 welch解决两个问题先回顾一个最朴素的做法把一整段信号直接做 FFT然后取模平方再除以频率分辨率得到的就是周期图periodogram。这种估计有两个很明显的毛病。第一是方差大。信号越长FFT 点数越多谱线上的起伏越剧烈。你可以把周期图理解成“只抽了一个样本就来估计整体分布”单次实现和真实功率谱之间可能差出好几个数量级而且这个误差不会随着 FFT 点数增加而自动收敛。第二是谱泄漏。因为 FFT 默认把信号当成周期延拓如果信号长度不是信号周期的整数倍首尾衔接处会产生虚假的间断能量就会从真实频率“漏”到旁边一堆频点上。welch 的核心思路很简单既然一整段信号做的周期图不可靠那就把长信号切成好多段每段单独做周期图最后求平均。分段让估计方差大致降为原来的 1/KK 是段数平均则把随机起伏抹平。代价是每段长度 nperseg 变短了频率分辨率变差所以 nperseg 本质上是在“分辨率”和“方差”之间做权衡。2.2 三个关键参数nperseg、noverlap、windownperseg每段的长度决定 FFT 点数和频率分辨率。分辨率约等于 fs / nperseg段数越多方差越小但每段太短会让两个靠近的频率峰分不开。noverlap相邻两段之间重叠的样本数。如果完全不重叠段数就少重叠起来段数变多方差进一步降低。但重叠太多会让相邻段高度相关新信息减少。window每段加窗的形状用来抑制谱泄漏。welch 要求必须加窗因为如果不加窗分段就等效于对原始信号用了一个矩形窗截断旁瓣泄漏会很难看。这三个参数不是独立的。工程默认值通常取 noverlap nperseg // 2也就是 50% 重叠配合 hann 窗。这个组合经验上性价比最高不是随便拍的。下面会专门讲原因。2.3 归一化为什么是除 fs 和窗能量复现 welch 最容易出错的地方就是归一化所以这里单独拎出来说清楚。SS 也就是功率谱密度物理含义是“单位频率宽度上的功率分布”。FFT 得到的是离散频谱每个 bin 的宽度是 fs / nperseg。如果你把一个周期图的模平方直接当成功率那它其实是“落在这个 bin 里的功率”而不是“单位频率上的密度”。要想变成密度必须再除以 bin 宽度也就是除以 fs / nperseg。但由于窗函数的存在信号能量会被窗改变所以还要把窗的能量贡献除掉。density 模式的归一化因子是scale 1.0 / (fs * sum(w^2))其中 sum(w^2) 就是窗的等效能量。除完以后单位是 V^2/Hz宽带信号比如噪声在这个模式下最能体现真实功率密度。spectrum 模式则不一样它的归一化因子是scale 1.0 / (sum(w^2))这里不除以 fs所以输出单位是 V^2代表“这段信号在这个频率附近的功率”。这个模式更适合测量离散谱线的幅度因为读数不依赖频率分辨率。很多教程只说“welch 要除以窗能量”不提为什么要除以 fs导致复现时把 density 和 spectrum 混为一谈。实际验证时只要对白噪声信号做 PSD 积分你会发现只有 density 模式下积分值约等于信号方差而 spectrum 模式下直接看单频正弦的峰值才更有意义。3. 手写复现完整代码与参数换算3.1 先搭分段索引别用循环复现的第一步不是 FFT而是把长信号切成重叠的段。这里的技巧是用“广播索引”一次性生成所有段的位置而不是写 for 循环一段一段切片。N x.shape[-1] step nperseg - noverlap K 1 (N - nperseg) // step idx np.arange(nperseg)[None, :] step * np.arange(K)[:, None] segments x[idx] # shape (K, nperseg)K 就是这么算的第一个段从 0 开始之后每步移动 step 个样本。如果剩下的样本凑不满一个 nperseg那就直接丢弃。这个行为跟 scipy 保持一致。注意 step 必须大于 0也就是 noverlap 必须小于 nperseg。如果你写 noverlap nperseg那 step 0K 会变成无穷大必须提前拦下来。scipy 会报 ValueError自己实现时也要加这个检查。3.2 去趋势、加窗、FFT一个都不能少分段完成后接下来是去趋势。默认是 constant也就是每一段减掉自己的均值。这一步对低频漂移很有效但如果你用的是高频信号影响很小。linear 模式则是每一段减去一条最小二乘拟合直线对缓慢变化的基线更有效。我自己实现的时候linear 去趋势没有用 polyfit而是用标准正交化的小技巧。对每段长度 M构造一个线性向量 t把它归一化然后从原始段里减去“投影到 t 方向”的分量。这样数值稳定性好而且可以直接向量化跑在多段上。加窗这一步更关键。窗函数必须和分段后的长度一致。scipy 内部对字符串窗名会用 get_window 生成默认是 periodic 模式不是 numpy 里常见的 symmetric 模式。这两个模式的区别在 5.1 节会细说你现在只要知道复现 scipy 的 welch必须用 periodic 窗否则对拍就对不上。然后每一段乘以窗再做 rfft实信号的单边 FFT。对二维 segments 来说axis-1 是段内的时间轴axis0 是段序号轴。3.3 单边谱的 2 倍补偿和频率轴实信号的 FFT 结果是对称的负频率部分是正频率部分的镜像。welch 默认 return_onesidedTrue只保留 0 到 Nyquist 频率但镜像部分的能量不能丢掉所以要把除了直流和 Nyquist 之外的所有频点乘以 2。这里有个非常容易踩的边界当 nperseg 是偶数时Nyquist 频点恰好存在但它没有镜像频率所以不能被乘 2。用代码表示就是pxx np.abs(fft_seg) ** 2 * scale if nperseg % 2 0: pxx[..., 1:-1] * 2 else: pxx[..., 1:] * 2奇数长度时没有 Nyquist 频点所以从索引 1 开始全部乘 2。频率轴直接用 np.fft.rfftfreq(nperseg, d1/fs) 生成干净利落。3.4 完整的 my_welch 实现最后把上面所有步骤收拢成一个函数。为了不把复杂度拉太高这个版本先处理一维信号。import numpy as np def my_welch(x, fs1.0, windowhann, npersegNone, noverlapNone, scalingdensity, return_onesidedTrue, detrendconstant, averagemean): x np.asarray(x, dtypefloat) if x.ndim ! 1: raise ValueError(当前实现只处理一维信号) N x.shape[-1] if nperseg is None: nperseg min(256, N) if nperseg N: raise ValueError(nperseg 不能大于信号长度) if noverlap is None: noverlap nperseg // 2 step nperseg - noverlap if step 1: raise ValueError(noverlap 必须小于 nperseg) # 窗函数字符串时用 periodc hann保证和 scipy 一致 if isinstance(window, str): if window hann: n np.arange(nperseg) win 0.5 - 0.5 * np.cos(2.0 * np.pi * n / nperseg) else: # 其他窗可以调用 scipy.signal.get_window 或手写对应公式 raise NotImplementedError(示例只实现 hann) else: win np.asarray(window, dtypefloat) if len(win) ! nperseg: raise ValueError(窗长度必须等于 nperseg) # 分段 K 1 (N - nperseg) // step idx np.arange(nperseg)[None, :] step * np.arange(K)[:, None] seg x[idx] # 去趋势 if detrend constant: seg seg - seg.mean(axis-1, keepdimsTrue) elif detrend linear: t np.linspace(-1.0, 1.0, nperseg, endpointTrue) t / np.linalg.norm(t) coef seg t seg seg - coef[..., None] * t elif detrend is not False and detrend is not None: seg detrend(seg, axis-1) # 归一化因子 win_energy np.sum(win ** 2) if scaling density: scale 1.0 / (fs * win_energy) elif scaling spectrum: scale 1.0 / win_energy else: raise ValueError(scaling 只支持 density 或 spectrum) # FFT 与单边谱补偿 if return_onesided: fft_seg np.fft.rfft(seg * win, axis-1) freqs np.fft.rfftfreq(nperseg, d1.0 / fs) pxx np.abs(fft_seg) ** 2 * scale if nperseg % 2 0: pxx[..., 1:-1] * 2.0 else: pxx[..., 1:] * 2.0 else: fft_seg np.fft.fft(seg * win, axis-1) freqs np.fft.fftfreq(nperseg, d1.0 / fs) pxx np.abs(fft_seg) ** 2 * scale # 平均 if average mean: pxx pxx.mean(axis0) elif average median: pxx np.median(pxx, axis0) else: raise ValueError(average 只支持 mean 或 median) return freqs, pxx这段代码把核心链路都串起来了。你要是只想拿一版能跑的实现复制上面这段就够了。但如果你要放进生产环境或者继续改造我建议再往下看几节因为边界条件才是真正省时间的地方。4. 与 scipy.signal.welch 对拍验证4.1 用正弦加噪声构造测试信号写代码是一回事能不能跟 scipy 对上是另一回事。我习惯的做法是构造一个“已知答案”的测试信号一个幅度 1.0 的 55 Hz 正弦叠加一个小白噪声。噪声既能考验 PSD 的密度归一化也顺便看看宽带底噪是否一致。import scipy.signal import numpy as np fs 1000 rng np.random.default_rng(42) t np.arange(0, 2.0, 1.0 / fs) x np.sin(2.0 * np.pi * 55.0 * t) 0.01 * rng.standard_normal(t.size) f1, p1 scipy.signal.welch(x, fsfs, nperseg512, noverlap256) f2, p2 my_welch(x, fsfs, nperseg512, noverlap256) print(np.max(np.abs(p1 - p2))) print(np.allclose(p1, p2))我在本机跑的结果是 max abs diff 在 1e-15 量级allclose 直接通过。这说明最基本的 density 模式已经对齐。4.2 逐点对比与误差分析为什么要看最大绝对误差而不是只调 allclose因为 allclose 有默认的 atol有时候差异很大也可能被相对容差放过。直接看 max abs diff 更直观能暴露出“大部分点对得上某几个点差很多”的问题。我之前踩过一个典型坑把窗函数用成了 numpy.hanning也就是 symmetric 模式。结果对拍时发现大多数频点误差很小唯独谱峰附近出现明显偏差最大误差到了 1e-2 量级。symmetric 窗和 periodic 窗在两端相差一点点但就是这一点点会让谱峰的能量分配出现细微差异。所以验证的时候别只看整体误差也要单独看峰值区域的误差。我最后把窗函数写死成 periodic 形式n np.arange(nperseg) win 0.5 - 0.5 * np.cos(2.0 * np.pi * n / nperseg)这才是 scipy 在用字符串 hann 时实际采用的窗。4.3 边界条件奇数 nperseg、detrend、spectrum 模式基础用例对齐之后还要过一组边界用例。下面这几个是我每次改完代码都会跑一遍的“回归清单”场景参数组合预期结果奇数段长nperseg255, noverlap127单边谱补偿逻辑走奇数分支无重叠nperseg256, noverlap0段数为 (N-nperseg)//256 1线性去趋势detrendlinear低频基线被移除低频分量衰减spectrum 模式scalingspectrum正弦峰值更接近线谱功率非 hann 窗windowblackman窗能量按 blackman 计算这些用例不一定每次都打印出 1e-15但只要逻辑对误差都会在 1e-12 以下。如果你在某一项上突然冒出 1e-3 级别的差异先别急着改 FFT多半是归一化因子或者乘 2 补偿写错了位置。还有一个很多人忽略的点scipy 的 detrend 参数支持传可调用对象比如你自己定义的去趋势函数。如果你传了函数它会被应用到每一段上。复现的时候不要丢掉这个能力因为实际信号预处理经常需要自定义去趋势。5. 复现之后踩过的坑和调参心得5.1 窗函数长度与 periodic / symmetric 的坑这个坑值得单独写一小节。scipy.signal.get_window 的默认行为是 fftbinsTrue也就是生成 periodic周期延拓友好窗而 numpy.hanning、numpy.hamming 这类函数返回的是 symmetric 窗。两者的公式差一个分母periodic: w[n] 0.5 - 0.5 * cos(2πn / N) symmetric: w[n] 0.5 - 0.5 * cos(2πn / (N-1))N 相同时两种窗在绝大多数点上的值几乎一样只有两端有差异。但这个差异在频谱上会产生可以观察到的泄漏差别尤其是对高动态范围的信号比如强单音旁边有个微弱信号spectrum 细节就会暴露问题。所以我的建议是手动实现 welch 时干脆统一用 periodic 公式别用 numpy 自带的窗函数。这样和 scipy 对拍最稳也最贴合 FFT 周期延拓的假设。5.2 重叠率不是越高越好很多初学者觉得重叠率越高段数越多方差越小于是直接把 noverlap 拉到 nperseg - 1。实际上这样做会让相邻段之间高度相关平均带来的方差改善会迅速饱和而且计算量成倍增加。以 hann 窗为例50% 重叠已经能利用窗函数两端的衰减信息是工程上公认的“甜点”。如果信号长度允许我更倾向于保持 50% 重叠优先增加 nperseg 或信号总时长来降低方差。只有当信号很短、没法延长时才会考虑把重叠率提到 75%。再往上就没什么性价比了。你可以自己做个实验固定信号和 nperseg分别用 0%、50%、75%、90% 重叠跑几遍拿相同随机种子对比 PSD 曲线的抖动幅度。实测下来0% 到 50% 的改善最明显50% 到 75% 还有一点75% 以上基本是心理安慰。5.3 从 PSD 还原幅度谱时别少乘 2复现 welch 完成后很容易顺手拿来测正弦幅度然后发现读数对不上。这里要给一个提醒density 模式下的 PSD 峰值并不直接等于信号幅度的平方它是经过了窗函数展宽和频率分辨率影响的密度值。如果你的目标是测单频信号幅度直接用 scalingspectrum 模式更合理因为它不做除以 fs 的密度归一化。想看幅度谱还要注意你拿的是单边谱还是双边谱至少把 2 倍补偿那一步想清楚。这块我早期吃过亏总把 density 模式的 PSD 峰值当成幅度平方结果总是差一个和 nperseg 有关的系数后来才明白是窗函数和谱密度共同作用的缘故。我个人现在的做法是测线谱用 spectrum测宽带噪声底噪用 density两者不混着看。这也是生产环境里最常见的分工。复现 scipy 的 welch 这件事技术上不难难的是把所有看似“隐式”的默认行为都复刻到位。写完这版纯 numpy 实现之后我自己对分段索引、窗能量归一化和单边谱补偿的理解明显比原来只调库的时候深了一个层次。如果你也在做类似的事情建议在动手前先列一个参数对照表把 nperseg、noverlap、window、detrend、scaling、return_onesided 这些选项全部钉死再去写代码能少走不少弯路。