简介这份资源面向机械故障诊断、信号处理方向的研究生与工程技术人员围绕凯斯西储大学CWRU轴承故障数据集展开时频分析实践。内容涵盖数据集实验台构成与驱动端、风扇端、基座三路振动信号的物理含义并系统对比短时傅里叶变换与连续小波变换两类方法STFT部分以正常信号与0.021英寸内圈、滚珠、外圈故障信号为对象在0.5重叠比例下比较16、32、64三种尺度最终确定尺度32CWT部分则对比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, fsfs, windowhann, npersegnperseg, noverlapnoverlap, nfftnfft, boundaryNone) # 转成 dB 并画图 mag_db 20 * np.log10(np.abs(Zxx) 1e-12) plt.pcolormesh(t_stft, f, mag_db, shadinggouraud, cmapjet) plt.ylabel(Frequency (Hz)) plt.xlabel(Time (s)) plt.title(STFT Spectrogram) plt.colorbar(labelMagnitude (dB)) plt.show()参数说明nperseg是每段长度直接决定时间-频率分辨率的权衡noverlap是段间重叠点数75% 重叠意味着每段移动nperseg - noverlap个点nfft是 FFT 点数补零不增加真实分辨率但让频率轴更密画图更好看boundaryNone避免在两端补零造成虚假边缘效应。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, axis0) # 对能量曲线做 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_period1/fs) # 画时频图 plt.pcolormesh(t, frequencies, np.abs(coefficients), shadinggouraud, cmapjet) 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_period1/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, :], axis0) # 第三步对时间序列做 FFT N_c len(cwt_band) yf_cwt np.abs(fft(cwt_band - np.mean(cwt_band))) xf_cwt fftfreq(N_c, d1/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(fSTFT 时频熵: {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 数据后先算时频熵做快速筛查再用频带能量比做趋势跟踪最后用模板匹配做分类确认。这套组合拳打下来比单纯画图看效率高一个数量级。希望帮到你。本文还有配套的精品资源点击获取
CWRU轴承故障时频分析：STFT与CWT实战指南
简介这份资源面向机械故障诊断、信号处理方向的研究生与工程技术人员围绕凯斯西储大学CWRU轴承故障数据集展开时频分析实践。内容涵盖数据集实验台构成与驱动端、风扇端、基座三路振动信号的物理含义并系统对比短时傅里叶变换与连续小波变换两类方法STFT部分以正常信号与0.021英寸内圈、滚珠、外圈故障信号为对象在0.5重叠比例下比较16、32、64三种尺度最终确定尺度32CWT部分则对比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, fsfs, windowhann, npersegnperseg, noverlapnoverlap, nfftnfft, boundaryNone) # 转成 dB 并画图 mag_db 20 * np.log10(np.abs(Zxx) 1e-12) plt.pcolormesh(t_stft, f, mag_db, shadinggouraud, cmapjet) plt.ylabel(Frequency (Hz)) plt.xlabel(Time (s)) plt.title(STFT Spectrogram) plt.colorbar(labelMagnitude (dB)) plt.show()参数说明nperseg是每段长度直接决定时间-频率分辨率的权衡noverlap是段间重叠点数75% 重叠意味着每段移动nperseg - noverlap个点nfft是 FFT 点数补零不增加真实分辨率但让频率轴更密画图更好看boundaryNone避免在两端补零造成虚假边缘效应。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, axis0) # 对能量曲线做 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_period1/fs) # 画时频图 plt.pcolormesh(t, frequencies, np.abs(coefficients), shadinggouraud, cmapjet) 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_period1/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, :], axis0) # 第三步对时间序列做 FFT N_c len(cwt_band) yf_cwt np.abs(fft(cwt_band - np.mean(cwt_band))) xf_cwt fftfreq(N_c, d1/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(fSTFT 时频熵: {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 数据后先算时频熵做快速筛查再用频带能量比做趋势跟踪最后用模板匹配做分类确认。这套组合拳打下来比单纯画图看效率高一个数量级。希望帮到你。本文还有配套的精品资源点击获取
