Python轴承故障诊断 (一)短时傅里叶变换STFT

03-27 1766阅读

目录

前言

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 np
import matplotlib.pyplot as plt
from 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 np from scipy.io import loadmat import numpy as np from scipy.signal import stft import matplotlib.pyplot as plt import matplotlib matplotlib.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%，继续调参可以进一步提高分类准确率
        
        参考文献
        
        [1]《小波分析及其工程应用》.机械工业出版社

VPS购买请点击我

文章版权声明：除非注明，否则均为主机测评原创文章，转载或复制请以超链接形式并注明出处。

相关阅读

目录[+]