简介这份文档围绕凯斯西储大学CWRU轴承故障数据集展开信号时频分析面向从事旋转机械故障诊断、健康监测与振动信号处理的学习者和研究人员帮助其理解如何从实测振动数据中提取可用于故障辨识的时频特征。内容先说明实验台组成与驱动端、风扇端、基座三处加速度计数据的差异再以正常信号和0.021英寸内圈、滚珠、外圈故障信号为对象比较短时傅里叶变换在重叠比例0.5、尺度16/32/64下的时频图表现并给出选择尺度32的依据连续小波变换部分则对比morl、cmor1-1、cmor1.5-2、cgau8等小波函数在尺度128下筛选出cmor1.5-2并进一步讨论32、64、128、256尺度对高低频特征的取舍。文档还附有Python读取mat数据、绘制原始信号与STFT时频图的代码片段便于对照复现和调整参数。资源包为1个docx文件大小约1.05MB已有1584人学习适合作为故障诊断入门与实验参数选择的参考材料。1. 一张时频图能分辨什么轴承故障数据集里被忽略的 85 毫秒很多人拿到某高校公开的轴承故障数据集第一反应是把0_0.mat和21_1.mat读进来跑一行stft撒一把 jet 色图然后在四个子图之间找差别。真正决定这张图能不能分辨内圈、滚珠、外圈的不是用 STFT 还是 CWT而是窗口长度与故障冲击节拍之间的比例。数据集里一个 0.021 英寸约 0.53 毫米的缺陷在 30 Hz 上下的轴转频下会激出周期仅几毫秒的冲击串窗口一旦长过几个冲击周期冲击就被平均成一条平滑亮带故障之间的差别随之消失。反过来窗口短到能咬住每一次冲击频率轴就粗到看不清特征频率。这篇把窗长、尺度、小波函数这三层参数各自在开什么价讲清楚也讲清那份时频分析文档里的参数是怎么被一步步逼出来的。2. STFT 的窗长为什么分辨率这个词在这份数据集上会骗人2.1 短时傅里叶变换在切的到底是什么短时傅里叶变换的做法可以一句话说完拿一个中心位于时间 α 的窗 g(t−α) 在信号上滑动在每个位置只对窗内那一小段做傅里叶变换。窗长固定窗的形状固定于是用多长的尺子去量信号这件事被一次性钉死了。钉死的代价就是不确定性原理时域窗越窄频域主瓣越宽两者乘积有下界。落到代码里这个下界表现为两个可以直接算出来的量——频率分辨率 Δf fs / nperseg和时间分辨率 Δt nperseg / fs。它们不是两个独立可调的旋钮而是同一根杠杆的两端动一个必然动另一个。这份数据集的驱动端加速度计采样率为 12 kHz单次分析截取 1024 点约 85.3 毫秒。85 毫秒里发生了什么值得先算清楚轴转频约 30 Hz一个转频周期 33 毫秒1024 点大约覆盖 2.5 个转频周期内圈故障特征频率的经验系数是 5.4152 倍转频约 162 Hz对应冲击重复周期 6.2 毫秒85 毫秒里大约能装下 13 到 14 次冲击外圈系数 3.5848 倍约 107 Hz周期 9.3 毫秒滚珠系数 4.7135 倍约 141 Hz周期 7.1 毫秒。这三个周期都落在 6 到 10 毫秒这个区间而窗长的候选值从 2.67 毫秒一路跨到 21.3 毫秒正好横跨能分辨冲击和糊掉冲击的分界线。2.2 窗长 16、32、64、256 的账Δf、Δt 与冲击节拍把四个候选窗长的账摆在一起很多结论就不用猜了。npersegΔf fs/npersegΔt nperseg/fs帧移重叠 50%频点个数06 kHz16750 Hz1.33 ms0.67 ms932375 Hz2.67 ms1.33 ms1764187.5 Hz5.33 ms2.67 ms3325646.875 Hz21.33 ms10.67 ms129窗长 16 时只剩 9 个频点整张图几乎是空的谈不上分辨率高低是信息量不够。窗长 256 的频率分辨率最好46.875 Hz理论上能把 162 Hz 和 107 Hz 的谱线分开——但它的时间窗是 21.33 毫秒一个窗里塞进 3 个以上的内圈冲击周期冲击串被卷积成一团时频图上只剩一条持续亮带四种状态长得越来越像。窗长 32 的处境刚好相反375 Hz 的频率分辨率连 162 Hz 与 107 Hz 都分不开可它的时间窗 2.67 毫秒比任何一个故障冲击周期都短每一次冲击都能被单独框住图上呈现的是一串等间隔的竖条间隔 6.2 毫秒还是 9.3 毫秒一眼可辨。注意这里有一个容易想反的地方。窗长 32 之所以在辨识度上胜出靠的不是频率分辨率而是它抓住的是时域冲击节拍不是频域谱线位置。想看清特征频率本身靠 STFT 图本来就不现实得走包络谱那条路。我一般会这么理解这个过程STFT 在这份数据集上实际充当的是一个冲击节奏计它回答每隔多久被敲一下至于敲在哪个频率上那是另一套工具的事。把这个前提想清楚参数的取舍就不再是玄学而是一道有明确约束的算术题。2.3 重叠 0.5 与帧移为什么不是 0.9重叠比例决定帧移帧移决定时间轴上的采样密度。重叠 0.5 意味着帧移等于窗长的一半窗长 32 时帧移 1.33 毫秒1024 点能得到约 64 帧。这个密度够用冲击周期 6.2 毫秒每 1.33 毫秒一帧一个冲击周期里有 4 到 5 帧冲击的起落能被描出来。重叠拉到 0.9帧数涨到 5 倍图更平滑但每一帧用的还是同样长度的窗信息量并没有增加多出来的只是插值密度计算量却翻了五倍。更麻烦的是重叠过高会让相邻帧高度相关时频图上出现横向的拖影看着连续实际是重复计数的假象。重叠降到 0帧移等于窗长冲击很可能正好落在两帧的交界处被两个窗各切一半幅度塌下去一半反而看不清。重叠 0.5 是这份分析的合理落点它让帧移小于最短的冲击周期又不至于制造冗余。import numpy as np from scipy.io import loadmat from scipy.signal import stft FS 12000 # 驱动端加速度计采样率 N 1024 # 单次分析窗长1024 点 ≈ 85.3 ms def load_de(path, key): CWRU 风格的 mat 文件读出来是 dict键名随文件变化必须显式指定 return loadmat(path)[key].reshape(-1)[:N] x load_de(21_1.mat, X209_DE_time) # 0.021 英寸内圈 for nperseg in (16, 32, 64, 256): noverlap nperseg // 2 # 重叠 50% → 帧移 nperseg/2 f, t, Z stft(x, fsFS, npersegnperseg, noverlapnoverlap, windowhann, boundaryNone, paddedFalse) df FS / nperseg # 频率分辨率 dt nperseg / FS # 时间分辨率单帧覆盖时长 print(fnperseg{nperseg:3d} df{df:7.2f} Hz dt{dt*1e3:6.2f} ms bins{len(f)})这段代码只做一件事把窗长扫描的代价量化出来。nperseg是每次傅里叶变换的点数直接决定dfnoverlap只影响时间轴上帧的密度不影响df和dtboundaryNone, paddedFalse关掉 scipy 默认的两端补零避免时间轴被悄悄延长半个窗。跑完打印出来的四行数字就是上面那张表的来源也是判断某个窗长到底行不行的第一手依据不用靠看图猜。3. 连续小波变换cparam 2·fc·totalscal 这行为什么决定了尺度可比性3.1 尺度-平移坐标与时间-频率坐标的换算连续小波变换把 STFT 里固定窗长、滑动位置的两参数结构换成了伸缩尺度、平移位置。小波函数 ψ(t) 被拉伸成 ψ((t−b)/a)a 是尺度b 是平移。尺度 a 变大的时候小波在时域被拉长、在频域被压窄并整体往低频挪于是一个参数同时改变了时间窗和频率窗——STFT 里那根固定的杠杆在这里变成了随频率自适应移动的杠杆。代价是坐标变了。小波变换天然活在尺度-平移平面上画图时要换回频率-时间平面换算关系是 f fc / (a · dt)fc 是小波的中心频率dt 是采样间隔。这个式子看起来平淡却是后面所有坑的源头尺度 a 本身没有物理意义只有乘上 fc 和 dt 之后才有赫兹。尺度与频率成反比这件事直接推出小波变换的双尺度性质。要关注高频特征选小尺度小波被压短时间定位准要关注低频特征选大尺度小波被拉长频率定位准。原文里那句尺度长度越大对低频特征关注越多说的正是这个方向性但如果只停在这句定性描述上调参还是靠手感。3.2 Morlet 的时频盒cmor 为什么能给出平滑包络Morlet 小波是一个复指数乘高斯包络写成 ψ(t) exp(−t²/2)·exp(jω₀t)。高斯包络保证时域局部化复指数保证频域是单峰。这两条性质合起来使它成为振动信号分析里最常用的小波之一——振动是准周期的冲击串既需要时间定位也需要频率定位。复数这个细节值得单独说。实小波比如实 Morlet、墨西哥帽的频谱关于零频对称正负频率各占一半对实信号做变换后取模正负频率的响应会叠在一起互相干涉幅值图上出现细密的振荡纹路。复小波是解析小波频谱只在正频率一侧有支撑负频部分几乎为零取模得到的是干净的包络不会自干涉。这就是为什么同样做 CWTcmor系列的幅值图看起来比morl干净冲击的位置和强度更直观。PyWavelets 里复数 Morlet 的参数写成cmorB-CB 是带宽参数C 是中心频率。cmor1.5-2意味着带宽 1.5、中心频率 2cmor1-1意味着带宽 1、中心频率 1。带宽参数控制高斯包络在频域的宽窄B 越大频域越宽、时域越短频率选择性下降但时间定位变好中心频率 C 决定小波的主频位置。cgau8是 8 阶复数高斯小波它是高斯函数取 8 阶导数没有单独的带宽参数阶数越高频率选择性越好、时域支撑越长。3.3 cparam 2·fc·totalscal 这行代码在做什么原文的尺度序列是这样生成的cparam 2 * fc * totalscal scales cparam / np.arange(totalscal, 0, -1)np.arange(totalscal, 0, -1)产生从 totalscal 递减到 1 的整数序列除法之后 scales 从大到小排列最小值是 cparam/totalscal 2·fc最大值是 cparam 2·fc·totalscal。代入频率换算 f fc/(a·dt)最小尺度 2·fc 对应最高频率 f_max fc / (2·fc·dt) 1/(2·dt) fs/2正好是奈奎斯特频率最大尺度 2·fc·totalscal 对应最低频率 f_min fc / (2·fc·totalscal·dt) fs / (2·totalscal)。两个结果里 fc 都被约掉了。这意味着无论选cmor1-1、cmor1.5-2还是cgau8只要各自的cparam用同一个2*fc*totalscal算频率覆盖区间就完全相同恒为 [fs/(2·totalscal), fs/2]。12 kHz 采样、totalscal128 时区间是 [46.875, 6000] Hz。import numpy as np, pywt FS, totalscal 12000, 128 dt 1.0 / FS for name in (cgau8, morl, cmor1-1, cmor1.5-2): fc pywt.central_frequency(name) # 数值积分得到的中心频率 cparam 2 * fc * totalscal scales cparam / np.arange(totalscal, 0, -1) freqs fc / (scales * dt) # 尺度 → 频率 df np.abs(np.diff(freqs)) print(f{name:11s} fc{fc:6.3f} scales[{scales.min():7.2f},{scales.max():7.2f}] ff[{freqs.min():7.2f},{freqs.max():7.2f}] Hz f相邻尺度频率差: 低端 {df.min():.3f} Hz / 高端 {df.max():.1f} Hz)跑一遍会看到两件事。第一四种小波打印出的f[46.88, 6000.00]完全一致fc列各不相同但都被cparam吸收掉了这正是这行代码存在的目的把中心频率这个自由参数从频率覆盖范围里剔除让换小波只改变时频盒形状不改变分析的频段比较才公平。第二df在低频端只有零点几赫兹在高频端上千赫兹——尺度是等间隔的频率是 1/a 关系所以频点在低频密集、高频稀疏。这是 CWT 的固有特性也是它相比 STFT 更适合看低频故障特征的原因。提示scales必须保持从大到小对应的freqs才会从小到大单调递增。contourf的 y 轴依赖这个单调性顺序反了图会上下颠倒低频跑到顶部看起来和别人的结果不一样。4. 小波函数与尺度实测cgau8、morl、cmor1-1、cmor1.5-2 差在哪4.1 复小波、实小波与辨识度的关系原文在尺度固定为 128 的前提下横向比了四种小波结论是cgau8和cmor1.5-2对不同故障的辨识度较高。这个结论背后的机制可以从三个角度拆。第一个是负频泄漏。morl在 PyWavelets 里是实 Morlet频谱对称取模后有自干涉纹路冲击之间的低能量区被填上振荡条纹三种复小波没有这个问题冲击与冲击之间的背景更干净对比度更高。第二个是时频盒形状。cgau8因为是 8 阶高斯导数时域支撑长、频域尖对窄带成分敏感冲击串出现时能量集中得很紧cmor1.5-2的带宽 1.5 比cmor1-1的 1 更宽时域更短对每次冲击的起落描得更细代价是频率方向的涂抹多一点。在这份数据上冲击节拍是主要判别信息时间定位更值钱所以cmor1.5-2的宽带宽反而成了优点。第三个是参数一致性。cgau8没有可调的带宽参数它的时频盒形状是被阶数锁死的cmor系列可以通过 B 和 C 调形状这意味着cmor的结论能往别的工况迁移cgau8的结论更依赖当前数据的冲击形态。综合下来选cmor1.5-2不是因为它在某张图上最漂亮而是因为它同时满足复小波无负频泄漏和时频盒可调这两条可迁移的性质。4.2 totalscal 从 32 到 256尺度长度变大到底改变了什么原文比较了 totalscal 取 32、64、128、256 时滚珠故障的时频图观察到尺度越长对低频特征关注越多。用上一节的公式可以把这个观察量化f_min fs/(2·totalscal)f_max 恒为 fs/2。totalscalf_minf_max低频端相邻尺度频率差尺度序列长度32187.5 Hz6000 Hz约 0.36 Hz326493.75 Hz6000 Hz约 0.18 Hz6412846.875 Hz6000 Hz约 0.09 Hz12825623.4 Hz6000 Hz约 0.05 Hz256两件事同时发生频率下限往下探同一频段的划分变密。totalscal 从 128 加到 256下限从 46.875 Hz 降到 23.4 Hz180 Hz 以下的低频区从十几个尺度变成三十几个尺度图上低频块的分层明显更细。但高频端一点也不受益6000 Hz 附近的尺度间隔几乎没变还是稀疏的。这解释了图上底部低频的区别最明显。同时也提醒一件事totalscal 不是越大越好。总尺度数翻倍计算量翻倍而新增的分辨率全部堆在低频。如果关心的频段在 1 kHz 以上把 totalscal 加到 256 只是白白多跑一倍时间。4.3 从 loadmat 到 contourf 的完整链路import numpy as np, pywt import matplotlib.pyplot as plt from scipy.io import loadmat FS, N 12000, 1024 dt 1.0 / FS def cwt_map(path, key, wname, totalscal): x loadmat(path)[key].reshape(-1)[:N] fc pywt.central_frequency(wname) cparam 2 * fc * totalscal scales cparam / np.arange(totalscal, 0, -1) # 大到小 → 频率低到高 coef, freqs pywt.cwt(x, scales, wname, sampling_perioddt) amp np.abs(coef) # 复小波取模 → 包络 db 20 * np.log10(amp / amp.max() 1e-12) # 归一化 dB压掉量级差 return freqs, db fig, axes plt.subplots(1, 4, figsize(20, 5), shareyTrue) for ax, ts in zip(axes, (32, 64, 128, 256)): freqs, db cwt_map(21_2.mat, X222_DE_time, cmor1.5-2, ts) t np.arange(N) * dt edge N // 8 # 掐掉两端边界效应 ax.contourf(t[edge:-edge], freqs, db[:, edge:-edge], levels30, cmapjet, vmin-40, vmax0) ax.set_title(ftotalscal{ts}) ax.set_xlabel(t / s) axes[0].set_ylabel(f / Hz) plt.show()几个参数的作用需要交代清楚。sampling_perioddt是pywt.cwt换算频率轴必需的漏掉它默认按 1 秒采样处理频率轴整体错 12000 倍。np.abs(coef)取模是针对复小波的操作实小波取模意义不大。20*log10归一化很关键——原始幅值的动态范围常常跨越三四个数量级线性色标下只有最强的那几个点有颜色其余全被压成背景色看起来没有特征实际是显示问题。vmin-40, vmax0把显示窗口固定住四张图才能横向比较否则每张图各自归一化亮度的相对关系就丢了。edge N // 8裁掉两端因为边界处小波支撑越界会出现两条亮的竖边不裁掉会误以为是故障冲击。5. 避坑与排查五类把时频图做废的现场5.1 坐标轴类时间轴和频率轴最容易一起错第一种情况时间轴长度对不上或者刻度的物理含义完全错。现象是横轴标着秒但从 0 到 1024看不出任何冲击间隔。原因是把窗长当成了采样率——原文代码里sampling_rate 1024紧接着t np.linspace(0, 1.0/sampling_rate, sampling_rate)这个 1024 是单次分析的点数不是每秒的采样点数。12 kHz 采样下 1024 点只覆盖 85.3 毫秒时间轴应该是 0 到 0.0853 秒。解决方法是单独定义一个FS 12000时间轴用np.arange(N) / FS别让窗长和采样率共用一个变量名。第二种情况CWT 图上下颠倒低频跑到顶部。现象是和别人同样的参数图却像镜像。原因是scales没有按从大到小排列pywt.cwt返回的freqs随之从大到小contourf按数组顺序铺 y 轴图就翻了。解决方法是坚持用cparam / np.arange(totalscal, 0, -1)生成尺度画图前打印一下freqs[:3]和freqs[-3:]确认单调递增。5.2 参数可比性类换个小波尺度就不是原来的尺度第三种情况跑morl那一组时直接抛NameError: name fc2 is not defined。现象是前三组正常第四组或者第二组一运行就崩。原因是fc2 pywt.central_frequency(wavename2)那一行被注释掉了但下一行cparam2 2 * fc2 * totalscal还在用这个变量。照抄代码时最容易踩这个坑而且报错位置和真正的问题不在同一行。解决方法是把注释行放开或者在循环里统一用fc pywt.central_frequency(name)现算从结构上消灭fc1、fc2、fc3、fc4这种容易错位的命名。第四种情况直接把scales np.arange(1, 129)喂给pywt.cwt然后拿不同小波的结果横向比。现象是每种小波的图频率范围都不一样cmor1.5-2的图明显比cmor1-1更高没法比。原因是没有用cparam 2*fc*totalscal做归一化f fc/(a·dt)里的fc没有被吸收中心频率 2 的小波频率整体高一倍。解决方法是统一走cparam那条路让四种小波的频率区间都落在 [fs/(2·totalscal), fs/2]。还有一个相邻的概念坑cmor1.5-2里的 1.5 是带宽参数、2 是中心频率顺序是cmorB-C。把它读成中心频率 1.5会导致后面算cparam时对不上尺度序列整体偏一截。写完参数名顺手print(pywt.central_frequency(cmor1.5-2))验一下返回值应该接近 2。5.3 观感失真图好看不等于信号对第五种情况时频图很干净但和原始波形对不上时间。现象是图上能量块的位置比实际冲击晚了一截或者两端莫名出现亮条。原因有两个叠加scipy.signal.stft默认boundaryzeros会在信号两端各补半个窗长的零时间轴被悄悄拉长同时默认返回的是复数谱直接np.abs画图没做对数压缩弱成分全被强成分盖住。解决方法是显式写boundaryNone, paddedFalse画图前做20*log10(amp/amp.max())并固定vmin/vmax。顺带一个显示层的坑pcolormesh(..., shadinggouraud)会对网格做插值平滑图看起来更漂亮但短时冲击那几帧会被抹圆冲击间隔这种关键信息反而模糊。排查阶段用shadingauto或nearest等参数定下来再考虑出图样式。中文标题也一样matplotlib.rc(font, familyMicrosoft YaHei)在缺少该字体的机器上会退化成方块和警告换成系统里实际存在的字体名或者干脆用英文标题。6. 取舍判据先于图像参数先于结论回过头看这份时频分析里那几个被选中的数字——STFT 窗长 32、CWT 尺度 128、小波cmor1.5-2——它们都不是看起来最好挑出来的而是各自被一个约束条件逼出来的。窗长 32 是时间窗必须短于 6 到 10 毫秒的冲击周期同时频点不能少到没法成像totalscal 128 是频率下限要压到 50 Hz 以下、又不愿意为低频多付一倍计算cmor1.5-2是复小波消除负频泄漏、带宽可调保证结论可迁移这两条同时成立。每条都对应一个能在纸上算出来的不等式而不是一句效果更好。这也是我看这类分析文档时最先找的东西作者有没有把判据写出来。判据在参数就能复现换一份工况数据也能照着重推一遍判据不在参数就成了一组玄学数字换台设备立刻失效。有一个边界值得说清楚。这份分析用的是 STFT 加 CWT 的时频图路线它能回答冲击多久来一次、能量集中在哪个时间段但回答不了特征频率到底是多少赫兹。原因在前面算过窗长 32 的频率分辨率 375 Hz比内圈特征频率 162 Hz 还大谱线上根本立不住。真要定位特征频率走的是包络解调那条路——Hilbert 变换取包络、去直流、再做 FFT在包络谱上找 162 Hz、107 Hz、141 Hz 这几根谱线及其转频边带。时频图和包络谱是互补的两件事别指望一张图全包。我早年做这类分析时踩过一次典型的坑为了让时频图更好看把窗长一路加到 256图确实细腻了故障之间的差别也没了我还以为是数据集本身的问题翻来覆去换小波、换尺度折腾了两天才想起来回去看时间分辨率那一列数字。从那以后我每次调时频参数都强制先跑一遍Δf、Δt、冲击周期这三个数的对照表写进注释里再动手画图。参数对不对纸上就能判不用等图出来。这份资料把数据集读取、STFT 窗长扫描、四种小波横向对比、totalscal 尺度对比这几段串成了一条完整链路代码和参数判断都在里面拿过去对着自己的数据重跑一遍那条链路上的每个取舍点都能亲手验一次。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/31 18:18:26 · 阅读 1.2k 轴承故障诊断中的STFT窗长与CWT尺度参数选择
简介这份文档围绕凯斯西储大学CWRU轴承故障数据集展开信号时频分析面向从事旋转机械故障诊断、健康监测与振动信号处理的学习者和研究人员帮助其理解如何从实测振动数据中提取可用于故障辨识的时频特征。内容先说明实验台组成与驱动端、风扇端、基座三处加速度计数据的差异再以正常信号和0.021英寸内圈、滚珠、外圈故障信号为对象比较短时傅里叶变换在重叠比例0.5、尺度16/32/64下的时频图表现并给出选择尺度32的依据连续小波变换部分则对比morl、cmor1-1、cmor1.5-2、cgau8等小波函数在尺度128下筛选出cmor1.5-2并进一步讨论32、64、128、256尺度对高低频特征的取舍。文档还附有Python读取mat数据、绘制原始信号与STFT时频图的代码片段便于对照复现和调整参数。资源包为1个docx文件大小约1.05MB已有1584人学习适合作为故障诊断入门与实验参数选择的参考材料。1. 一张时频图能分辨什么轴承故障数据集里被忽略的 85 毫秒很多人拿到某高校公开的轴承故障数据集第一反应是把0_0.mat和21_1.mat读进来跑一行stft撒一把 jet 色图然后在四个子图之间找差别。真正决定这张图能不能分辨内圈、滚珠、外圈的不是用 STFT 还是 CWT而是窗口长度与故障冲击节拍之间的比例。数据集里一个 0.021 英寸约 0.53 毫米的缺陷在 30 Hz 上下的轴转频下会激出周期仅几毫秒的冲击串窗口一旦长过几个冲击周期冲击就被平均成一条平滑亮带故障之间的差别随之消失。反过来窗口短到能咬住每一次冲击频率轴就粗到看不清特征频率。这篇把窗长、尺度、小波函数这三层参数各自在开什么价讲清楚也讲清那份时频分析文档里的参数是怎么被一步步逼出来的。2. STFT 的窗长为什么分辨率这个词在这份数据集上会骗人2.1 短时傅里叶变换在切的到底是什么短时傅里叶变换的做法可以一句话说完拿一个中心位于时间 α 的窗 g(t−α) 在信号上滑动在每个位置只对窗内那一小段做傅里叶变换。窗长固定窗的形状固定于是用多长的尺子去量信号这件事被一次性钉死了。钉死的代价就是不确定性原理时域窗越窄频域主瓣越宽两者乘积有下界。落到代码里这个下界表现为两个可以直接算出来的量——频率分辨率 Δf fs / nperseg和时间分辨率 Δt nperseg / fs。它们不是两个独立可调的旋钮而是同一根杠杆的两端动一个必然动另一个。这份数据集的驱动端加速度计采样率为 12 kHz单次分析截取 1024 点约 85.3 毫秒。85 毫秒里发生了什么值得先算清楚轴转频约 30 Hz一个转频周期 33 毫秒1024 点大约覆盖 2.5 个转频周期内圈故障特征频率的经验系数是 5.4152 倍转频约 162 Hz对应冲击重复周期 6.2 毫秒85 毫秒里大约能装下 13 到 14 次冲击外圈系数 3.5848 倍约 107 Hz周期 9.3 毫秒滚珠系数 4.7135 倍约 141 Hz周期 7.1 毫秒。这三个周期都落在 6 到 10 毫秒这个区间而窗长的候选值从 2.67 毫秒一路跨到 21.3 毫秒正好横跨能分辨冲击和糊掉冲击的分界线。2.2 窗长 16、32、64、256 的账Δf、Δt 与冲击节拍把四个候选窗长的账摆在一起很多结论就不用猜了。npersegΔf fs/npersegΔt nperseg/fs帧移重叠 50%频点个数06 kHz16750 Hz1.33 ms0.67 ms932375 Hz2.67 ms1.33 ms1764187.5 Hz5.33 ms2.67 ms3325646.875 Hz21.33 ms10.67 ms129窗长 16 时只剩 9 个频点整张图几乎是空的谈不上分辨率高低是信息量不够。窗长 256 的频率分辨率最好46.875 Hz理论上能把 162 Hz 和 107 Hz 的谱线分开——但它的时间窗是 21.33 毫秒一个窗里塞进 3 个以上的内圈冲击周期冲击串被卷积成一团时频图上只剩一条持续亮带四种状态长得越来越像。窗长 32 的处境刚好相反375 Hz 的频率分辨率连 162 Hz 与 107 Hz 都分不开可它的时间窗 2.67 毫秒比任何一个故障冲击周期都短每一次冲击都能被单独框住图上呈现的是一串等间隔的竖条间隔 6.2 毫秒还是 9.3 毫秒一眼可辨。注意这里有一个容易想反的地方。窗长 32 之所以在辨识度上胜出靠的不是频率分辨率而是它抓住的是时域冲击节拍不是频域谱线位置。想看清特征频率本身靠 STFT 图本来就不现实得走包络谱那条路。我一般会这么理解这个过程STFT 在这份数据集上实际充当的是一个冲击节奏计它回答每隔多久被敲一下至于敲在哪个频率上那是另一套工具的事。把这个前提想清楚参数的取舍就不再是玄学而是一道有明确约束的算术题。2.3 重叠 0.5 与帧移为什么不是 0.9重叠比例决定帧移帧移决定时间轴上的采样密度。重叠 0.5 意味着帧移等于窗长的一半窗长 32 时帧移 1.33 毫秒1024 点能得到约 64 帧。这个密度够用冲击周期 6.2 毫秒每 1.33 毫秒一帧一个冲击周期里有 4 到 5 帧冲击的起落能被描出来。重叠拉到 0.9帧数涨到 5 倍图更平滑但每一帧用的还是同样长度的窗信息量并没有增加多出来的只是插值密度计算量却翻了五倍。更麻烦的是重叠过高会让相邻帧高度相关时频图上出现横向的拖影看着连续实际是重复计数的假象。重叠降到 0帧移等于窗长冲击很可能正好落在两帧的交界处被两个窗各切一半幅度塌下去一半反而看不清。重叠 0.5 是这份分析的合理落点它让帧移小于最短的冲击周期又不至于制造冗余。import numpy as np from scipy.io import loadmat from scipy.signal import stft FS 12000 # 驱动端加速度计采样率 N 1024 # 单次分析窗长1024 点 ≈ 85.3 ms def load_de(path, key): CWRU 风格的 mat 文件读出来是 dict键名随文件变化必须显式指定 return loadmat(path)[key].reshape(-1)[:N] x load_de(21_1.mat, X209_DE_time) # 0.021 英寸内圈 for nperseg in (16, 32, 64, 256): noverlap nperseg // 2 # 重叠 50% → 帧移 nperseg/2 f, t, Z stft(x, fsFS, npersegnperseg, noverlapnoverlap, windowhann, boundaryNone, paddedFalse) df FS / nperseg # 频率分辨率 dt nperseg / FS # 时间分辨率单帧覆盖时长 print(fnperseg{nperseg:3d} df{df:7.2f} Hz dt{dt*1e3:6.2f} ms bins{len(f)})这段代码只做一件事把窗长扫描的代价量化出来。nperseg是每次傅里叶变换的点数直接决定dfnoverlap只影响时间轴上帧的密度不影响df和dtboundaryNone, paddedFalse关掉 scipy 默认的两端补零避免时间轴被悄悄延长半个窗。跑完打印出来的四行数字就是上面那张表的来源也是判断某个窗长到底行不行的第一手依据不用靠看图猜。3. 连续小波变换cparam 2·fc·totalscal 这行为什么决定了尺度可比性3.1 尺度-平移坐标与时间-频率坐标的换算连续小波变换把 STFT 里固定窗长、滑动位置的两参数结构换成了伸缩尺度、平移位置。小波函数 ψ(t) 被拉伸成 ψ((t−b)/a)a 是尺度b 是平移。尺度 a 变大的时候小波在时域被拉长、在频域被压窄并整体往低频挪于是一个参数同时改变了时间窗和频率窗——STFT 里那根固定的杠杆在这里变成了随频率自适应移动的杠杆。代价是坐标变了。小波变换天然活在尺度-平移平面上画图时要换回频率-时间平面换算关系是 f fc / (a · dt)fc 是小波的中心频率dt 是采样间隔。这个式子看起来平淡却是后面所有坑的源头尺度 a 本身没有物理意义只有乘上 fc 和 dt 之后才有赫兹。尺度与频率成反比这件事直接推出小波变换的双尺度性质。要关注高频特征选小尺度小波被压短时间定位准要关注低频特征选大尺度小波被拉长频率定位准。原文里那句尺度长度越大对低频特征关注越多说的正是这个方向性但如果只停在这句定性描述上调参还是靠手感。3.2 Morlet 的时频盒cmor 为什么能给出平滑包络Morlet 小波是一个复指数乘高斯包络写成 ψ(t) exp(−t²/2)·exp(jω₀t)。高斯包络保证时域局部化复指数保证频域是单峰。这两条性质合起来使它成为振动信号分析里最常用的小波之一——振动是准周期的冲击串既需要时间定位也需要频率定位。复数这个细节值得单独说。实小波比如实 Morlet、墨西哥帽的频谱关于零频对称正负频率各占一半对实信号做变换后取模正负频率的响应会叠在一起互相干涉幅值图上出现细密的振荡纹路。复小波是解析小波频谱只在正频率一侧有支撑负频部分几乎为零取模得到的是干净的包络不会自干涉。这就是为什么同样做 CWTcmor系列的幅值图看起来比morl干净冲击的位置和强度更直观。PyWavelets 里复数 Morlet 的参数写成cmorB-CB 是带宽参数C 是中心频率。cmor1.5-2意味着带宽 1.5、中心频率 2cmor1-1意味着带宽 1、中心频率 1。带宽参数控制高斯包络在频域的宽窄B 越大频域越宽、时域越短频率选择性下降但时间定位变好中心频率 C 决定小波的主频位置。cgau8是 8 阶复数高斯小波它是高斯函数取 8 阶导数没有单独的带宽参数阶数越高频率选择性越好、时域支撑越长。3.3 cparam 2·fc·totalscal 这行代码在做什么原文的尺度序列是这样生成的cparam 2 * fc * totalscal scales cparam / np.arange(totalscal, 0, -1)np.arange(totalscal, 0, -1)产生从 totalscal 递减到 1 的整数序列除法之后 scales 从大到小排列最小值是 cparam/totalscal 2·fc最大值是 cparam 2·fc·totalscal。代入频率换算 f fc/(a·dt)最小尺度 2·fc 对应最高频率 f_max fc / (2·fc·dt) 1/(2·dt) fs/2正好是奈奎斯特频率最大尺度 2·fc·totalscal 对应最低频率 f_min fc / (2·fc·totalscal·dt) fs / (2·totalscal)。两个结果里 fc 都被约掉了。这意味着无论选cmor1-1、cmor1.5-2还是cgau8只要各自的cparam用同一个2*fc*totalscal算频率覆盖区间就完全相同恒为 [fs/(2·totalscal), fs/2]。12 kHz 采样、totalscal128 时区间是 [46.875, 6000] Hz。import numpy as np, pywt FS, totalscal 12000, 128 dt 1.0 / FS for name in (cgau8, morl, cmor1-1, cmor1.5-2): fc pywt.central_frequency(name) # 数值积分得到的中心频率 cparam 2 * fc * totalscal scales cparam / np.arange(totalscal, 0, -1) freqs fc / (scales * dt) # 尺度 → 频率 df np.abs(np.diff(freqs)) print(f{name:11s} fc{fc:6.3f} scales[{scales.min():7.2f},{scales.max():7.2f}] ff[{freqs.min():7.2f},{freqs.max():7.2f}] Hz f相邻尺度频率差: 低端 {df.min():.3f} Hz / 高端 {df.max():.1f} Hz)跑一遍会看到两件事。第一四种小波打印出的f[46.88, 6000.00]完全一致fc列各不相同但都被cparam吸收掉了这正是这行代码存在的目的把中心频率这个自由参数从频率覆盖范围里剔除让换小波只改变时频盒形状不改变分析的频段比较才公平。第二df在低频端只有零点几赫兹在高频端上千赫兹——尺度是等间隔的频率是 1/a 关系所以频点在低频密集、高频稀疏。这是 CWT 的固有特性也是它相比 STFT 更适合看低频故障特征的原因。提示scales必须保持从大到小对应的freqs才会从小到大单调递增。contourf的 y 轴依赖这个单调性顺序反了图会上下颠倒低频跑到顶部看起来和别人的结果不一样。4. 小波函数与尺度实测cgau8、morl、cmor1-1、cmor1.5-2 差在哪4.1 复小波、实小波与辨识度的关系原文在尺度固定为 128 的前提下横向比了四种小波结论是cgau8和cmor1.5-2对不同故障的辨识度较高。这个结论背后的机制可以从三个角度拆。第一个是负频泄漏。morl在 PyWavelets 里是实 Morlet频谱对称取模后有自干涉纹路冲击之间的低能量区被填上振荡条纹三种复小波没有这个问题冲击与冲击之间的背景更干净对比度更高。第二个是时频盒形状。cgau8因为是 8 阶高斯导数时域支撑长、频域尖对窄带成分敏感冲击串出现时能量集中得很紧cmor1.5-2的带宽 1.5 比cmor1-1的 1 更宽时域更短对每次冲击的起落描得更细代价是频率方向的涂抹多一点。在这份数据上冲击节拍是主要判别信息时间定位更值钱所以cmor1.5-2的宽带宽反而成了优点。第三个是参数一致性。cgau8没有可调的带宽参数它的时频盒形状是被阶数锁死的cmor系列可以通过 B 和 C 调形状这意味着cmor的结论能往别的工况迁移cgau8的结论更依赖当前数据的冲击形态。综合下来选cmor1.5-2不是因为它在某张图上最漂亮而是因为它同时满足复小波无负频泄漏和时频盒可调这两条可迁移的性质。4.2 totalscal 从 32 到 256尺度长度变大到底改变了什么原文比较了 totalscal 取 32、64、128、256 时滚珠故障的时频图观察到尺度越长对低频特征关注越多。用上一节的公式可以把这个观察量化f_min fs/(2·totalscal)f_max 恒为 fs/2。totalscalf_minf_max低频端相邻尺度频率差尺度序列长度32187.5 Hz6000 Hz约 0.36 Hz326493.75 Hz6000 Hz约 0.18 Hz6412846.875 Hz6000 Hz约 0.09 Hz12825623.4 Hz6000 Hz约 0.05 Hz256两件事同时发生频率下限往下探同一频段的划分变密。totalscal 从 128 加到 256下限从 46.875 Hz 降到 23.4 Hz180 Hz 以下的低频区从十几个尺度变成三十几个尺度图上低频块的分层明显更细。但高频端一点也不受益6000 Hz 附近的尺度间隔几乎没变还是稀疏的。这解释了图上底部低频的区别最明显。同时也提醒一件事totalscal 不是越大越好。总尺度数翻倍计算量翻倍而新增的分辨率全部堆在低频。如果关心的频段在 1 kHz 以上把 totalscal 加到 256 只是白白多跑一倍时间。4.3 从 loadmat 到 contourf 的完整链路import numpy as np, pywt import matplotlib.pyplot as plt from scipy.io import loadmat FS, N 12000, 1024 dt 1.0 / FS def cwt_map(path, key, wname, totalscal): x loadmat(path)[key].reshape(-1)[:N] fc pywt.central_frequency(wname) cparam 2 * fc * totalscal scales cparam / np.arange(totalscal, 0, -1) # 大到小 → 频率低到高 coef, freqs pywt.cwt(x, scales, wname, sampling_perioddt) amp np.abs(coef) # 复小波取模 → 包络 db 20 * np.log10(amp / amp.max() 1e-12) # 归一化 dB压掉量级差 return freqs, db fig, axes plt.subplots(1, 4, figsize(20, 5), shareyTrue) for ax, ts in zip(axes, (32, 64, 128, 256)): freqs, db cwt_map(21_2.mat, X222_DE_time, cmor1.5-2, ts) t np.arange(N) * dt edge N // 8 # 掐掉两端边界效应 ax.contourf(t[edge:-edge], freqs, db[:, edge:-edge], levels30, cmapjet, vmin-40, vmax0) ax.set_title(ftotalscal{ts}) ax.set_xlabel(t / s) axes[0].set_ylabel(f / Hz) plt.show()几个参数的作用需要交代清楚。sampling_perioddt是pywt.cwt换算频率轴必需的漏掉它默认按 1 秒采样处理频率轴整体错 12000 倍。np.abs(coef)取模是针对复小波的操作实小波取模意义不大。20*log10归一化很关键——原始幅值的动态范围常常跨越三四个数量级线性色标下只有最强的那几个点有颜色其余全被压成背景色看起来没有特征实际是显示问题。vmin-40, vmax0把显示窗口固定住四张图才能横向比较否则每张图各自归一化亮度的相对关系就丢了。edge N // 8裁掉两端因为边界处小波支撑越界会出现两条亮的竖边不裁掉会误以为是故障冲击。5. 避坑与排查五类把时频图做废的现场5.1 坐标轴类时间轴和频率轴最容易一起错第一种情况时间轴长度对不上或者刻度的物理含义完全错。现象是横轴标着秒但从 0 到 1024看不出任何冲击间隔。原因是把窗长当成了采样率——原文代码里sampling_rate 1024紧接着t np.linspace(0, 1.0/sampling_rate, sampling_rate)这个 1024 是单次分析的点数不是每秒的采样点数。12 kHz 采样下 1024 点只覆盖 85.3 毫秒时间轴应该是 0 到 0.0853 秒。解决方法是单独定义一个FS 12000时间轴用np.arange(N) / FS别让窗长和采样率共用一个变量名。第二种情况CWT 图上下颠倒低频跑到顶部。现象是和别人同样的参数图却像镜像。原因是scales没有按从大到小排列pywt.cwt返回的freqs随之从大到小contourf按数组顺序铺 y 轴图就翻了。解决方法是坚持用cparam / np.arange(totalscal, 0, -1)生成尺度画图前打印一下freqs[:3]和freqs[-3:]确认单调递增。5.2 参数可比性类换个小波尺度就不是原来的尺度第三种情况跑morl那一组时直接抛NameError: name fc2 is not defined。现象是前三组正常第四组或者第二组一运行就崩。原因是fc2 pywt.central_frequency(wavename2)那一行被注释掉了但下一行cparam2 2 * fc2 * totalscal还在用这个变量。照抄代码时最容易踩这个坑而且报错位置和真正的问题不在同一行。解决方法是把注释行放开或者在循环里统一用fc pywt.central_frequency(name)现算从结构上消灭fc1、fc2、fc3、fc4这种容易错位的命名。第四种情况直接把scales np.arange(1, 129)喂给pywt.cwt然后拿不同小波的结果横向比。现象是每种小波的图频率范围都不一样cmor1.5-2的图明显比cmor1-1更高没法比。原因是没有用cparam 2*fc*totalscal做归一化f fc/(a·dt)里的fc没有被吸收中心频率 2 的小波频率整体高一倍。解决方法是统一走cparam那条路让四种小波的频率区间都落在 [fs/(2·totalscal), fs/2]。还有一个相邻的概念坑cmor1.5-2里的 1.5 是带宽参数、2 是中心频率顺序是cmorB-C。把它读成中心频率 1.5会导致后面算cparam时对不上尺度序列整体偏一截。写完参数名顺手print(pywt.central_frequency(cmor1.5-2))验一下返回值应该接近 2。5.3 观感失真图好看不等于信号对第五种情况时频图很干净但和原始波形对不上时间。现象是图上能量块的位置比实际冲击晚了一截或者两端莫名出现亮条。原因有两个叠加scipy.signal.stft默认boundaryzeros会在信号两端各补半个窗长的零时间轴被悄悄拉长同时默认返回的是复数谱直接np.abs画图没做对数压缩弱成分全被强成分盖住。解决方法是显式写boundaryNone, paddedFalse画图前做20*log10(amp/amp.max())并固定vmin/vmax。顺带一个显示层的坑pcolormesh(..., shadinggouraud)会对网格做插值平滑图看起来更漂亮但短时冲击那几帧会被抹圆冲击间隔这种关键信息反而模糊。排查阶段用shadingauto或nearest等参数定下来再考虑出图样式。中文标题也一样matplotlib.rc(font, familyMicrosoft YaHei)在缺少该字体的机器上会退化成方块和警告换成系统里实际存在的字体名或者干脆用英文标题。6. 取舍判据先于图像参数先于结论回过头看这份时频分析里那几个被选中的数字——STFT 窗长 32、CWT 尺度 128、小波cmor1.5-2——它们都不是看起来最好挑出来的而是各自被一个约束条件逼出来的。窗长 32 是时间窗必须短于 6 到 10 毫秒的冲击周期同时频点不能少到没法成像totalscal 128 是频率下限要压到 50 Hz 以下、又不愿意为低频多付一倍计算cmor1.5-2是复小波消除负频泄漏、带宽可调保证结论可迁移这两条同时成立。每条都对应一个能在纸上算出来的不等式而不是一句效果更好。这也是我看这类分析文档时最先找的东西作者有没有把判据写出来。判据在参数就能复现换一份工况数据也能照着重推一遍判据不在参数就成了一组玄学数字换台设备立刻失效。有一个边界值得说清楚。这份分析用的是 STFT 加 CWT 的时频图路线它能回答冲击多久来一次、能量集中在哪个时间段但回答不了特征频率到底是多少赫兹。原因在前面算过窗长 32 的频率分辨率 375 Hz比内圈特征频率 162 Hz 还大谱线上根本立不住。真要定位特征频率走的是包络解调那条路——Hilbert 变换取包络、去直流、再做 FFT在包络谱上找 162 Hz、107 Hz、141 Hz 这几根谱线及其转频边带。时频图和包络谱是互补的两件事别指望一张图全包。我早年做这类分析时踩过一次典型的坑为了让时频图更好看把窗长一路加到 256图确实细腻了故障之间的差别也没了我还以为是数据集本身的问题翻来覆去换小波、换尺度折腾了两天才想起来回去看时间分辨率那一列数字。从那以后我每次调时频参数都强制先跑一遍Δf、Δt、冲击周期这三个数的对照表写进注释里再动手画图。参数对不对纸上就能判不用等图出来。这份资料把数据集读取、STFT 窗长扫描、四种小波横向对比、totalscal 尺度对比这几段串成了一条完整链路代码和参数判断都在里面拿过去对着自己的数据重跑一遍那条链路上的每个取舍点都能亲手验一次。希望帮到你。本文还有配套的精品资源点击获取