目录
前言
1 短时傅里叶变换STFT原理介绍
1.1 傅里叶变换的本质
1.2 STFT概述
1.3 STFT的原理和过程
1.3.1 时间分割
1.3.2 傅里叶变换
1.3.3 时频图：
1.4 公式表示
2 基于Python的STFT实现与参数对比
2.1 代码示例
2.2 参数选择和对比
2.2.1 nperseg（窗口长度）：
2.2.2 noverlap（重叠长度）：
2.2.3 选择策略：
2.3 凯斯西储大学轴承数据的加载
2.4 STFT与参数选择
2.4.1 基于重叠比例为0.5，选择内圈数据比较 STFT 的不同尺度：16、32 、64、128
2.4.1 根据正常数据和三种故障数据，对比不同尺度的辨识度
3 基于时频图像的轴承故障诊断分类
3.1 生成时频图像数据集
3.2 定义数据加载器和VGG网络模型
3.3 设置参数，训练模型
往期精彩内容：
Python-凯斯西储大学（CWRU）轴承数据解读与分类处理
Python轴承故障诊断 (一)短时傅里叶变换STFT
Python轴承故障诊断 (二)连续小波变换CWT_pyts 小波变换 故障-CSDN博客
Python轴承故障诊断 (三)经验模态分解EMD_轴承诊断 pytorch-CSDN博客
Pytorch-LSTM轴承故障一维信号分类(一)_cwru数据集pytorch训练-CSDN博客
Pytorch-CNN轴承故障一维信号分类(二)-CSDN博客
Pytorch-Transformer轴承故障一维信号分类(三)-CSDN博客
Python轴承故障诊断 (四)基于EMD-CNN的故障分类-CSDN博客
Python轴承故障诊断 (五)基于EMD-LSTM的故障分类-CSDN博客
Python轴承故障诊断 (六)基于EMD-Transformer的故障分类-CSDN博客
Python轴承故障诊断 (七)基于EMD-CNN-LSTM的故障分类-CSDN博客
Python轴承故障诊断 (八)基于EMD-CNN-GRU并行模型的故障分类-CSDN博客
基于FFT + CNN - BiGRU-Attention 时域、频域特征注意力融合的轴承故障识别模型-CSDN博客
基于FFT + CNN - Transformer 时域、频域特征融合的轴承故障识别模型-CSDN博客
大甩卖-（CWRU）轴承故障诊数据集和代码全家桶-CSDN博客
Python轴承故障诊断 (九)基于VMD+CNN-BiLSTM的故障分类-CSDN博客
Python轴承故障诊断 (十)基于VMD+CNN-Transfromer的故障分类-CSDN博客
Python轴承故障诊断 (11)基于VMD+CNN-BiGRU-Attenion的故障分类-CSDN博客
交叉注意力融合时域、频域特征的FFT + CNN -BiLSTM-CrossAttention轴承故障识别模型-CSDN博客
交叉注意力融合时域、频域特征的FFT + CNN-Transformer-CrossAttention轴承故障识别模型-CSDN博客
轴承故障诊断 (12)基于交叉注意力特征融合的VMD+CNN-BiLSTM-CrossAttention故障识别模型-CSDN博客
Python轴承故障诊断入门教学-CSDN博客
Python轴承故障诊断 (13)基于故障信号特征提取的超强机器学习识别模型-CSDN博客
Python轴承故障诊断 (14)高创新故障识别模型-CSDN博客
Python轴承故障诊断 (15)基于CNN-Transformer的一维故障信号识别模型-CSDN博客
前言
本文基于凯斯西储大学（CWRU）轴承数据，进行短时傅里叶变换的介绍与参数选择，最后通过Python实现对故障数据的时频图像分类。凯斯西储大学（CWRU）轴承数据的详细介绍可以参考下文：
Python-凯斯西储大学（CWRU）轴承数据解读与分类处理_cwru轴承数据集-CSDN博客
1 短时傅里叶变换STFT原理介绍
1.1 傅里叶变换的本质
傅里叶变换的基本思想：
将信号分解成一系列不同频率的连续正弦波的叠加；
或者说，将信号从时间域转换到频率域。
由于傅里叶变换是对整个信号进行变换，将整个信号从时域转换到频域，得到一个整体的频谱；丢掉了时间信息，无法根据傅立叶变换的结果判断一个特定信号在什么时候发生；所以傅里叶变换缺乏时频分析能力、多分辨率分析能力，难以分析非平稳信号。
1.2 STFT概述
短时傅里叶变换（Short-Time Fourier Transform，STFT）是一种将信号分解为时域和频域信息的时频分析方法。它通过将信号分成短时段，并在每个短时段上应用傅里叶变换来捕捉信号的瞬时频率。即采用中心位于时间α的时间窗g(t-α)在时域信号上滑动，在时间窗g(t-α)限定的范围内进行傅里叶变换，这样就使短时傅里叶变换具有了时间和频率的局部化能力，兼顾了时间和频率的分析[1]。
使用窄窗，时间分辨率提高而频率分辨率降低；
使用宽窗，频率分辨率提高而时间分辨率降低。
比如用利用高斯窗STFT对非平稳信号进行分析：
1.3 STFT的原理和过程
1.3.1 时间分割
在短时傅里叶变换（STFT）中，首先将信号分割成短时段。这个过程通常通过使用窗口函数来实现，窗口函数是一个在有限时间内非零，而在其他时间上为零的函数。常见的窗口函数有矩形窗、汉明窗、汉宁窗等。通过将窗口函数应用于信号，可以将信号分成许多短时段。
1.3.2 傅里叶变换
对于每个短时段，都会进行傅里叶变换。傅里叶变换是一种将信号从时域（时间域）转换为频域（频率域）的方法。在这个上下文中，它用于分析每个短时段内信号的频率成分。傅里叶变换将信号表示为不同频率的正弦和余弦函数的组合。
1.3.3 时频图：
将每个短时段的傅里叶变换结果排列成一个矩阵，构成了时频图。时频图的横轴表示时间，纵轴表示频率，而每个点的强度表示对应频率在对应时刻的幅度。时频图提供了一种直观的方式来观察信号在时间和频率上的变化。
1.4 公式表示
STFT的数学表达式如下：
其中，
2 基于Python的STFT实现与参数对比
在 Python 中，使用 scipy 库来实现短时傅里叶变换（STFT）
2.1 代码示例
import numpy as npimport matplotlib.pyplot as pltfrom scipy.signal import stft# 生成三个不同频率成分的信号fs = 1000 # 采样率t = np.linspace(0, 1, fs, endpoint=False) # 时间# 第一个频率成分signal1 = np.sin(2 * np.pi * 50 * t)# 第二个频率成分signal2 = np.sin(2 * np.pi * 120 * t)# 第三个频率成分signal3 = np.sin(2 * np.pi * 300 * t)# 合并三个信号signal = np.concatenate((signal1, signal2, signal3))# 进行短时傅里叶变换 f, t, spectrum = stft(signal, fs, nperseg=100, noverlap=50)# 绘制时频图plt.pcolormesh(t, f, np.abs(spectrum), shading='gouraud')plt.title('STFT Magnitude')plt.ylabel('Frequency [Hz]')plt.xlabel('Time [sec]')plt.show()参数解释
nperseg：表示每个段的长度，即窗口大小
noverlap：表示相邻段的重叠部分
f：是频率数组，表示频谱的频率分量
t：是时间数组，表示每个时间段的起始时间
spectrum：是频谱矩阵，spectrum[f, t] 表示在频率 f 处于时间 t 时的频谱强度
可以使用 np.abs(spectrum) 获取频谱的幅度，表示在不同频率和时间上的信号强度
2.2 参数选择和对比
2.2.1 nperseg（窗口长度）：
nperseg 定义了每个时间段内的信号长度，也就是说，它规定了窗口的大小。
较小的 nperseg 可以提供更高的时间分辨率，因为窗口更短，可以捕捉到信号的更快变化。但频率分辨率较差。
较大的 nperseg 提供更高的频率分辨率，但可能失去对信号快速变化的捕捉。
选择适当的窗口长度要根据信号特性，如频率变化的速率和信号长度等。
2.2.2 noverlap（重叠长度）：
noverlap 定义了相邻时间段之间的重叠部分。
较大的 noverlap 可以提高时间分辨率，因为相邻时间段之间有更多的共享信息。但可能导致频谱图中的泄漏（leakage）。
较小的 noverlap 可以提高频率分辨率，减少泄漏，但时间分辨率可能下降。
2.2.3 选择策略：
参数的选择需要考虑到对信号分析的具体需求，平衡时间和频率分辨率，尝试不同的 nperseg 和 noverlap 值，观察结果的变化。
对于 nperseg，选择较小的值可能需要根据信号的特定频率内容进行调整。确保窗口足够短，以捕捉到频率变化。
对于 noverlap，通常选择 50% 的重叠是一个常见的起点
2.3 凯斯西储大学轴承数据的加载
选择正常信号和 0.021英寸内圈、滚珠、外圈故障信号数据来做对比
第一步，导入包，读取数据
import numpy as npfrom scipy.io import loadmatimport numpy as npfrom scipy.signal import stftimport matplotlib.pyplot as pltimport matplotlibmatplotlib.rc("font", family='Microsoft YaHei')# 读取MAT文件 data1 = loadmat('0_0.mat') # 正常信号data2 = loadmat('21_1.mat') # 0.021英寸 内圈data3 = loadmat('21_2.mat') # 0.021英寸 滚珠data4 = loadmat('21_3.mat') # 0.021英寸 外圈# 注意，读取出来的data是字典格式，可以通过函数type(data)查看。第二步，数据集中统一读取 驱动端加速度数据，取一个长度为1024的信号进行后续观察和实验
# DE - drive end accelerometer data 驱动端加速度数据data_list1 = data1['X097_DE_time'].reshape(-1)data_list2 = data2['X209_DE_time'].reshape(-1) data_list3 = data3['X222_DE_time'].reshape(-1)data_list4 = data4['X234_DE_time'].reshape(-1)# 划窗取值（大多数窗口大小为1024）data_list1 = data_list1[0:1024]data_list2 = data_list2[0:1024]data_list3 = data_list3[0:1024]data_list4 = data_list4[0:1024]第三步，进行数据可视化
plt.figure(figsize=(20,10))plt.subplot(2,2,1)plt.plot(data_list1)plt.title('正常')plt.subplot(2,2,2)plt.plot(data_list2)plt.title('内圈')plt.subplot(2,2,3)plt.plot(data_list3)plt.title('滚珠')plt.subplot(2,2,4)plt.plot(data_list4)plt.title('外圈')plt.show()
2.4 STFT与参数选择
2.4.1 基于重叠比例为0.5，选择内圈数据比较 STFT 的不同尺度：16、32 、64、128
from scipy.signal import stft# 设置STFT参数window_size = 32 # 窗口大小overlap = 0.5 # 重叠比例# 计算重叠的样本数overlap_samples = int(window_size * overlap)frequencies1, times1, magnitude1 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples)# 设置STFT参数window_size = 64 # 窗口大小overlap = 0.5 # 重叠比例# 计算重叠的样本数overlap_samples = int(window_size * overlap)frequencies2, times2, magnitude2 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples)# 设置STFT参数window_size = 128 # 窗口大小overlap = 0.5 # 重叠比例# 计算重叠的样本数overlap_samples = int(window_size * overlap)frequencies3, times3, magnitude3 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples)# 设置STFT参数window_size = 256 # 窗口大小overlap = 0.5 # 重叠比例# 计算重叠的样本数overlap_samples = int(window_size * overlap)frequencies4, times4, magnitude4 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples数据可视化
plt.figure(figsize=(20,10), dpi=100)plt.subplot(2,2,1)plt.pcolormesh(times1, frequencies1, np.abs(magnitude1), shading='gouraud')plt.title('尺度 16-内圈')plt.subplot(2,2,2)plt.pcolormesh(times2, frequencies2, np.abs(magnitude2), shading='gouraud')plt.title('尺度 32-内圈')plt.subplot(2,2,3)plt.pcolormesh(times3, frequencies3, np.abs(magnitude3), shading='gouraud')plt.title('尺度 64-内圈')plt.subplot(2,2,4)plt.pcolormesh(times4, frequencies4, np.abs(magnitude4), shading='gouraud')plt.title('尺度 128-内圈')plt.show()
对比不同尺度的变化，大尺度有着较高的频率分辨率，小的尺度有着较高的时间分辨率，要权衡故障的分类辨识度，需要进一步做对比。
2.4.1 根据正常数据和三种故障数据，对比不同尺度的辨识度
综合考虑频率分辨率和时间分辨率，选择尺度为32。
3 基于时频图像的轴承故障诊断分类
下面以短时傅里叶变换（Short Time Fourier Transform，STFT）作为轴承故障数据的处理方法进行讲解：
数据介绍，凯斯西储大学（CWRU）轴承数据10分类数据集
train_set、val_set、test_set 均为按照7：2：1划分训练集、验证集、测试集，最后保存数据
3.1 生成时频图像数据集
如图所示为生成的时频图像数据集
3.2 定义数据加载器和VGG网络模型
制作数据标签，保存数据
定义VGG网络模型
3.3 设置参数，训练模型
100个epoch，准确率将近92%，继续调参可以进一步提高分类准确率
代码、数据如下：
对数据集和代码感兴趣的，可以关注最后一行
# 加载数据import torchfrom joblib import dump, loadimport torch.utils.data as Dataimport numpy as npimport pandas as pdimport torchimport torch.nn as nn# 参数与配置torch.manual_seed(100) # 设置随机种子，以使实验结果具有可重复性device = torch.device("cuda" if torch.cuda.is_available() else "cpu") #代码和数据集：https://mbd.pub/o/bread/ZZeXm5tv
参考文献
[1]《小波分析及其工程应用》.机械工业出版社
