轴承故障诊断时频分析：STFT与连续小波参数标定实践

简介：面向机械故障诊断初学者与信号处理方向研究者的时频分析学习文档，围绕凯斯西储大学轴承故障数据集展开，适合理解正常、内圈、滚珠与外圈故障振动信号的时频特征差异。文档结合驱动端、风扇端和基座加速度计数据，介绍短时傅里叶变换与连续小波变换的基本流程，讨论分段加窗、重叠比例、尺度参数及小波函数选择对时间分辨率和频率分辨率的影响，并通过不同尺度时频图直观对比故障辨识度，给出尺度32的STFT、cmor1.5-2小波等参数取舍思路。文中附有可参考的代码片段，便于复现关键实验与理解参数设置。资源包仅含1个docx文件，约1.05MB，便于集中阅读。已有1584人学习，可作为轴承故障诊断入门、课程实验或论文前期调研的参考资料。

1. 从一个「换参数就变脸」的时频图说起：这套方案要解决什么

拿到一份轴承振动数据集，多数人的第一反应是把波形画出来，随手调个窗长跑一遍 STFT，图一出来就觉得活儿干完了。真正被问住的地方在后面：把窗长从 32 改成 64，整张时频图的纹理全变了，内圈和滚珠之间的差别也跟着糊掉，那到底哪一个才算数？这套方案要解决的就是这个问题——在一份公开的轴承故障数据集上，把时频分析每一步参数为什么这么定讲清楚，让短时傅里叶变换与连续小波两条路线各自有可复现的取值依据，而不是靠肉眼看图反复试。它适合正在做故障诊断方向课程设计、毕业设计的人，也适合要给设备健康监测出一份分析方案、需要把参数标定过程写进文档的从业者。

2. 数据与实验台拆解：三个测点、四种状态与一段 1024 点窗

动手写代码之前，先把数据的来路摸清楚。时频分析最怕的不是算法写错，而是拿着不同工况、不同通道的信号硬凑成一组对比图，图看着挺漂亮，结论一问就塌。这套数据集的价值恰恰在于它的实验条件被记录得比较完整：故障类型、故障尺寸、测点位置、负载转速都能对上号，只要你愿意在动手前多花二十分钟读一遍说明。

2.1 实验台构成与三类加速度计测点的物理含义

数据采集实验台由四部分组成：一台 1.5 kW（2 马力）的驱动电机、一段连接用的扭矩传感器与译码器、一套功率测试计，以及一台未在台架上直接露面的电子控制器。电机带动轴系旋转，扭矩传感器夹在中间既测扭矩也提供转速参考，功率测试计吸收负载，控制器负责把转速稳定在设定值。轴承就装在这条轴系上，人为地在内圈、滚珠、外圈上加工出不同尺寸的点蚀，然后让电机带着它跑。

真正决定你能分析出什么的是加速度计贴在哪。台架上共有三个测点：驱动端、风扇端、基座。驱动端的信号主要受转子和传动系统激励，轴承故障的冲击成分在这里最明显；风扇端受叶片旋转与风扇系统激励主导，适合看叶片失衡一类的问题；基座测的是整机振动，反映的是不平衡、底座松动这类全局性状态。三个测点测的是同一台机器的不同物理侧面，不是三份互为备份的数据。

提示：做轴承故障分类时，优先用驱动端数据。风扇端和基座里的轴承冲击成分被结构传递路径衰减得厉害，同样的故障在图上可能只差一点点，初学者很容易误判成「这个方法不灵」。

2.2 四种状态的文件命名与 MAT 键名对照

本轮分析用到四组信号：一组正常，三组是 0.021 英寸的故障，分别对应内圈、滚珠、外圈。数据以 MAT 格式存放，每组一个文件，文件里又按通道存了不同的变量名。这里最容易翻车的地方是变量名不统一——正常数据的键是X097_DE_time，而三组故障数据的键分别是X209_DE_time、X222_DE_time、X234_DE_time，编号完全不连续。

我一般会在方案文档里先把这张对照表列出来，后面所有代码都引用它，避免脚本里散落着一堆魔法字符串：

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, lw=0.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, fs=FS, nperseg=w, noverlap=noverlap) ax.pcolormesh(t, f, np.abs(Z), shading="gouraud", cmap="jet") 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

渲染方式看着是美化问题，实际会直接影响你从图上读出来的东西。shading='gouraud'会对相邻网格做插值，图看起来平滑，但弱能量区域容易被插值抹掉，本来能看见的一条细边带就消失了。做定量分析时我一般改用默认的 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是复 Morlet，cmor1-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, totalscal+1)）不会报错，但整张时频图会上下颠倒，低频跑到顶部。这种错误在图上不容易一眼看出来，因为图依旧花花绿绿很好看。

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_period=1.0 / FS) ax.contourf(np.arange(WIN) / FS, freqs, np.abs(coef), levels=32, cmap="jet") 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_period=1.0 / FS) ax.contourf(np.arange(WIN) / FS, freqs, np.abs(coef), levels=32, cmap="jet") ax.set_title(f"totalscal={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_period=1.0/FS，别依赖外部变量。

现象四：时频图上下颠倒。前文提过，尺度序列的方向反了。判断方法是找一条已知的特征频率——比如外圈故障的特征频率，如果它出现在图的上方而不是下方，就是方向错了。改法是确保np.arange从大到小生成，并用scales[0] > scales[-1]做一次断言。

现象五：四张对比图全是一片深红，什么结构都看不出来。原因是小波系数的动态范围跨了好几个数量级，线性映射下弱成分全被压没了。解决办法是对模值取对数再画（np.log10(np.abs(coef) + 1e-12)），或者用等值线分层填充并把层数调低，让弱能量区域也能显示出来。加的那个极小量是为了防止取对数时出现负无穷。

提示：这五条里，前两条会导致结果完全错误却不一定报错，后三条会报错或者图明显异常。真正危险的是那些「能跑通、图也好看、结论却不对」的情况，所以每改一次参数，都拿正常信号跑一遍当基准，它应该始终是相对平缓的图，一旦正常信号都变得花里胡哨，说明参数或者数据读取出了问题。

6. 验收与答辩：用一张参数对照表和一次冷启动重跑自证

时频分析这类工作，最容易在答辩现场被问倒的地方不是算法原理，而是「你这些参数凭什么这么定」。我一般的做法是在报告最后放一张参数对照表，把每个环节的取值、候选范围和定下来的依据并排列出来，让人一眼看到选择过程，而不是只看到一个结果：

这张表还有一个隐性作用：它逼着你把「为什么不是另一个」讲清楚。答辩时最忌讳的答法是「试出来的，这个看着最好」，换成「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, fs=FS, nperseg=STFT_WIN, noverlap=int(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_period=1.0 / FS) tt = np.arange(WIN) / FS

这段结构的价值在于把「方法」和「数据」解耦成两个开关，两种变换共用同一份切片和同一套时间轴，出图函数只写一次。现场被追问某个参数时，改一个常量、重跑一遍，几十秒内就能给出新图，比翻着代码解释有说服力得多。

最后还有一件事：把「结论的边界」写出来。这套方案能说明的是 0.021 英寸这一档故障在特定负载转速下的时频特征，换成更大的故障尺寸或者别的转速，最优尺度未必还是这两个数。承认这一点不会削弱方案，反而会让评审觉得你知道自己在做什么。从那以后我每次交时频分析的结果，都会强制走一遍「参数表写全、冷启动重跑、留一个改一个常量就能重出的入口」这三步，被问住过太多次，这算是唯一管用的后悔药。希望帮到你。

本文还有配套的精品资源，点击获取