做现场调试那会儿我遇到过一次很典型的折腾。一套旋转机械的振动监测系统传感器贴好了数据采回来了时域波形看着跟心电图似的挺正常。可设备就是时不时发出一种沉闷的“呜噜”声领导急着要结论。我随手把振动信号丢进FFT出来的频谱图毛刺多得跟刺猬一样密密麻麻全是尖峰根本分不清哪些是真实特征哪些是统计波动。后来换成Welch方法做功率谱密度估计谱线一下子就平顺了轴承外圈故障特征频率清清楚楚摆在眼前。从那以后Welch方法就成了我处理工程信号时的默认选项。这篇文章就把Welch方法从头到尾拆一遍它到底解决了周期图法的什么问题参数怎么选代码怎么写实际调试中会遇到哪些坑。内容适合正在做信号处理相关工作的工程师、研究生以及任何需要从一段复杂信号里提取频域特征的人。就算你之前没接触过功率谱密度估计只要跟完这篇文章也能上手跑出可用的谱。1. 为什么周期图法不够用Welch方法解决的核心问题1.1 直接用FFT做功率谱谱线为什么像刺猬在那个现场之前我一直习惯对整段信号直接做FFT这其实就是所谓的周期图法。它的原理很简单对一个N点信号做傅里叶变换然后取模的平方再除以N就得到了功率谱Pxx(f) (1/N) * |X(f)|²这个公式好懂也好写但它有一个很麻烦的统计缺陷。一个N点信号能做N/2根谱线数据越长谱线越密但每根谱线周围的统计波动并不随N增大而减小。换句话说你把采样时间从10秒延长到100秒频率分辨率是提高了但谱线周围的毛刺还是那么多只是被压缩到了更窄的频率区间里。专业一点的说法是周期图法不是一致估计它的方差不随数据长度增加而趋近于零。我这样理解这个问题的每个频点好比一个射击选手他单次射击的成绩波动很大你让他多打几枪每一枪的成绩还是波动那么大问题是他多打几枪并不能降低单枪的随机误差。除非你把这几次成绩取平均才能得到一个更稳定的评估。这个直觉正好就是Welch方法改进的核心逻辑。1.2 Welch的三板斧分段、加窗、平均Welch方法做的事情等于是把上面那个射击例子落到了实地上。它不直接对整段信号做FFT而是先把长信号切成若干小段允许相邻段之间有重叠然后对每一段加一个非矩形窗函数再做FFT得到各自的周期图最后把所有段的功率谱取平均作为最终的功率谱密度估计。拆开来看这三步各有各的目的分段加平均把一段长数据变成K段独立或有相关性的子段分别做周期图再平均估计方差压到原来的1/K左右。这是解决毛刺问题的根本手段。重叠直接分段会有一个副作用——窗函数把每段两端的数据压扁了造成数据浪费。允许相邻段重叠可以让上一段衰减掉的部分正好由下一段的头部补回来提高了数据利用率。加窗直接切段等于默认用了矩形窗旁瓣泄漏问题非常严重。换成汉宁窗这类旁瓣低的窗函数谱线之间的泄漏会大幅减少代价是主瓣变宽一点。这三个动作合在一起就是Welch方法比周期图法实用的根本原因。它牺牲了一部分频率分辨率换来了谱估计的稳定性和抗泄漏能力。工程上稳定性有时候比分辨率更重要——你谱峰值再准如果周围毛刺导致根本找不准峰的位置那也是白搭。1.3 方差到底降了多少一个粗略的量化我估计不少人对“方差降低”这个概念比较模糊总觉得这是个定性说法。这里给一个粗略的量化假设你把N点信号分成K个互不重叠的段每段长度是N/K各段周期图的平均会让估计方差降到单段周期图方差的1/K。这就是Bartlett方法的核心思想Welch方法相当于它的加强版。Welch的加强体现在重叠加上加窗之后可以切出更多的段。比如你的一段信号总长10000点nperseg设为1000点如果不重叠能切10段如果改成50%重叠每500点滑动一次能切出19段。虽然重叠段之间并非完全独立方差改善不会严格到1/19但实际效果非常可观。对有经验的人来说50%重叠配合汉宁窗等效的独立段数大约是不重叠段数的1.8到1.9倍几乎是翻倍的效果。这也是为什么工程上默认推荐50%重叠成本低收益大。2. 参数怎么定窗口、分段长度、重叠率的取舍逻辑2.1 窗函数选型汉宁、汉明、布莱克曼怎么挑Welch方法里窗函数的选择核心还是那句老话主瓣宽度和旁瓣衰减之间的权衡。主瓣越宽频率分辨率越差旁瓣越高能量泄漏越严重。下面是几种常用窗的对比窗函数主瓣宽度相对矩形窗旁瓣水平典型用途矩形窗基准-13 dB瞬态信号、频率分辨率优先汉宁窗Hann约2倍-31 dB通用默认兼顾分辨率和泄漏汉明窗Hamming约2倍第一旁瓣约-43 dB音频窄带信号但远端旁瓣衰减慢布莱克曼窗Blackman约3倍-58 dB动态范围要求高、需要极低泄漏凯塞窗Kaiser可调可调对不同需求灵活调节我自己的经验是不知道选什么就选汉宁窗。它主瓣比矩形窗宽一倍但旁瓣从-13dB压到-31dB这个交换在大多数工程场景下非常划算。布莱克曼窗虽然旁瓣更低但主瓣太宽两个频率接近的谱峰会糊成一片得不偿失。汉明窗偶尔在对幅值精度要求高、且已知信号频带很窄的时候用但它的远端旁瓣衰减不如汉宁干净所以通用性弱一些。这里要说一个容易忽略的点加窗之后单频正弦在谱上不再是理论上的一根孤立线而是变成了窗函数频谱的形状——中间一个峰两边带旁瓣。这是必然的能量展宽不是算法出错了。你看到的是一个峰它的宽度就是窗函数主瓣它的底部的花纹就是窗函数旁瓣。理解了这一点再去调窗函数的参数就顺理成章了。2.2 重叠率50%起步为什么推荐重叠重叠率的设计本质上是为了应对窗函数带来的边缘衰减。以汉宁窗为例它两端接近0如果各段首尾相连不重叠每段两端的约25%数据几乎不贡献有效信息整段数据的实际利用率只有一半左右。这就太浪费了。Welch把相邻两段的起始位置拉开一定距离比如noverlap设为nperseg的一半也就是50%重叠。这时候上一段被压扁的区域正好是下一段的高权重区域数据利用率能提升到85%以上。这是一个接近“白捡”的收益计算量只增加了不到一倍方差改善却接近翻倍。所以我的建议是如果没有特殊原因重叠率不要低于50%。那是不是越高越好不完全是。75%重叠还能再榨出一些方差改善等效独立段数大约到不重叠段的2.6到2.8倍比50%再提升约50%但计算量是50%重叠的1.5倍。到87.5%以上新增的计算量换来的方差改善就非常有限了因为相邻段之间的相关性太高信息冗余严重。工程上我的默认路线是先试50%谱线已经够干净就维持谱线还是粗糙就提高到75%实时系统扛不住计算量就降到30%到50%之间取一个平衡点。2.3 频率分辨率nperseg说了算nfft只是插值很多初学者特别容易被一个参数绕晕nfft。看到scipy.signal.welch里有个nfft参数以为把nfft调大就能提高频率分辨率结果发现谱线变密了但两个相邻的谱峰还是分不开。这不奇怪因为真实的分辨率根本不取决于nfft它取决于每一段参与FFT的有效长度nperseg。频率分辨率的公式是Δf fs / nperseg这个公式的含义是如果采样率固定一段信号能分辨的最小频率间隔取决于这段信号的实际时长。nperseg1000点、fs1000Hz那Δf就是1Hz你想分辨两个相隔0.5Hz的峰至少要2秒的段。如果信号有效采样率下做不到那再大的nfft也只是在已有的谱线上做插值谱线之间画出来的“峰”全是假象。我做过的真实案例可以说明问题现场采了一段轴承振动信号fs512Hz采样时长8秒怀疑有两个故障特征频率相差约1.8Hz。我一开始把nperseg设成256Δf2Hz两个峰完全糊成一个后来把nperseg增大到1024Δf0.5Hz两个峰才清晰分开。但代价是段数从31段降到了8段多一点谱线略糙。这个权衡你就得根据实际情况拍板分辨率优先还是稳定性优先。2.4 段数与总时长数据不够平均无从谈起Welch的“平均”是需要素材的。一段信号能切出多少段有个简单公式K floor((N - nperseg) / step) 1其中step nperseg * (1 - overlap_ratio)overlap_ratio是重叠比例。举个例子总长度N10000点nperseg1000重叠50%step500K (10000-1000)/500 1 19段。19段平均出来的谱相当平稳。但如果你手里的数据总共只有2000点nperseg还非要取1024那最多只能切2到3段。你是在用大量重叠去制造“伪段数”实际上等效的独立段数很低方差压不下来。这时候我的建议是要么接受粗糙一点的分辨率把nperseg降下来多换几段要么去更长时间的数据要么就考虑其他谱估计手段比如自回归AR类方法或者Multitaper方法它们在短数据下的表现通常比Welch好。记住Welch方法的前提是数据量足够支持你“挥霍”。3. 实操流程从原始信号到功率谱密度的完整实现3.1 数据预处理去趋势、去直流、剔除异常尖刺第一步是去趋势项。传感器温漂、电路基线漂移都是低频干扰的大户如果不处理谱图低频段会出现一个巨大而宽的假峰严重时能把整个谱图的动态范围压扁。scipy.signal.welch本身提供了一个detrend参数默认是去均值但这只够对付直流分量。如果信号基线明显上下漂我习惯先单独做一个线性趋势去除甚至多项式拟合去除再送进Welch。第二步是剔除异常尖刺。采集卡受到电磁干扰、数据传输丢包可能产生个别幅度极大的点。这些点在时域上看着只是一个小毛刺但在频域上会抬高整个频带的底噪掩盖小信号的真实幅度。一个简单粗暴的办法是先查看时域波形的最大值如果某几个点比整体RMS高一个数量级以上用中值滤波或者直接标红删除再进Welch。第三步是确认采样率。这里的坑主要是在做信号回放或者数据格式转换时fs写错了导致整个频率轴都偏掉。先喊停一秒钟在跑Welch之前花30秒确认fs这个数字是从原始硬件配置里来的而不是拍脑袋写的。3.2 scipy.signal.welch参数逐项拆解下面这段代码是我调试现场振动信号时随手整理的模板简化后放在这里每个参数都配上解释import numpy as np from scipy import signal import matplotlib.pyplot as plt # 采样率与时间轴 fs 1000 t np.arange(0, 10, 1/fs) # 构造一个含50Hz和120Hz正弦、叠加噪声的测试信号 x 1.0 * np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 120 * t) x 0.2 * np.random.randn(len(t)) # Welch功率谱密度估计 f, Pxx signal.welch( x, fsfs, # 采样率决定频率轴的范围 windowhann, # 窗函数默认建议汉宁 nperseg1000, # 每段长度决定频率分辨率Δffs/nperseg noverlap500, # 重叠点数取nperseg的一半即50%重叠 nfftNone, # 默认等于nperseg一般不额外设置 detrendconstant, # 去均值若基线漂移改为linear return_onesidedTrue, # 实信号只返回0~fs/2的单边谱 scalingdensity # 返回功率谱密度单位V^2/Hz ) # 用对数坐标绘制PSD plt.semilogy(f, Pxx) plt.xlabel(Frequency (Hz)) plt.ylabel(PSD (V^2/Hz)) plt.grid(True) plt.show()参数背后我的取值逻辑是这样信号总长10秒fs1000HzN10000点取nperseg1000点意味着Δf1Hz能清晰区分相隔70Hz的50Hz和120Hz两个峰同时还能切出19段做平均。这个组合既有分辨率又有稳定性在试跑阶段非常稳妥。等你对信号的频段特征更了解了再去针对性地调整nperseg和noverlap。noverlap取500就是50%重叠前面讲过这是汉宁窗下的经典推荐值。detrend我默认选constant如果看到低频段有鼓包再回去改linear。nfft保持None它跟nperseg一样除非你是为了画图的平滑性想插值否则没有理由单独设置它。3.3 从PSD提取工程指标频段能量与RMS跑出PSD只是第一步工程上往往还需要从谱里提取具体的指标。最常用的一个指标是某个频段的总能量。因为PSD是“单位频率上的功率”对一个频段积分就得到该频带内的总功率。离散化处理时积分变成求和再乘频率间隔# 计算1~200Hz频段内的总RMS band (f 1) (f 200) df f[1] - f[0] rms_1_200 np.sqrt(np.sum(Pxx[band]) * df) print(f1-200Hz频段总RMS {rms_1_200:.4f})这里的重点在于理解单位。PSD的单位是信号幅值单位的平方再除以Hz比如电压信号就是V²/Hz。求和乘df之后单位回到V²开方就是V。如果信号是加速度信号那单位就是(m/s²)²/Hz开方后得到的RMS就是某个频段的振动加速度有效值。很多刚接触PSD的人搞不清这个链条结果在报告里单位写错非常尴尬。另外一个高频操作是提取某个频带内的能量占比。比如在轴承故障诊断中我常常把故障特征频率附近若干阶谐波的能量加起来跟整个高频段的能量做比值作为一个无量纲的退化指标。这个比值的变化趋势比单根谱线的绝对幅值稳定得多适合做趋势监测。3.4 画图显示为什么用对数坐标这是一个没什么技术含量、但非常影响工作体验的细节。功率谱密度经常跨越好几个数量级比如一个50Hz的大峰可能有1e-2 V²/Hz量级而底噪只有1e-6 V²/Hz量级。用线性坐标画底噪那条线会直接贴到横轴上啥也看不见只有一个个大尖峰。用semilogy或者loglog把纵轴换成对数刻度底噪起伏一目了然小峰也能看清楚。频率轴也建议按需截取。你如果关心的是0到500Hz的宽频段那就全画如果关心某个窄带故障频段比如2000到3000Hz就把横轴切到这一段。有经验的工程师从来不会把一张完整的0到fs/2谱图直接甩给领导看而是先把关心的频段放大、标注好特征峰值再交出去。这跟用不用Welch没有直接关系但你做出来的图清晰别人对你的结果信任度自然会高不少。4. 常见问题与排查技巧实录4.1 谱线毛糙得没法看先查这两项毛刺多最直接的原因就是参与平均的段数太少。这时候第一反应不是换算法而是先检查nperseg是不是设得太高了。比如你总共有20000点数据nperseg直接取8192只能切3到4段谱线粗得跟砂纸一样。把nperseg调到2048段数能到10段左右效果立刻改善。第二反应是检查重叠率noverlap是不是设成了0或者特别低。如果是从别处拷贝来的代码特别容易碰到noverlap忘记改成合适的值一直用默认设置跑。我踩过一个具体的坑从实验室旧代码里继承了nperseg1024、noverlap0的配置跑出来的谱峰位置是对的但背景特别糙一度以为传感器坏了或者环境干扰太大。后来把重叠率提到512点情况立刻改善底噪平顺了两三倍。所以遇到毛糙先动段数和重叠率不要急着怀疑硬件。4.2 分辨率不够怎么办物理约束与方法选型频率分辨率不足的时候谱上的两个峰可能只是一个鼓包旁边还拖着一条往下的尾巴。优化方向分两步走先看物理上能不能延长每段长度。只要总时长允许把nperseg加大Δf就降下来这是最正统的思路。如果总时长就那么多加大nperseg会直接导致段数减少谱更糙——这时候你面临的是Welch方法固有的分辨率和方差的跷跷板。如果两侧都调不动那就只能考虑换方法了。比如Burg法自回归模型谱估计在短数据下分辨率明显优于Welch但阶次选择需要谨慎阶数太高容易出现虚假谱峰Multitaper方法在频谱同时包含尖峰和平滑底噪时表现均衡但实现复杂参数也多。我在实际项目里遇到过非要用Welch做1秒短数据、还想分辨1Hz之内的两个峰的需求最后物理上做不到劝用户换了采集方案延长了采样时间才彻底解决。方法选型再怎么折腾也抵不过数据本身更合适。4.3 低频大鼓包趋势项和泄漏的误判低频段出现超级大的宽峰是现场数据最常见的状况。你从时域图看信号基线可能根本没有固定在0附近而是缓慢漂移的。这些低频成分的能量主要集中在0到几赫兹在谱图上表现为一个巨大的鼓包把其他频段的谱线都压扁了。碰到这种情况先别急着怀疑轴承、齿轮有什么低频故障。先做两件事一是把detrend设为linear去除线性趋势二是给信号加一个高通滤波器比如在0.5Hz或1Hz以下切掉。做完这两步再跑Welch低频鼓包通常会大幅缩小。如果鼓包还在那才值得怀疑存在真实的物理低频激励。这个排查顺序我已经用了很多年基本百试百灵。4.4 等间距尖峰先确认是外部干扰还是信号特征谱图上出现一系列等间距的尖峰第一反应是看基频是不是50Hz或60Hz。如果是那大概率是工频干扰及其谐波不是设备本身的振动特征。判断方法很直接让设备停机或者让传感器处于静止状态重新采一段背景噪声再跑一次Welch。如果等间距尖峰在背景谱里依然存在那就是环境干扰或采集系统自身的电磁耦合需要在硬件层面做屏蔽和接地处理。反过来如果等间距尖峰的基频跟设备的某个转频或齿轮啮合频率对得上那才是信号的真实特征。比如一个轴的转频是25Hz谱上在25、50、75、100Hz处都有尖峰而且依次衰减这大概率是转频的各次谐波常见于不对中、松动等故障。这种问题的关键是千万别只盯着第一根峰谐波结构才是更可靠的判据。下面整理一份排查速查表方便现场直接翻症状可能原因排查方向整条谱线毛糙、起伏大段数太少或重叠率过低减小nperseg、增大noverlap到0.5或0.75两个相邻峰糊成一个频率分辨率不够增大nperseg必要时延长采样时长低频段出现巨大鼓包趋势项、基线漂移开detrendlinear或加高通滤波等间距尖峰工频干扰或其谐波采背景信号对照确认是否外部干扰谱峰幅值明显偏低窗函数能量归一化差异确认scaling参数以及窗能量sum(window^2)某个峰旁边出现对称小峰窗函数旁瓣泄漏换旁瓣更低的窗如Blackman5. 应用场景与其它谱估计方法对比5.1 机械振动故障诊断中的实际用法滚动轴承故障是Welch方法用得最频繁的领域之一。轴承外圈、内圈、滚动体故障时会在特征频率上产生能量集中的振动调制。我这里说的是包络谱的做法先对振动信号做带通滤波通常在轴承结构共振频段选带然后对这个带通信号取包络再对包络信号做Welch谱估计。因为故障信号往往被高频载波调制了直接对原始信号做Welch可能什么都看不出来但一旦取包络再谱分析故障特征频率就清晰出现了。具体流程是确定轴承几何尺寸和转频计算BPFO、BPFI这些特征频率对信号做带通滤波比如500到5000Hz用Hilbert变换取包络对包络用Welch估计功率谱在谱里找特征频率的峰值和谐波结构。每一步单独看都不复杂但串起来就需要对信号处理有整体把握。Welch在中间承担的角色是“提供一个稳定、平滑的谱图”让后续的峰值搜索算法不必跟毛刺搏斗。5.2 生物信号与音频分析中的标准操作生物信号领域Welch几乎是脑电和心率变异性分析的默认工具。脑电的alpha波8到13Hz、beta波13到30Hz这些频带能量估计标准做法就是对多导联脑电信号做Welch后取各频段积分。心率变异性研究里的低频功率和高频功率之比本质上也是先对RR间期序列做Welch再算对应频带的功率比例。区别只在于信号本身的频率低很多通常到0.5Hz以内nperseg得按秒级别的时长来设。音频和声学测量里Welch被用来估计噪声频谱、房间声学响应。一个常见应用是环境噪声的1/3倍频程分析它需要把连续频段按倍频程分成多个频带再统计每个频带的能量。Welch的平滑谱线非常适合这种频带积分操作不会因为个别谱线的随机起伏导致相邻频带之间出现莫名其妙的数据跳变。5.3 Welch、周期图、Burg、Multitaper怎么选我把不同谱估计方法放在一起对比过很多次各有各的适用场景方法优势劣势典型场景周期图法简单、快方差大、谱线毛糙快速预览信号频段组成Barttlett法比周期图平稳灵活性差、不能重叠教学示例Welch方法平稳性好、参数灵活、开源库成熟短数据下分辨率受限工程默认首选Burg/AR法短数据下分辨率高模型阶数敏感、易出假峰短数据高分辨率研究Multitaper方差与分辨率兼顾参数多、使用门槛高科研、弱信号提取工程上的决策顺序我一般是这样数据够长默认Welch数据太短且必须高分辨率优先Burg法但一定要做阶数敏感性分析频谱动态范围特别大、强线谱旁边埋着微弱的窄带信号Multitaper更合适如果只是看一眼大概频谱那直接周期图法也无妨省事。选型有个大前提谱估计方法只是工具决定上限的是数据质量和对信号的物理理解。你连信号本身是什么成分都不知道换任何高级方法都是花架子。先花时间搞清信号的物理来源、采样参数、干扰类型再选合适的方法这才是正经工程师的做事顺序。我在实际项目里的最大体会是Welch不是那种看着高深但用不上的数学工具它是信号频域分析的第一道通用工序。参数这东西拿真实数据多试几次比背一百遍文档都管用。如果只记住一个参数组合我推荐汉宁窗、50%重叠、nperseg取总长度的1/10到1/20之间先跑一版看谱形再根据谱线形态去细调。这样做十有八九走不太偏剩下的就靠你对数据和物理场景的理解。另外如果你后面要做实时谱分析Welch的逐段计算结构很适合做流式更新把新来的数据段增量式纳入平均能省下大量重复计算这块之后有机会再展开聊。