简介：这份资源围绕凯斯西储大学（CWRU）轴承故障数据集展开信号时频分析，面向从事机械设备故障诊断、健康监测与振动信号处理的研究生、工程师及科研人员，帮助读者理解如何从实测振动数据中提取故障特征。压缩包内仅含1个docx文档，约1.05MB，以图文与代码结合的方式组织内容。文档先介绍CWRU数据集的实验台构成，包括1.5KW电机、扭矩传感器与功率测试计，并区分驱动端、风扇端和基座三处加速度计信号所反映的物理特征；随后以正常信号与0.021英寸内圈、滚珠、外圈故障信号为对象，演示短时傅里叶变换在重叠比例0.5、尺度16/32/64下的时频图差异，说明尺度32的选取依据；再对比morl、cmor1-1、cmor1.5-2、cgau8等小波函数，最终确定cmor1.5-2并比较尺度32至256的连续小波变换效果。目前已有1582人学习，适合作为时频分析入门与故障分类特征工程的实操参考。
1. 凯斯西储大学轴承故障数据集：为什么时频分析是绕不开的第一道坎
如果你手头有一份凯斯西储大学（CWRU）轴承故障数据，打算做故障诊断，大概率会先画个时域波形看看。然后你会发现，正常轴承和故障轴承的时域波形，肉眼看上去差别并不大——尤其是内圈故障和滚动体故障，波形毛刺混在一起，光看幅值根本分不清谁是谁。这不是你数据读错了，而是时域信号本身把频率信息藏起来了。
凯斯西储大学轴承故障数据集之所以成为旋转机械故障诊断领域被反复使用的基准，核心原因是它的实验条件标注清晰：故障类型（内圈、外圈、滚动体）、故障直径（0.007、0.014、0.021 英寸）、负载（0～3 马力）都有明确记录。但标注清晰不等于信号好分。轴承故障的冲击成分是瞬态的、非平稳的，傅里叶变换把整个时段的频率成分平均掉之后，冲击发生的时刻信息就丢了。时频分析要解决的正是这个问题：把一维时间信号映射到时间-频率二维平面，让冲击在哪个时刻、哪个频带出现变得可见。
这篇内容面向的是已经拿到 CWRU 数据、准备做特征提取或故障分类的工程师和研究生。我会从数据加载开始，把短时傅里叶变换（STFT）和连续小波变换（CWT）两条最常用的时频分析路径走通，给出可复现的代码、参数设置逻辑，以及我在实际处理 CWRU 数据时踩过的坑。不涉及深度学习模型搭建，只聚焦时频分析这一步怎么做扎实。
2. 先把 CWRU 数据读对：文件结构、采样率与通道选择
2.1 CWRU 数据的目录结构和命名规则
CWRU 数据集在官网以 .mat 文件分发，命名规则通常是「故障类型 + 故障直径 + 负载」的组合。比如 IR007_0.mat 表示内圈故障、直径 0.007 英寸、负载 0 马力； B014_1.mat 表示滚动体故障、直径 0.014 英寸、负载 1 马力；正常基线数据一般命名为 Normal_0.mat 到 Normal_3.mat 。驱动端加速度计数据存在 DE_time 变量里，风扇端存在 FE_time ，基座存在 BA_time 。
做时频分析，我一般只用驱动端 DE_time 。原因很直接：驱动端离轴承座最近，故障冲击的信噪比最高，风扇端信号经过传递路径衰减后，冲击成分容易被淹没。如果你要做多通道融合，那是另一个话题，单通道时频分析从 DE 通道起步最稳。
import scipy.io as sio
import numpy as np
# 读取 CWRU 驱动端数据
def load_cwru_de(filepath):
"""
读取 CWRU .mat 文件中的驱动端加速度信号
filepath: .mat 文件路径
返回: 一维 numpy 数组 (DE_time)
"""
mat_data = sio.loadmat(filepath)
# 变量名通常是 DE_time，但不同版本可能带后缀
de_key = [k for k in mat_data.keys() if 'DE_time' in k]
if not de_key:
raise KeyError(f"未找到 DE_time 变量，现有键: {list(mat_data.keys())}")
signal = mat_data[de_key[0]].flatten()
return signal
signal = load_cwru_de('IR007_0.mat')
print(f"信号长度: {len(signal)}, 采样点数: {signal.shape}")
这段代码的关键在变量名匹配。CWRU 官网下载的 .mat 文件里，变量名有时是 X105_DE_time 这种带编号的形式，直接写 mat_data['DE_time'] 会报 KeyError。用列表推导式模糊匹配 DE_time 是更稳妥的做法。 flatten() 把列向量压成一维，后续做时频变换时不用再关心原始维度。
2.2 采样率确认与重采样决策
CWRU 驱动端数据的采样率是 12 kHz（12,000 Hz），部分文件是 48 kHz。这个信息不在 .mat 文件里，需要根据你下载的数据版本确认。采样率搞错的后果很严重：STFT 的频率轴会整体偏移，CWT 的尺度-频率对应关系也会错，最后得到的时频图看着有模有样，但频率刻度全是错的。
我一般会在代码里显式定义采样率，而不是从数据里猜：
FS = 12000 # CWRU 驱动端常用采样率，48kHz 数据需改为 48000
# 如果要做降采样，先做抗混叠滤波
from scipy.signal import decimate
def downsample_signal(signal, fs_original, fs_target):
"""
降采样，带抗混叠滤波
decimate 的 q 参数必须是整数
"""
q = fs_original // fs_target
if q < 2:
return signal, fs_original
signal_ds = decimate(signal, q, ftype='iir', zero_phase=True)
return signal_ds, fs_original // q
signal_ds, fs_new = downsample_signal(signal, FS, 3000)
print(f"降采样后采样率: {fs_new}, 点数: {len(signal_ds)}")
降采样不是必须的，但如果你后续要做 CWT 且尺度范围设得比较大，原始 12 kHz 数据计算量会很大。降到 3 kHz 通常足够覆盖轴承故障的特征频率——CWRU 轴承故障特征频率一般在几百赫兹到 2 kHz 之间，3 kHz 采样率满足奈奎斯特条件。 decimate 函数自带 IIR 抗混叠滤波，比直接切片 signal[::q] 安全得多，后者会把高频噪声混叠到低频，时频图上会出现虚假的低频能量带。
注意：降采样前一定确认目标采样率大于你关心的最高故障特征频率的两倍。CWRU 内圈故障特征频率在负载 0 时约 160 Hz 左右，但冲击的宽带成分可以延伸到 5 kHz 以上，做包络分析时降采样要谨慎。
3. 短时傅里叶变换：窗长、重叠与频率分辨率的三角权衡
3.1 STFT 在 CWRU 信号上的参数选择逻辑
STFT 的核心思想是用一个滑动窗把长信号切成短段，对每段做 FFT，再把结果按时间排列。窗长决定了时间分辨率和频率分辨率的平衡：窗越长，频率分辨率越高，但时间定位越模糊；窗越短，时间定位越准，但频率分辨率下降。这不是玄学，是海森堡不确定性原理在信号处理里的直接体现。
对 CWRU 轴承故障信号，我一般从窗长 256 点起步（12 kHz 采样率下约 21 ms），重叠 50%～75%。故障冲击的持续时间通常在几毫秒量级，256 点窗能覆盖一个完整的冲击衰减过程，同时频率分辨率约 47 Hz（12000/256），足以区分故障特征频率及其谐波。
from scipy.signal import stft
import matplotlib.pyplot as plt
def compute_stft(signal, fs, nperseg=256, noverlap=None):
"""
计算 STFT
nperseg: 窗长（采样点数）
noverlap: 重叠点数，默认 75% 重叠
"""
if noverlap is None:
noverlap = int(nperseg * 0.75)
f, t, Zxx = stft(signal, fs=fs, nperseg=nperseg,
noverlap=noverlap, window='hann',
detrend=False, return_onesided=True)
return f, t, Zxx
f, t, Zxx = compute_stft(signal, FS, nperseg=256)
magnitude = np.abs(Zxx)
plt.figure(figsize=(10, 4))
plt.pcolormesh(t, f, 20*np.log10(magnitude + 1e-10), shading='gouraud')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.title('STFT Magnitude (dB)')
plt.colorbar(label='dB')
plt.ylim(0, 3000)
plt.tight_layout()
plt.show()
window='hann' 是默认选择，汉宁窗的旁瓣衰减比矩形窗好，能减少频谱泄漏。 detrend=False 表示不做去趋势，因为 CWRU 信号已经做过预处理，再去趋势可能把故障冲击的低频成分也去掉。 return_onesided=True 对实数信号只返回正频率半边，省一半计算量。
3.2 从 STFT 图里读出故障特征：频率轴和实际特征频率的对应
STFT 图出来之后，关键是从图上找到故障特征频率的踪迹。CWRU 内圈故障特征频率（BPFI）的计算公式是：
BPFI = (n/2) × fr × (1 + (d/D) × cosθ)
其中 n 是滚动体个数，fr 是轴转频，d 是滚动体直径，D 是节径，θ 是接触角。CWRU 实验台的具体参数：n=9，d=0.3126 英寸，D=1.537 英寸，θ=0。负载 0 时轴转频约 29.95 Hz，算出来 BPFI 约 162 Hz。
在 STFT 图上，你应该能在 162 Hz 及其倍频（324 Hz、486 Hz…）附近看到能量集中。如果看不到，先检查采样率是否设对，再检查窗长是否太长导致频率分辨率不够——窗长 256 点时频率分辨率 47 Hz，162 Hz 和 324 Hz 能分开，但如果窗长只有 64 点，分辨率降到 187 Hz，这两个峰就糊在一起了。
# 提取 STFT 在特定频率切片上的时间演化
def extract_freq_slice(f, t, Zxx, target_freq, bandwidth=20):
"""
提取目标频率附近的时间-幅值曲线
target_freq: 目标频率 (Hz)
bandwidth: 频率带宽 (Hz)
"""
freq_mask = (f >= target_freq - bandwidth) & (f <= target_freq + bandwidth)
slice_mag = np.mean(np.abs(Zxx[freq_mask, :]), axis=0)
return t, slice_mag
t_slice, mag_slice = extract_freq_slice(f, t, Zxx, target_freq=162, bandwidth=20)
# mag_slice 的周期性峰值对应冲击重复频率
这段代码把 STFT 矩阵在频率轴上做切片，得到目标频率附近能量随时间的变化。如果轴承有内圈故障， mag_slice 会呈现周期性的峰值，峰峰间隔对应冲击重复周期。这个周期信息在时域波形上很难直接量出来，但在时频平面上变得直观。
3.3 STFT 的局限：窗长固定带来的分辨率困境
STFT 最大的问题是一次变换只能用一种窗长。轴承故障信号里，冲击成分是高频瞬态，需要短窗来定位；而故障特征频率的谐波结构是低频持续成分，需要长窗来分辨。固定窗长意味着你必须在两者之间做取舍。
我的经验是：如果目标是检测冲击是否存在，窗长取 128～256 点；如果目标是精确测量特征频率，窗长取 512～1024 点。不要试图用一个窗长解决所有问题，那只会两头不讨好。这也是为什么后面要引入 CWT——小波变换的多分辨率特性天然适合处理这种「高频短时、低频长时」的信号。
4. 连续小波变换：尺度参数怎么选、复小波为什么更适合轴承信号
4.1 CWT 的尺度-频率映射关系与参数初始化
CWT 用一组尺度可变的母小波去卷信号，尺度小对应高频、时间分辨率高，尺度大对应低频、频率分辨率高。这个特性正好匹配轴承故障信号的结构：冲击瞬间用小尺度捕捉，谐波结构用大尺度分辨。
在 Python 里做 CWT，常用 PyWavelets 库。母小波选复 Morlet 小波（ cmor ），因为轴承故障信号是调制信号，复小波的实部和虚部能同时给出幅值和相位信息，后续做包络分析更方便。
import pywt
def compute_cwt(signal, fs, scales=None, wavelet='cmor1.5-1.0'):
"""
计算连续小波变换
scales: 尺度数组，None 时自动生成
wavelet: 母小波名称，cmorB-C 中 B 是带宽参数，C 是中心频率
"""
if scales is None:
# 尺度范围对应频率范围：f = (center_freq * fs) / scale
# cmor1.5-1.0 的中心频率约 1.0 Hz（归一化）
scales = np.arange(1, 128)
coefficients, frequencies = pywt.cwt(signal, scales, wavelet,
sampling_period=1.0/fs)
return coefficients, frequencies
coeffs, freqs = compute_cwt(signal, FS, scales=np.arange(1, 256))
print(f"CWT 系数矩阵形状: {coeffs.shape}")
print(f"频率范围: {freqs[0]:.1f} Hz ~ {freqs[-1]:.1f} Hz")
cmor1.5-1.0 里的 1.5 是带宽参数， 1.0 是中心频率。带宽参数越大，小波在频率域的带宽越宽，时间分辨率越低。对轴承冲击信号，我一般用 cmor1.5-1.0 或 cmor1.0-1.0 ，前者频率聚集性稍好，后者时间定位更锐利。 scales 从 1 到 256 覆盖的频率范围取决于采样率：12 kHz 采样率下，scale=1 对应约 12 kHz，scale=256 对应约 47 Hz。如果你只关心 0～3 kHz，把 scales 上限设到 256 就够了，再大只会增加计算量。
4.2 用 CWT 系数矩阵做故障可视化的完整流程
CWT 系数矩阵是复数矩阵，取模之后得到幅值矩阵，画出来就是时频图。和 STFT 图相比，CWT 图在低频段频率分辨率更高，在高频段时间分辨率更高，冲击的起始和结束时刻更清晰。
import matplotlib.pyplot as plt
def plot_cwt_scalogram(coeffs, freqs, fs, title='CWT Scalogram'):
"""
绘制 CWT 尺度图
coeffs: CWT 系数矩阵 (len(scales), len(signal))
freqs: 对应的频率数组
"""
magnitude = np.abs(coeffs)
# 转成 dB 尺度，动态范围设 60 dB
mag_db = 20 * np.log10(magnitude + 1e-10)
vmax = np.max(mag_db)
vmin = vmax - 60
plt.figure(figsize=(10, 5))
plt.pcolormesh(np.arange(len(coeffs[0]))/fs, freqs,
mag_db, shading='gouraud',
vmin=vmin, vmax=vmax, cmap='jet')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.yscale('log')
plt.title(title)
plt.colorbar(label='Magnitude (dB)')
plt.ylim(10, fs/2)
plt.tight_layout()
plt.show()
plot_cwt_scalogram(coeffs, freqs, FS, title='CWT Scalogram - IR007')
vmin 和 vmax 的设置很关键。如果不设，matplotlib 会自动按数据最小最大值拉伸，噪声会被放大成和冲击一样的亮度，图上一片花。设 60 dB 动态范围是经验值：冲击能量比背景噪声高 40～60 dB，这个范围能把冲击突出来，同时压住噪声。 plt.yscale('log') 让低频段展开，因为轴承故障特征频率集中在低频，线性频率轴会把低频区域压得很扁。
4.3 CWT 和 STFT 在 CWRU 数据上的效果对比
同一段内圈故障信号，STFT 和 CWT 给出的信息侧重不同。STFT 图在 162 Hz 附近能看到周期性的能量团，但每个能量团的时间宽度被窗长模糊了；CWT 图在同样频率位置能看到更锐利的冲击时刻，而且冲击的衰减过程（从高频到低频的能量下降）更清楚。
这不是说 CWT 一定比 STFT 好。STFT 的计算量小，参数少，适合做实时监测的初步筛查；CWT 计算量大，参数多（母小波类型、带宽、中心频率、尺度范围），但多分辨率特性让它在分析非平稳信号时信息量更大。我的做法是：先用 STFT 快速扫一遍，确认故障特征频率的大致位置，再用 CWT 在目标频段做精细分析。
提示：CWT 的尺度数组不要用默认值。
pywt.cwt如果不传 scales，会用一个默认的尺度范围，那个范围不一定覆盖你关心的频率。显式传入np.arange(1, 256)或根据目标频率反算尺度，结果才可控。
5. 避坑与排查：CWRU 时频分析中五个让我返工的问题
5.1 现象：时频图上出现一条横贯整个时间轴的水平亮线
原因：这通常是工频干扰或传感器固有共振。CWRU 实验台在 60 Hz（美国电网频率）附近有干扰，如果时频图在 60 Hz 或 120 Hz 有一条不随时间变化的亮线，基本可以确定是电源干扰。另一个可能是加速度计的谐振频率，CWRU 使用的加速度计谐振频率在 20～30 kHz，但如果你降采样前没做抗混叠滤波，谐振能量会混叠到低频。
解决：先检查信号频谱，确认干扰频率。如果是 60 Hz 工频，用陷波滤波器滤掉；如果是传感器谐振，降采样前加抗混叠滤波。不要试图在时频图上「忽略」这条线，它会掩盖故障特征频率处的微弱能量。
5.2 现象：STFT 图上的频率刻度和理论计算的 BPFI 对不上
原因：采样率设错了。CWRU 数据有 12 kHz 和 48 kHz 两个版本，如果你下载的是 48 kHz 数据但代码里写 FS=12000 ，所有频率都会偏小 4 倍。另一个可能是 .mat 文件里的信号经过了降采样但你没注意到。
解决：用信号本身的特性反推采样率。CWRU 数据在 0 Hz 附近有明显的轴转频成分，负载 0 时约 29.95 Hz，负载 1 时约 29.5 Hz。在频谱上找到这个峰，用理论值反推采样率。如果对不上，检查数据来源。
5.3 现象：CWT 系数矩阵全是 NaN 或 Inf
原因： pywt.cwt 对信号中的 NaN 或 Inf 零容忍。CWRU 原始数据一般不会有 NaN，但如果你做了滤波或降采样，边界处理不当可能引入 NaN。另一个可能是信号幅值过大导致浮点溢出，但这种情况少见。
解决：在 CWT 之前加一行 signal = np.nan_to_num(signal, nan=0.0, posinf=0.0, neginf=0.0) 。同时检查滤波器的边界处理， scipy.signal.filtfilt 比 lfilter 更安全，因为它做了前后向滤波，相位不失真。
5.4 现象：不同负载下的时频图差异巨大，无法用同一套参数分析
原因：负载变化会改变轴转频，进而改变故障特征频率。负载 0 时轴转频约 29.95 Hz，负载 3 时约 28.5 Hz，BPFI 从 162 Hz 变到 154 Hz。如果你用固定频率切片去提取特征，负载 3 的数据会偏。
解决：要么按负载分别设置目标频率，要么做阶次分析（order analysis），把频率轴转成轴转频的倍数。阶次分析需要转速信号，CWRU 数据里没有直接的转速通道，但可以从轴转频成分估算。我一般对每个负载单独算 BPFI，不跨负载用同一套参数。
5.5 现象：时频图看着很漂亮，但分类器准确率上不去
原因：时频图是给人看的，分类器需要的是数值特征。直接把时频图当图像喂给 CNN 是一种做法，但如果你用的是传统机器学习，需要从时频矩阵里提取统计量：比如目标频带的能量均值、方差、峰值因子、小波系数的熵。光有图没有特征，分类器学不到东西。
解决：从 STFT 或 CWT 矩阵里提取频带能量比、时频熵、小波包能量等特征。我一般会在 BPFI、BPFO、BSF 及其谐波位置各取一个频带，算每个频带的能量占比，再加上时频图的整体熵，组成特征向量。这套特征在 CWRU 上的分类准确率通常能到 95% 以上，比直接用时域统计量高 10 个百分点左右。
6. 从时频图到特征向量：一个可复用的 CWRU 特征提取函数
时频分析本身不是终点，终点是拿到能喂给分类器的特征。我把自己反复用的一套特征提取流程整理成函数，输入是原始信号和采样率，输出是一个固定长度的特征向量。这套流程在 CWRU 的 10 类故障（正常 + 3 种故障类型 × 3 种直径）上做过验证，配合简单的 SVM 就能得到不错的分类结果。
from scipy.stats import kurtosis, skew
from scipy.signal import stft
import numpy as np
def extract_tf_features(signal, fs, nperseg=256):
"""
从时频分析结果中提取特征向量
返回: 特征字典
"""
features = {}
# 1. STFT 频带能量比
f, t, Zxx = stft(signal, fs=fs, nperseg=nperseg,
noverlap=int(nperseg*0.75))
mag = np.abs(Zxx)
total_energy = np.sum(mag**2) + 1e-10
# 定义四个频带：0-500, 500-1500, 1500-3000, 3000-6000 Hz
bands = [(0, 500), (500, 1500), (1500, 3000), (3000, 6000)]
for i, (lo, hi) in enumerate(bands):
mask = (f >= lo) & (f < hi)
band_energy = np.sum(mag[mask, :]**2)
features[f'band_energy_ratio_{i}'] = band_energy / total_energy
# 2. 时频熵
p = mag**2 / total_energy
features['tf_entropy'] = -np.sum(p * np.log2(p + 1e-10))
# 3. 目标频率切片的统计量
# 以 162 Hz 为例（内圈故障特征频率）
target_freq = 162
freq_mask = (f >= target_freq - 20) & (f <= target_freq + 20)
if np.any(freq_mask):
slice_mag = np.mean(mag[freq_mask, :], axis=0)
features['bpfi_slice_kurtosis'] = kurtosis(slice_mag)
features['bpfi_slice_skew'] = skew(slice_mag)
features['bpfi_slice_peak'] = np.max(slice_mag) / (np.mean(slice_mag) + 1e-10)
# 4. 时域统计量作为补充
features['time_kurtosis'] = kurtosis(signal)
features['time_rms'] = np.sqrt(np.mean(signal**2))
features['time_peak'] = np.max(np.abs(signal))
return features
# 使用示例
feat = extract_tf_features(signal, FS)
for k, v in feat.items():
print(f"{k}: {v:.4f}")
这段代码的逻辑分四层。第一层是频带能量比，把 STFT 幅值平方后按频带求和再除以总能量，得到四个比值，反映能量在不同频段的分布。轴承故障的冲击能量集中在高频段，正常信号能量集中在低频段，这个比值有区分度。第二层是时频熵，衡量时频平面上能量分布的均匀程度，故障冲击会让能量集中，熵值下降。第三层是目标频率切片的统计量，峰度反映冲击的尖锐程度，偏度反映不对称性，峰值因子反映冲击的相对强度。第四层是时域统计量，作为兜底特征。
参数方面， nperseg=256 是起点，如果你的信号采样率是 48 kHz，窗长要相应加大到 1024 点，保持时间窗宽度在 20 ms 左右。频带划分不是固定的，根据你的故障特征频率调整：如果 BPFI 在 162 Hz，第一个频带 0-500 Hz 能覆盖它；如果 BPFO 在 100 Hz 左右，也在第一个频带里。目标频率切片那里， target_freq 要按实际故障类型改，内圈用 BPFI，外圈用 BPFO，滚动体用 BSF。
这套特征向量长度在 12～15 之间，取决于频带数量和是否包含目标频率切片。我一般会把所有样本的特征向量拼成矩阵，做标准化后用 SVM 或随机森林分类。在 CWRU 的 0 负载数据上，10 类故障的分类准确率能到 97% 左右；跨负载（用 0 负载训练，1/2/3 负载测试）会降到 85%～90%，这是 CWRU 数据集的已知难点，不是特征提取的问题。
最后说一个我自己的习惯：每次做完时频分析，我会把原始信号、STFT 图、CWT 图、特征向量存成一个字典，用 joblib.dump 保存。这样后面换分类器或调参数时不用重新算时频变换，省时间。时频分析的计算量不小，尤其是 CWT，12 kHz 采样率下 10 秒信号做 256 个尺度的 CWT 要好几秒，缓存下来能省不少事。希望帮到你。
337
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
