简介面向机械故障诊断初学者与信号处理方向研究者的时频分析学习文档围绕凯斯西储大学轴承故障数据集展开适合理解正常、内圈、滚珠与外圈故障振动信号的时频特征差异。文档结合驱动端、风扇端和基座加速度计数据介绍短时傅里叶变换与连续小波变换的基本流程讨论分段加窗、重叠比例、尺度参数及小波函数选择对时间分辨率和频率分辨率的影响并通过不同尺度时频图直观对比故障辨识度给出尺度32的STFT、cmor1.5-2小波等参数取舍思路。文中附有可参考的代码片段便于复现关键实验与理解参数设置。资源包仅含1个docx文件约1.05MB便于集中阅读。已有1584人学习可作为轴承故障诊断入门、课程实验或论文前期调研的参考资料。1. 从一个「换参数就变脸」的时频图说起这套方案要解决什么拿到一份轴承振动数据集多数人的第一反应是把波形画出来随手调个窗长跑一遍 STFT图一出来就觉得活儿干完了。真正被问住的地方在后面把窗长从 32 改成 64整张时频图的纹理全变了内圈和滚珠之间的差别也跟着糊掉那到底哪一个才算数这套方案要解决的就是这个问题——在一份公开的轴承故障数据集上把时频分析每一步参数为什么这么定讲清楚让短时傅里叶变换与连续小波两条路线各自有可复现的取值依据而不是靠肉眼看图反复试。它适合正在做故障诊断方向课程设计、毕业设计的人也适合要给设备健康监测出一份分析方案、需要把参数标定过程写进文档的从业者。2. 数据与实验台拆解三个测点、四种状态与一段 1024 点窗动手写代码之前先把数据的来路摸清楚。时频分析最怕的不是算法写错而是拿着不同工况、不同通道的信号硬凑成一组对比图图看着挺漂亮结论一问就塌。这套数据集的价值恰恰在于它的实验条件被记录得比较完整故障类型、故障尺寸、测点位置、负载转速都能对上号只要你愿意在动手前多花二十分钟读一遍说明。2.1 实验台构成与三类加速度计测点的物理含义数据采集实验台由四部分组成一台 1.5 kW2 马力的驱动电机、一段连接用的扭矩传感器与译码器、一套功率测试计以及一台未在台架上直接露面的电子控制器。电机带动轴系旋转扭矩传感器夹在中间既测扭矩也提供转速参考功率测试计吸收负载控制器负责把转速稳定在设定值。轴承就装在这条轴系上人为地在内圈、滚珠、外圈上加工出不同尺寸的点蚀然后让电机带着它跑。真正决定你能分析出什么的是加速度计贴在哪。台架上共有三个测点驱动端、风扇端、基座。驱动端的信号主要受转子和传动系统激励轴承故障的冲击成分在这里最明显风扇端受叶片旋转与风扇系统激励主导适合看叶片失衡一类的问题基座测的是整机振动反映的是不平衡、底座松动这类全局性状态。三个测点测的是同一台机器的不同物理侧面不是三份互为备份的数据。提示做轴承故障分类时优先用驱动端数据。风扇端和基座里的轴承冲击成分被结构传递路径衰减得厉害同样的故障在图上可能只差一点点初学者很容易误判成「这个方法不灵」。2.2 四种状态的文件命名与 MAT 键名对照本轮分析用到四组信号一组正常三组是 0.021 英寸的故障分别对应内圈、滚珠、外圈。数据以 MAT 格式存放每组一个文件文件里又按通道存了不同的变量名。这里最容易翻车的地方是变量名不统一——正常数据的键是X097_DE_time而三组故障数据的键分别是X209_DE_time、X222_DE_time、X234_DE_time编号完全不连续。我一般会在方案文档里先把这张对照表列出来后面所有代码都引用它避免脚本里散落着一堆魔法字符串状态文件名MAT 键名故障尺寸说明正常0_0.matX097_DE_time无基准对照内圈故障21_1.matX209_DE_time0.021 英寸驱动端通道滚珠故障21_2.matX222_DE_time0.021 英寸驱动端通道外圈故障21_3.matX234_DE_time0.021 英寸驱动端通道loadmat读出来是一个字典除了数据键还塞了一堆__header__、__version__之类的元信息所以拿到字典后一定要先打印一遍键名再取值。用.reshape(-1)把列向量压平成一维也是必要的否则后续stft会按二维数组处理返回的形状和你想的完全不一样。2.3 1024 点窗与后续分析的输入约定原始信号长度动辄十几万点直接丢进时频变换既慢又没必要。工程上的通行做法是切一段固定长度的窗把窗内的信号当作一次观测。这里选 1024 点一是够覆盖若干个完整的旋转周期二是四种状态用同一个长度切片时频图的横纵坐标刻度一致肉眼对比才成立——如果正常信号切 1024 点、故障信号切 2048 点两张图的时间轴对不上任何「看起来不一样」的结论都没有说服力。采样频率同样要统一。驱动端通道的采样率按数据说明文件标注的取值来这套数据常见的是 12 kHz具体以你手上的说明为准不要图省事把窗长 1024 直接当成采样率填进去那样画出来的频率轴刻度全是错的答辩时被问一句「这个 3000 Hz 是怎么来的」就露馅了。统一约定之后后面 STFT 和小波两条路线共用同一份切片才谈得上参数对比。3. STFT 参数标定窗长 16/32/64 对比与 0.5 重叠的实现短时傅里叶变换的思路很朴素让一个时间窗在信号上滑动窗滑到哪儿就在哪儿做一次傅里叶变换把每次变换的频谱按时间顺序摞起来就得到一张时频图。它之所以能同时给出时间和频率信息靠的就是这个窗——窗把信号切成了局部片段傅里叶变换只看得见窗内那一段。代价也随之而来窗开得宽频率看得细但时间糊窗开得窄时间看得准但频率糊。这不是实现问题是这类变换的固有约束。3.1 窗口滑动模型与时间-频率分辨率的取舍把窗函数记作 g(t−α)α 就是窗中心的位置。窗在时间轴上平移每平移一步就在窗限定的范围内做一次傅里叶变换于是原本只在一维上展开的信号被摊到了时间-频率二维平面上。可分辨的最小时间间隔大致跟着窗长走可分辨的最小频率间隔大致与窗长成反比两者乘起来是个常数你没法同时把两边都压小。理解这一点之后参数选择就变成了一个明确的工程问题你要看的是冲击的瞬间位置还是故障特征频率附近的谱线轴承的局部点蚀在时域上表现为周期性的冲击脉冲冲击本身很窄但在频域上会激起一串以特征频率为间隔的边带。想看清冲击发生在哪一刻窗要短想看清边带间隔窗要长。方案里把 16、32、64 三档都跑一遍就是为了把这条权衡曲线画出来给人看而不是拍脑袋定一个数。3.2 读取与原始波形核对读数据这一步要先把四种状态的原始波形摊在一张图上核对确认故障信号里确实能看到周期性的冲击、正常信号相对平缓。这一步不做后面时频图出了异常你根本分不清是参数问题还是数据本身就没读对。import numpy as np from scipy.io import loadmat import matplotlib.pyplot as plt # 中文字体避免标题变成方块 plt.rcParams[font.family] Microsoft YaHei plt.rcParams[axes.unicode_minus] False FS 12000 # 驱动端采样频率按数据说明文件标注取值 WIN 1024 # 统一切片长度四种状态保持一致 def load_de(path, key): 读取 MAT 文件中的驱动端通道压平并截取前 WIN 点 mat loadmat(path) return mat[key].reshape(-1)[:WIN] sig { 正常: load_de(0_0.mat, X097_DE_time), 内圈: load_de(21_1.mat, X209_DE_time), 滚珠: load_de(21_2.mat, X222_DE_time), 外圈: load_de(21_3.mat, X234_DE_time), } fig, axes plt.subplots(2, 2, figsize(16, 8)) for ax, (name, x) in zip(axes.ravel(), sig.items()): ax.plot(np.arange(WIN) / FS, x, lw0.8) ax.set_title(f{name} 原始波形) ax.set_xlabel(时间 / s) ax.set_ylabel(加速度 / g) plt.tight_layout() plt.show()这段代码做了三件事统一读取四个文件、统一截取 1024 点、统一用np.arange(WIN) / FS生成时间轴。时间轴除以采样频率这一步看着不起眼但它决定了横坐标是「第几个采样点」还是「第几秒」写进报告里前者是不专业的。load_de把「文件名 键名」封装成一次调用后面换数据只需要改字典不用动绘图逻辑。3.3 四组窗长的对比实验与 32 的定标理由对比实验固定重叠比例为 0.5只让窗长变化。重叠的作用是让相邻两窗之间有一半的数据被重复利用滑动的步长变成窗长的一半时频图在时间方向上不会出现明显的条带断裂。窗长 16 时步长只有 8 点图很密但每一列的频率分辨率都很粗窗长 128 时步长 64 点频率方向细腻了时间方向开始糊成一片。from scipy.signal import stft wins [16, 32, 64, 128] overlap 0.5 # 重叠比例 x sig[内圈] # 先用内圈数据做窗长对比 fig, axes plt.subplots(2, 2, figsize(16, 9)) for ax, w in zip(axes.ravel(), wins): noverlap int(w * overlap) # 重叠样本数 f, t, Z stft(x, fsFS, npersegw, noverlapnoverlap) ax.pcolormesh(t, f, np.abs(Z), shadinggouraud, cmapjet) ax.set_title(f窗长 {w}重叠 {noverlap} 点) ax.set_xlabel(时间 / s) ax.set_ylabel(频率 / Hz) plt.tight_layout() plt.show()nperseg是每段长度noverlap是相邻段的重复点数两者之差就是滑动步长。重叠比例取 0.5 意味着步长恒为窗长的一半这样四张图只有分辨率这一个变量在动对比才有意义。np.abs(Z)取的是复数频谱的模pcolormesh把它映射成颜色。跑完四张图会发现一个规律窗长 16 时间轴最细冲击的起跳点看得清清楚楚但纵向的频率条纹糊成一团窗长 64 频率方向的边带结构开始分得开时间上冲击的边界已经有点晕开。再往下比较四种状态在同一窗长下的辨识度——16 时四种状态的图差别不大128 时故障之间开始互相长得像32 和 64 之间是辨识度最好的区间。综合两端方案里把窗长定在 32时间分辨率足以框住冲击出现的时刻频率分辨率也还留得住故障特征频率附近的谱结构。这个「定在 32」不是算出来的最优点而是在可解释的区间里选的一个折中报告里必须把 16 和 64 的对比图一并放上说清楚为什么两端都被排除。只放一张 32 的图等于把选择过程藏起来了。3.4 时频图渲染细节pcolormesh 与 gouraud渲染方式看着是美化问题实际会直接影响你从图上读出来的东西。shadinggouraud会对相邻网格做插值图看起来平滑但弱能量区域容易被插值抹掉本来能看见的一条细边带就消失了。做定量分析时我一般改用默认的 flat 渲染或者干脆对幅值取对数再画把动态范围压下来低频弱成分才浮得出来。另一个细节是配色。jet对比强烈适合放在答辩 PPT 里让人一眼看到差异但它的亮度不是单调变化的看图的人容易把某个中间色误当成峰值。如果图要进正式文档换成viridis这类亮度单调的配色更稳妥。这些都属于出图规范不影响算法结论但会影响别人对你结论的信任度。4. 连续小波变换选型四种小波函数对比与 cmor1.5-2 尺度序列构造STFT 的窗长一旦定死整张图上所有时刻、所有频率用的都是同一个分辨率。轴承信号恰恰是「低频段要细看、高频段要粗看」的典型故障特征频率落在低频冲击激起的高频共振带又需要足够的时间精度去定位。这种需求天然适合小波——它用可伸缩的窗低频处窗自动拉长、高频处窗自动缩短不用你在两者之间二选一。4.1 从 STFT 到 CWT为什么要换一套基底连续小波变换把信号与一族经过伸缩和平移的小波函数做内积。尺度参数控制小波的伸缩尺度大小波被拉宽对上的是低频成分尺度小小波被压窄对上的是高频成分。平移参数控制小波在时间轴上的位置扫一遍就得到整条时间轴上的系数。这和 STFT 的差别不只是「更灵活」。STFT 的基是固定的复指数CWT 的基可以换——换成 Morlet 系、高斯系时频局部化的性格完全不同。对一个具体的故障信号选哪个小波本身就是方案的一部分需要拿出对比证据。方案里准备了四个候选morl、cmor1-1、cmor1.5-2、cgau8尺度统一设为 128 先做函数对比把变量控制住。morl是实 Morlet 小波时频域局部化都不错对振动信号这类瞬态成分友好cmor是复 Morletcmor1-1和cmor1.5-2里的前一个数字是中心频率、后一个数字是带宽带宽越大频率方向越散、时间方向越紧cgau8是 8 阶复高斯小波对高斯型包络的近似能力强冲击成分在它下面往往显得格外干净。4.2 尺度序列怎么算cparam 与 arange 的构造小波变换里最容易写错的一行不是pywt.cwt调用而是尺度序列的构造。直接把np.arange(1, 129)丢进去得到的频率轴和真实频率对不上因为尺度与实际频率之间的换算还差一个中心频率的因子。通行做法是先取小波的中心频率fc用它和一个常数项把尺度序列归一化import pywt totalscal 128 # 尺度序列长度即频率方向的分辨档数 wavename cmor1.5-2 fc pywt.central_frequency(wavename) # 小波的中心频率 cparam 2 * fc * totalscal # 尺度序列的常数项 scales cparam / np.arange(totalscal, 0, -1) # 由大到小排列 print(fc, cparam, scales[:3], scales[-3:])cparam 2 * fc * totalscal这个式子的作用是让尺度序列与采样频率解耦2 * fc把中心频率折算成周期乘上总尺度数之后最大尺度对应到的频率大致落在奈奎斯特频率附近最小尺度落到最低关注频率。np.arange(totalscal, 0, -1)生成的是从大到小的整数序列除下来得到的scales就是由大到小——尺度大对应低频排在上面画出来的时频图低频在下方符合直觉。注意np.arange的方向写反改成np.arange(1, totalscal1)不会报错但整张时频图会上下颠倒低频跑到顶部。这种错误在图上不容易一眼看出来因为图依旧花花绿绿很好看。4.3 四种小波函数的对比实验固定totalscal 128让四个小波函数各跑一遍用的还是内圈数据。这样四张图只有小波基在变其余条件完全一致。wavelets [cgau8, morl, cmor1-1, cmor1.5-2] x sig[内圈] fig, axes plt.subplots(2, 2, figsize(16, 9)) for ax, wn in zip(axes.ravel(), wavelets): fc pywt.central_frequency(wn) scales (2 * fc * 128) / np.arange(128, 0, -1) coef, freqs pywt.cwt(x, scales, wn, sampling_period1.0 / FS) ax.contourf(np.arange(WIN) / FS, freqs, np.abs(coef), levels32, cmapjet) ax.set_title(f{wn}尺度长度 128) ax.set_xlabel(时间 / s) ax.set_ylabel(频率 / Hz) plt.tight_layout() plt.show()sampling_period传的是采样周期的倒数关系——这里填1.0 / FS返回的freqs才是真实的物理频率。这个参数漏填pywt会按 1 处理频率轴整体错位。np.abs(coef)取小波系数模值contourf用 32 层等值线填充比imshow更容易控制动态范围。四张图摆在一起看cgau8和cmor1.5-2对三种故障的区分度明显更好内圈、滚珠、外圈的能量带位置和形态能拉开差距而morl和cmor1-1下三种故障的图长得比较接近尤其是滚珠和外圈容易混。再权衡一次时间与频率的分配——cgau8时间方向更利落但频率方向的谱线略显发散cmor1.5-2两边的表现更均衡最终定为后续分析的主用小波。4.4 尺度 32/64/128/256 的对比与低频关注小波函数定下来之后尺度长度成了唯一的变量。把这四个档位各跑一遍重点看每张图底部低频区域的区别totals [32, 64, 128, 256] wn cmor1.5-2 fc pywt.central_frequency(wn) fig, axes plt.subplots(2, 2, figsize(16, 9)) for ax, n in zip(axes.ravel(), totals): scales (2 * fc * n) / np.arange(n, 0, -1) # 常数项随总尺度同步变化 coef, freqs pywt.cwt(x, scales, wn, sampling_period1.0 / FS) ax.contourf(np.arange(WIN) / FS, freqs, np.abs(coef), levels32, cmapjet) ax.set_title(ftotalscal{n}) ax.set_xlabel(时间 / s) ax.set_ylabel(频率 / Hz) plt.tight_layout() plt.show()这里有个容易被忽略的一致性要求cparam里的totalscal必须和尺度序列长度用同一个数改成 32 就把常数项一起改。如果常数项还用 128 的值、序列只取 32 个点频率轴的范围会跟着缩水四张图的纵坐标量程都不一样横向对比直接失效。结果印证了小波变换的双尺度性质尺度长度越大对低频特征的刻画越充分四张图的差别集中在下半部分——256 时低频区域的能量分布明显比 32 时更细腻故障之间的低频差异被摊开了。反过来如果关注的是冲击瞬间的高频共振应该往小尺度走。方案里最终取 128 作为默认档理由是它在低频细节和时间定位之间留有余量同时计算量还没有到 256 那档的负担。5. 避坑与排查MAT 键名、尺度序列、归一化这几处最容易翻车上面两套流程跑通一次不难难的是换个人、换台机器、隔两周再跑还能得到同一张图。下面这几处是我在这类时频分析里见过翻车最多的基本都发生在「代码没报错但结果是错的」这个区间。5.1 数据读取环节的两个坑现象一脚本跑一半抛KeyError提示某个X***_DE_time不存在。原因通常是照抄了别人脚本里的键名而不同文件的编号根本不连续正常数据是 097故障数据跳到 209、222、234没有规律可循。解决办法是别硬编码先打印mat.keys()看一眼或者写一个按后缀匹配的取键函数[k for k in mat if k.endswith(DE_time)]取到唯一一个再往下走。如果返回多个说明这个文件里同时存了驱动端和风扇端得指定位置再筛。现象二能读出数据画出来的波形却是一条几乎平的线或者满屏削顶。原因是loadmat返回的是二维列向量直接丢进stft会被当成多通道信号逐列处理输出的形状和预期完全不同另一种情况是误取了别的通道比如把基座数据当成驱动端。解决办法是先reshape(-1)压平再打印shape和前十个值确认量级正常振动的加速度值应该在正负零点几到几个 g 之间如果全是几千上万多半取错了变量。5.2 参数设置环节的三个坑现象三NameError: name sampling_period is not defined。这个几乎每个人都踩过一次——pywt.cwt的采样周期参数在某些写法里必须显式传而代码里那句sampling_period是从别处抄来的变量名自己脚本里压根没定义。解决办法是老老实实写sampling_period1.0/FS别依赖外部变量。现象四时频图上下颠倒。前文提过尺度序列的方向反了。判断方法是找一条已知的特征频率——比如外圈故障的特征频率如果它出现在图的上方而不是下方就是方向错了。改法是确保np.arange从大到小生成并用scales[0] scales[-1]做一次断言。现象五四张对比图全是一片深红什么结构都看不出来。原因是小波系数的动态范围跨了好几个数量级线性映射下弱成分全被压没了。解决办法是对模值取对数再画np.log10(np.abs(coef) 1e-12)或者用等值线分层填充并把层数调低让弱能量区域也能显示出来。加的那个极小量是为了防止取对数时出现负无穷。提示这五条里前两条会导致结果完全错误却不一定报错后三条会报错或者图明显异常。真正危险的是那些「能跑通、图也好看、结论却不对」的情况所以每改一次参数都拿正常信号跑一遍当基准它应该始终是相对平缓的图一旦正常信号都变得花里胡哨说明参数或者数据读取出了问题。6. 验收与答辩用一张参数对照表和一次冷启动重跑自证时频分析这类工作最容易在答辩现场被问倒的地方不是算法原理而是「你这些参数凭什么这么定」。我一般的做法是在报告最后放一张参数对照表把每个环节的取值、候选范围和定下来的依据并排列出来让人一眼看到选择过程而不是只看到一个结果分析环节关键参数候选范围最终取值定标依据切片窗内点数—1024覆盖多个旋转周期四种状态统一STFT窗口长度16 / 32 / 64 / 12832时间与频率分辨率的折中点STFT重叠比例—0.5步长恒定抑制时间向条带断裂CWT小波函数cgau8 / morl / cmor1-1 / cmor1.5-2cmor1.5-2三种故障区分度与两向均衡性CWT尺度长度32 / 64 / 128 / 256128低频细节与计算量的平衡CWT尺度序列常数—2·fc·totalscal尺度与物理频率对齐这张表还有一个隐性作用它逼着你把「为什么不是另一个」讲清楚。答辩时最忌讳的答法是「试出来的这个看着最好」换成「16 和 128 分别在时间、频率方向失效32 是两端都能接受的区间」问题的性质就从主观变成了可讨论。第二步是冷启动重跑。把脚本里所有硬编码的路径、变量名清掉从头到尾跑一遍中间不修任何东西。这一步专治「笔记本里能跑、换台机器就崩」——环境里残留的变量、IDE 缓存、上一次调试时改过没存盘的参数都会在这时候暴露。跑通之后把四张图导出来存成固定文件名和报告里的图逐张对一遍确认不是手滑贴错了。第三步是留一个最小的复现入口。我通常在脚本开头放一个参数区把所有可调的量集中成几行常量答辩现场被要求「换成外圈看看」时只需要改一个键名不用翻遍全文找变量# 参数区答辩现场改这里就够了 CASE 外圈 # 正常 / 内圈 / 滚珠 / 外圈 METHOD cwt # stft / cwt WIN 1024 # 切片长度 STFT_WIN 32 # STFT 窗长 WAVELET cmor1.5-2 # CWT 小波函数 TOTALSCAL 128 # CWT 尺度长度 # data sig[CASE] if METHOD stft: f, t, Z stft(data, fsFS, npersegSTFT_WIN, noverlapint(STFT_WIN * 0.5)) coef, freqs, tt Z, f, t else: fc pywt.central_frequency(WAVELET) scales (2 * fc * TOTALSCAL) / np.arange(TOTALSCAL, 0, -1) coef, freqs pywt.cwt(data, scales, WAVELET, sampling_period1.0 / FS) tt np.arange(WIN) / FS这段结构的价值在于把「方法」和「数据」解耦成两个开关两种变换共用同一份切片和同一套时间轴出图函数只写一次。现场被追问某个参数时改一个常量、重跑一遍几十秒内就能给出新图比翻着代码解释有说服力得多。最后还有一件事把「结论的边界」写出来。这套方案能说明的是 0.021 英寸这一档故障在特定负载转速下的时频特征换成更大的故障尺寸或者别的转速最优尺度未必还是这两个数。承认这一点不会削弱方案反而会让评审觉得你知道自己在做什么。从那以后我每次交时频分析的结果都会强制走一遍「参数表写全、冷启动重跑、留一个改一个常量就能重出的入口」这三步被问住过太多次这算是唯一管用的后悔药。希望帮到你。本文还有配套的精品资源点击获取
