简介：这份资源面向机械故障诊断、信号处理方向的研究生与工程技术人员，围绕凯斯西储大学（CWRU）轴承故障数据集展开时频分析实践。内容涵盖数据集实验台构成与驱动端、风扇端、基座三路振动信号的物理含义，并系统对比短时傅里叶变换与连续小波变换两类方法：STFT部分以正常信号与0.021英寸内圈、滚珠、外圈故障信号为对象，在0.5重叠比例下比较16、32、64三种尺度，最终确定尺度32；CWT部分则对比morl、cmor1-1、cmor1.5-2、cgau8等小波函数，选定cmor1.5-2并进一步比较32至256尺度下的时频表现。资源包为1个docx文档，约1.05MB，内含完整分析流程与可运行代码片段，便于读者复现实验、理解时频分辨率权衡并迁移到自身故障分类任务。目前已有1582人学习，适合作为入门时频分析与轴承故障识别的实操参考。
1. 凯斯西储大学轴承故障数据集做时频分析：为什么你画的频谱图总像一团糊
如果你从凯斯西储大学（CWRU）轴承故障数据集里随便抽一段振动信号，直接做 FFT，大概率会得到一张让你怀疑人生的图：故障特征频率的边带被淹没在噪声里，1x、2x 转频和它们的调制成分糊成一片。这不是你代码写错了，而是因为轴承故障产生的冲击是典型的非平稳瞬态信号——冲击的重复频率、衰减包络、共振频带都随时间变化，而 FFT 把整个时间轴上的信息平均掉了，时间分辨率直接归零。
时频分析要解决的就是这个问题：把一维振动信号映射到时间-频率二维平面，让你同时看到「什么时候发生了冲击」和「冲击激起了哪个频带」。短时傅里叶变换（STFT）和连续小波变换（CWT）是两条最常用的路径，前者适合看整体节奏，后者适合抓瞬态冲击的细节。这套流程适合做旋转机械故障诊断的工程师、研究生，以及想把 CWRU 数据集用起来的算法开发者。读完你应该能自己跑通从数据加载、时频变换、参数调优到故障特征提取的完整链路，并且知道每一步哪里容易翻车。
2. 先搞清楚 CWRU 数据集的信号长什么样，再谈时频分析
2.1 采样率、故障频率与数据文件命名规则
CWRU 数据集的核心文件是.mat格式，每个文件包含振动信号和对应的转速信息。驱动端加速度计数据采样率有 12 kHz 和 48 kHz 两档，风扇端只有 12 kHz。文件命名规则通常是「转速_故障类型_故障尺寸」，比如X105_DE_time表示 105 号实验、驱动端（Drive End）振动信号。故障类型包括内圈（IR）、外圈（OR）、滚动体（B）三种，故障尺寸从 0.007 英寸到 0.040 英寸不等。
做时频分析之前，你必须先算清楚几个特征频率，否则后面看到时频图上的亮线你也不知道对应什么。轴承型号是 6205-2RS JEM SKF，节径约 39.04 mm，滚动体直径约 7.94 mm，接触角 0 度，滚动体数量 9 个。外圈故障特征频率 BPFO、内圈故障特征频率 BPFI、滚动体故障特征频率 BSF 都有标准公式，代入转速就能算出来。
import scipy.io as sio import numpy as np # 加载 CWRU 数据文件，常见做法是直接读 .mat data = sio.loadmat('X105_DE_time.mat') # 不同版本的文件 key 名可能不同，先看 keys print(data.keys()) # 假设 key 是 'X105_DE_time' signal = data['X105_DE_time'].flatten() fs = 12000 # 驱动端 12 kHz 采样 N = len(signal) t = np.arange(N) / fs print(f"信号长度: {N}, 采样率: {fs} Hz, 时长: {N/fs:.2f} s")这段代码做了三件事：加载.mat文件、展平成一维数组、生成时间轴。参数说明：fs必须和数据采集时一致，用错采样率会导致所有频率轴刻度偏移，这是新手最常踩的坑之一。flatten()是因为.mat读出来通常是二维数组，不展平后面做变换会报维度错误。
2.2 为什么不能直接 FFT：从稳态假设到非平稳现实
FFT 的数学前提是信号在观测窗口内是平稳的，也就是说频率成分不随时间变化。但轴承故障的冲击信号完全违背这个前提：每次滚动体碾过缺陷，产生一个短促的冲击，冲击激发轴承座和传感器的共振，然后快速衰减。这个过程的持续时间可能只有几毫秒，而重复周期取决于转速和故障位置。
如果你对 10 秒数据做一次 FFT，得到的是所有冲击的平均频谱。故障特征频率的幅值会被背景噪声和其他频率成分稀释，尤其是早期微弱故障，故障冲击的能量可能比转频谐波低一个数量级。时频分析的价值就在于把「平均」换成「分段」，让你看到冲击发生的时刻和它激起的频带。
常见做法是先做一个简单的包络谱分析作为对照：对信号做带通滤波，然后 Hilbert 变换取包络，再对包络做 FFT。如果包络谱里能看到 BPFO 或 BPFI 的谐波，说明故障特征存在，接下来用 STFT 或 CWT 去定位它。
3. 短时傅里叶变换：窗长选不对，时频图白做
3.1 STFT 的窗函数、重叠率与频率分辨率三角关系
STFT 的思路很简单：把长信号切成很多短段，每段做 FFT，然后把结果按时间排列成二维矩阵。但这里有一个绕不开的矛盾——时间分辨率和频率分辨率不能同时最优。窗长越长，频率分辨率越高，但时间定位越模糊；窗长越短，时间定位越准，但频率分辨率下降。
对于 CWRU 驱动端 12 kHz 采样数据，故障冲击的持续时间通常在 1-5 ms 量级。如果你用 1024 点窗，对应约 85 ms，时间分辨率太差，冲击会被抹平。用 128 点窗，对应约 10.7 ms，时间分辨率够了，但频率分辨率只有约 94 Hz，对于 BPFO 在 100 Hz 左右的低频故障特征来说又不够细。
我的经验是：先确定你关心的频率范围。如果只看低频故障特征频率（通常 100-500 Hz），窗长取 512-1024 点；如果看高频共振带（2-5 kHz），窗长取 128-256 点。重叠率一般取 75%-90%，重叠越高，时频图越平滑，但计算量线性增长。
from scipy.signal import stft import matplotlib.pyplot as plt # STFT 参数 nperseg = 256 # 窗长 noverlap = 192 # 重叠 75% nfft = 512 # FFT 点数，补零到 512 提高频率轴密度 f, t_stft, Zxx = stft(signal, fs=fs, window='hann', nperseg=nperseg, noverlap=noverlap, nfft=nfft, boundary=None) # 转成 dB 并画图 mag_db = 20 * np.log10(np.abs(Zxx) + 1e-12) plt.pcolormesh(t_stft, f, mag_db, shading='gouraud', cmap='jet') plt.ylabel('Frequency (Hz)') plt.xlabel('Time (s)') plt.title('STFT Spectrogram') plt.colorbar(label='Magnitude (dB)') plt.show()参数说明：nperseg是每段长度，直接决定时间-频率分辨率的权衡；noverlap是段间重叠点数，75% 重叠意味着每段移动nperseg - noverlap个点；nfft是 FFT 点数，补零不增加真实分辨率，但让频率轴更密，画图更好看；boundary=None避免在两端补零造成虚假边缘效应。1e-12是防止 log10(0) 报错。
3.2 用 STFT 定位内圈故障的冲击周期
内圈故障的特点是：故障点随轴旋转，当它进入承载区时冲击更强，离开承载区时冲击减弱，所以时频图上会出现幅度调制。调制频率等于转频。如果你在时频图上看到一串等间隔的亮斑，间隔对应 BPFI 的倒数，而且亮斑幅度有周期性起伏，基本可以确认内圈故障。
实际操作时，建议先把时频图转成灰度图，然后沿时间轴对某个高频带（比如 2-4 kHz）做能量积分，得到一条「冲击包络曲线」。对这条曲线做 FFT，就能提取出故障特征频率。这比直接对原始信号做 FFT 的信噪比高得多。
# 取高频带 2000-4000 Hz 的能量积分 freq_mask = (f >= 2000) & (f <= 4000) energy_band = np.sum(np.abs(Zxx[freq_mask, :])**2, axis=0) # 对能量曲线做 FFT 找故障特征频率 from scipy.fft import fft, fftfreq N_e = len(energy_band) yf = np.abs(fft(energy_band - np.mean(energy_band))) xf = fftfreq(N_e, d=(t_stft[1] - t_stft[0])) # 只看正频率 pos_mask = xf > 0 plt.plot(xf[pos_mask], yf[pos_mask]) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.xlim(0, 500) plt.show()这段代码的逻辑是：STFT 已经帮你把高频共振带的能量随时间的变化提取出来了，对这条曲线做 FFT，相当于做了一次「包络谱分析」，但不需要手动设计带通滤波器。参数说明：频带范围 2000-4000 Hz 不是固定的，你需要根据时频图上能量集中的区域来调整。如果共振带在 3-5 kHz，就改这个范围。d=(t_stft[1] - t_stft[0])是 STFT 时间步长，必须用这个而不是1/fs，因为 STFT 输出已经降采样了。
4. 连续小波变换：抓瞬态冲击比 STFT 更顺手
4.1 小波基选择与尺度到频率的映射
CWT 和 STFT 的根本区别在于：STFT 用固定窗长，CWT 用可变窗长——高频处窗短，低频处窗长。这个特性天然适合轴承故障信号，因为冲击的瞬态部分在高频，需要短窗定位；而故障特征频率在低频，需要长窗分辨。
但 CWT 的参数比 STFT 更玄学。第一个要选的是小波基。对于冲击类信号，常用的有 Morlet 小波、Mexican hat 小波和 Daubechies 系列。Morlet 小波在时频聚集性上表现最好，也是 CWRU 相关文献里用得最多的。第二个要选的是尺度范围，尺度a和频率f的对应关系是f ≈ fc * fs / a，其中fc是小波中心频率。Morlet 小波的fc大约是 0.8125 Hz（归一化后）。
import pywt # 连续小波变换 scales = np.arange(1, 128) # 尺度范围 wavelet = 'cmor1.5-1.0' # 复 Morlet 小波，带宽 1.5，中心频率 1.0 coefficients, frequencies = pywt.cwt(signal, scales, wavelet, sampling_period=1/fs) # 画时频图 plt.pcolormesh(t, frequencies, np.abs(coefficients), shading='gouraud', cmap='jet') plt.ylabel('Frequency (Hz)') plt.xlabel('Time (s)') plt.title('CWT Scalogram') plt.colorbar() plt.show()参数说明：scales决定了分析的频率范围，尺度越小对应频率越高。cmor1.5-1.0是复 Morlet 小波，1.5是带宽参数，1.0是中心频率。sampling_period=1/fs让pywt.cwt直接返回实际频率轴，省去手动换算。注意pywt.cwt的输出是复数，取模得到幅值。
4.2 用 CWT 系数做故障特征增强的实操步骤
CWT 输出的系数矩阵和 STFT 类似，但频率轴不是线性的，低频处分辨率高，高频处分辨率低。这个特性对故障诊断有利：故障特征频率通常在低频，CWT 在低频的精细分辨能力正好用上。
具体操作分三步。第一步，对 CWT 系数取模，得到幅值矩阵。第二步，选择包含故障特征频率的尺度范围，对这个范围内的系数沿尺度轴求和，得到一条时间序列。第三步，对这条时间序列做 FFT，提取故障特征频率。
# 第一步：取模 cwt_mag = np.abs(coefficients) # 第二步：选择低频尺度范围，对应频率 100-500 Hz freq_mask_cwt = (frequencies >= 100) & (frequencies <= 500) cwt_band = np.sum(cwt_mag[freq_mask_cwt, :], axis=0) # 第三步：对时间序列做 FFT N_c = len(cwt_band) yf_cwt = np.abs(fft(cwt_band - np.mean(cwt_band))) xf_cwt = fftfreq(N_c, d=1/fs) pos_mask_cwt = xf_cwt > 0 plt.plot(xf_cwt[pos_mask_cwt], yf_cwt[pos_mask_cwt]) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.xlim(0, 500) plt.show()这段代码和 STFT 那段的逻辑完全一致，只是把 STFT 的频带换成了 CWT 的尺度带。关键区别在于：CWT 在低频段的频率分辨率更高，所以 100-500 Hz 范围内的故障特征频率会更清晰。但代价是计算量比 STFT 大，尤其是尺度范围取很宽的时候。
提示：CWT 的尺度范围不要盲目取太大。尺度 1 到 128 已经覆盖了 12 kHz 采样下的大部分有用频带。尺度超过 256 后，频率低于 50 Hz，转频和故障特征频率可能混在一起，反而不好分辨。
5. 避坑与排查：时频分析里那些让你白干一天的坑
5.1 采样率用错导致频率轴整体偏移
现象：时频图上看到的亮线频率和理论计算的 BPFO、BPFI 对不上，偏差比例固定。原因：CWRU 数据集里驱动端有 12 kHz 和 48 kHz 两种采样率，风扇端只有 12 kHz。如果你加载的是 48 kHz 文件但按 12 kHz 处理，所有频率会变成实际值的 1/4。解决：加载数据后先检查文件来源，确认采样率。如果不确定，可以看信号长度和实验时长的对应关系，或者直接查数据集说明文档。
5.2 STFT 窗长选得太长导致冲击被平均掉
现象：时频图上看不到明显的垂直亮线，故障冲击的瞬态特征完全消失，整张图看起来像稳态信号的频谱。原因：窗长过大，时间分辨率太低，每个窗内包含了多个冲击周期，FFT 把冲击平均成了稳态成分。解决：先估算冲击持续时间，窗长取冲击持续时间的 2-3 倍。对于 CWRU 数据，驱动端 12 kHz 采样下，窗长 128-256 点是比较稳妥的起点。
5.3 CWT 尺度范围没覆盖故障特征频率
现象：CWT 时频图在低频区域一片模糊，看不到故障特征频率的亮线。原因：尺度范围上限太小，对应的最低频率高于故障特征频率。比如尺度只取到 64，对应最低频率约 150 Hz，而 BPFO 可能只有 100 Hz。解决：先算清楚故障特征频率，然后根据f ≈ fc * fs / a反推需要的尺度上限。保险做法是尺度上限取到fc * fs / f_min，其中f_min是你关心的最低频率。
5.4 忘记去均值导致时频图出现 0 Hz 亮带
现象：时频图最底部有一条极亮的水平线，掩盖了低频故障特征。原因：振动信号通常有直流偏置，STFT 和 CWT 都会把这个直流分量映射到 0 Hz 附近。解决：做时频变换前先减去信号均值。一行代码的事，但不做的话低频分析基本废掉。
5.5 用错时间轴导致故障特征频率提取偏差
现象：对 STFT 能量曲线做 FFT 时，提取出的故障特征频率和理论值有固定比例偏差。原因：STFT 输出已经是降采样后的时间序列，时间步长是nperseg - noverlap除以fs，而不是1/fs。如果你用1/fs做 FFT 的频率轴，所有频率会偏大。解决：用t_stft[1] - t_stft[0]作为 FFT 的时间步长，这是最稳妥的做法。
6. 把时频图变成可量化的故障指标：我的三个私藏技巧
时频图好看归好看，但如果你只是画出来看一眼就扔掉，那等于白做。真正有价值的是把时频图变成可量化的指标，用来做故障分类或退化趋势跟踪。我一般用三个技巧。
第一个技巧：时频熵。对 STFT 或 CWT 的幅值矩阵做归一化，然后计算 Shannon 熵。正常轴承的时频分布比较集中，熵值低；故障轴承的时频分布更分散，熵值高。这个指标对早期故障很敏感，而且不需要知道故障特征频率的具体数值。
def time_freq_entropy(mag_matrix): # 归一化 p = mag_matrix / np.sum(mag_matrix) # 避免 log(0) p = p[p > 0] # Shannon 熵 entropy = -np.sum(p * np.log2(p)) return entropy # 对 STFT 幅值矩阵计算时频熵 entropy_stft = time_freq_entropy(np.abs(Zxx)) print(f"STFT 时频熵: {entropy_stft:.4f}")第二个技巧：频带能量比。选择两个频带，一个包含故障特征频率，一个作为参考（比如转频附近），计算能量比值。这个比值随时间的变化曲线可以反映故障的严重程度。正常状态下比值稳定，故障发展时比值上升。
第三个技巧：时频图模板匹配。如果你有已知故障类型的数据，可以把它们的平均时频图作为模板，然后用相关系数匹配新数据。这个方法在故障类型分类上比直接看时频图靠谱得多，因为人眼容易被颜色映射和对比度误导。
from scipy.stats import pearsonr # 假设 template 是已知故障的平均时频图 # new_mag 是新数据的时频图幅值 template_flat = template.flatten() new_flat = new_mag.flatten() corr, _ = pearsonr(template_flat, new_flat) print(f"与模板的相关系数: {corr:.4f}")这三个技巧的共同点是：把二维时频图压缩成一维或标量指标，方便做统计分析和自动化决策。我自己的习惯是，拿到一批 CWRU 数据后，先算时频熵做快速筛查，再用频带能量比做趋势跟踪，最后用模板匹配做分类确认。这套组合拳打下来，比单纯画图看效率高一个数量级。
希望帮到你。
本文还有配套的精品资源，点击获取
