
import numpy as np
 
import matplotlib.pyplot as plt
 
import pywt
 
 
# 生成三个不同频率成分的信号  3000个点
 
fs = 1000  # 采样率
 
time = np.linspace(0, 1, fs, endpoint=False)  # 时间
 
# 第一个频率成分
 
signal1 = np.sin(2 * np.pi * 30 * time)
 
# 第二个频率成分
 
signal2 = np.sin(2 * np.pi * 60 * time)
 
# 第三个频率成分
 
signal3 = np.sin(2 * np.pi * 120 * time)
 
# 合并三个信号
 
signal = np.concatenate((signal1, signal2, signal3))
 
 
# 连续小波变换参数
 
# 采样频率
 
sampling_rate = 3000
 
# 尺度长度
 
totalscal = 128    
 
# 小波基函数
 
wavename = 'morl'
 
# 小波函数中心频率
 
fc = pywt.central_frequency(wavename)
 
# 常数c
 
cparam = 2 * fc * totalscal  
 
# 尺度序列
 
scales = cparam / np.arange(totalscal, 0, -1)
 
 # 进行CWT连续小波变换
 
coefficients, frequencies = pywt.cwt(signal, scales, wavename, 1.0/1000)
 
# 小波系数矩阵绝对值
 
amp = abs(coefficients)
 
 
# 根据采样频率 sampling_period 生成时间轴 t
 
t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)
 
# 绘制时频图谱
 
plt.figure(figsize=(20,10))
 
plt.subplot(2,1,1)
 
plt.plot(signal)
 
plt.title('30Hz和60Hz和120Hz的分段波形')
 
plt.subplot(2,1,2)
 
plt.contourf(t, frequencies, amp, cmap='jet')
 
plt.title('对应时频图')
 
plt.show()

参数解释

totalscal：表示尺度长度，表征频率的参数，选择不同的尺度可以捕捉信号不同尺度上的特征；
wavename：表示小波函数的类型，不同的小波适用于不同类型的信号；
fc：表示小波函数的中心频率，小波函数在频率域中有一个中心频率；这个中心频率是与小波函数形状和性质有关的一个重要参数；
scales：表示尺度序列，CWT本质上是将你的信号与不同尺度的小波进行相关，scales 参数确定尺度范围；

coefficients：表示信号变换后的小波系数；
frequencies ：表示对应的频率信息

2.2 参数介绍和选择策略

2.2.1 尺度长度：

在连续小波变换（CWT）中，尺度参数是一个关键的选择，因为它决定了小波函数的宽度，从而影响了频率分辨率。尺度与频率成反比，尺度反映了分析的频率范围，尺度越小，小波函数衰减越快，频率越高；尺度越大，小波函数衰减越慢，频率越低[1]。

选择小波尺度的一般原则是：

高频特征：如果关注信号的高频特征，应该选择较小的尺度；
低频特征：如果关注信号的低频特征，应该选择较大的尺度；

覆盖感兴趣的频率范围：尺度参数的选择应该使小波函数能够覆盖感兴趣的频率范围；如果期望信号有很高的频率变化，可能需要选择较小的尺度；
频率分辨率： 较小的尺度提供更好的频率分辨率，因为小波函数较窄，可以更精细地定位频率；但是，这也意味着在时间上的分辨率较差，因此，需要权衡时间分辨率和频率分辨率；
信号持续时间：尺度参数的选择还应考虑信号的持续时间，如果信号是短暂的，可能需要较小的尺度；
尺度间隔： 在尺度参数上选择合适的间隔，以确保在整个频率范围内进行了适当的采样；这取决于具体的应用和信号特性。

2.2.2 小波函数（wavelet）：

小波函数（wavelet）的选择也连续小波变换中的一个重要参数，它决定了小波基函数的形状，不同的小波函数适用于不同类型的信号和应用。

打印 Python pywt 包中的所有小波函数类型：


import pywt
 
# 获取小波函数列表
 
wavelets = pywt.wavelist()
 
# 打印小波函数列表
 
# 按照12行，10列的形式打印数据
 
num_rows = 12
 
num_columns = 10
 
for i in range(0, len(wavelets), num_columns):
 
    row = wavelets[i:i + num_columns]
 
    print(row)

介绍常用的小波基函数：

'morl' ：Morlet小波是一种复杂的小波函数，它在频率域和时域都有较好的局部化性质；Morlet小波通常用于处理时频局部化要求较高的信号，比如处理振动信号或某些生物医学信号；
'cmor': Complex Morlet wavelet，复数 Morlet 小波，是 Morlet 小波的一个变种；
'cgau' ：Complex Gaussian wavelet，复数高斯小波，用于近似高斯信号；
'db1'：Daubechies小波是离散小波变换（Discrete Wavelet Transform, DWT）中常用的小波函数。Daubechies小波是紧支撑的小波，适用于处理有限长度的信号。
'haar'：Haar小波是最简单的小波函数之一，适用于对信号进行基本的低通和高通分解；
'mexh'：Mexican Hat小波，也称为Ricker小波，适用于处理具有尖峰或波包特性的信号；

'bior1.1'：Bior小波是一类双正交小波，其特点是具有对称和非对称两组滤波器，尤其适用于一些信号的多分辨率分析，如图像处理；
'sym2'：Symlet小波（Symmetric Wavelets），是一类对称的小波函数，它在某些方面类似于Daubechies小波，但是Symlet小波在设计上更加灵活。Symlet小波也是一种紧支撑小波，适用于有限长度的信号处理。

小波函数的选择通常取决于处理的信号类型以及分析的目标。在实际应用中，可以尝试不同的小波函数，观察它们在信号上的效果，然后根据实验结果选择最适合的小波函数。在选择小波函数时，也要考虑小波函数的性质，如平滑性、局部化等。

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 CWT与参数选择对比

本实验以某轴承故障诊断论文推荐'cgau8'小波，以及'morl'小波、'cmor1-1'小波、'cmor1.5-2'小波为实验做对比，尺度先设定为128，来对比不同小波函数的影响。

2.4.1 基于尺度为128，选择内圈数据比较 CWT 的不同小波函数


import numpy as np
 
import matplotlib.pyplot as plt
 
import pywt
 
import pandas as pd
 
 
# 连续小波变换参数
 
# 采样频率
 
sampling_rate = 1024
 
# 尺度长度
 
totalscal = 128  
 
 
wavename1 = 'cgau8'
 
fc1 = pywt.central_frequency(wavename1)
 
cparam1 = 2 * fc1 * totalscal  
 
scales1 = cparam1 / np.arange(totalscal, 0, -1)
 
 
wavename2 = 'morl' #
 
fc2 = pywt.central_frequency(wavename2)
 
cparam2 = 2 * fc2 * totalscal  
 
scales2 = cparam2 / np.arange(totalscal, 0, -1)
 
 
wavename3 = "cmor1-1"
 
fc3 = pywt.central_frequency(wavename3)
 
cparam3 = 2 * fc3 * totalscal  
 
scales3 = cparam3 / np.arange(totalscal, 0, -1)
 
 
wavename4 = 'cmor1.5-2'
 
fc4 = pywt.central_frequency(wavename4)
 
cparam4 = 2 * fc4 * totalscal  
 
scales4 = cparam4 / np.arange(totalscal, 0, -1)
 
 
# 进行连续小波变换
 
coefficients1, frequencies1 = pywt.cwt(data_list2, scales1, wavename1, sampling_period)
 
coefficients2, frequencies2 = pywt.cwt(data_list2, scales2, wavename2, sampling_period)
 
coefficients3, frequencies3 = pywt.cwt(data_list2, scales3, wavename3, sampling_period)
 
coefficients4, frequencies4 = pywt.cwt(data_list2, scales4, wavename4, sampling_period)
 
 
# 小波系数矩阵绝对值
 
amp1 = abs(coefficients1)
 
amp2 = abs(coefficients2)
 
amp3 = abs(coefficients3)
 
amp4 = abs(coefficients4)
 
 
# 根据采样频率 sampling_period 生成时间轴 t
 
t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)
 
数据可视化
 
plt.figure(figsize=(20,10))
 
plt.subplot(2,2,1)
 
plt.contourf(t, frequencies1, amp1, cmap='jet')
 
plt.title('内圈-cgau8')
 
plt.subplot(2,2,2)
 
plt.contourf(t, frequencies2, amp2, cmap='jet')
 
plt.title('内圈-morl')
 
plt.subplot(2,2,3)
 
plt.contourf(t, frequencies3, amp3, cmap='jet')
 
plt.title('内圈-cmor1-1')
 
plt.subplot(2,2,4)
 
plt.contourf(t, frequencies4, amp4, cmap='jet')
 
plt.title('内圈-cmor1.5-2')
 
plt.show()

对比不同小波函数，对于内圈故障信号来说，'cgau8'小波有着较高的频率分辨率，需要对比其他类型故障数据，进一步观察。

2.4.2 根据正常数据和三种故障数据，对比不同小波函数的辨识度

比较来看，'cgau8'小波和'cmor1.5-2'小波对于不同故障都有着较高的辨识度，综合考虑频率分辨率和时间分辨率，选择'cmor1.5-2'小波来进一步分析。

2.4.3 基于'cmor1.5-2'小波，选择滚珠故障数据比较 CWT 的不同尺度的变化：32、64、128、256；

同时尺度序列scales的常数项cparams均采用：cparam = 2 * fc * totalscal；'cmor1.5-2'小波中1 代表中心频率参数，1.5代表带宽参数；

注意，在连续小波变换中，影响小波系数的是尺度序列scales，但是仍能从上图中看出尺度长度越大，对低频特征关注越多（注意每幅小图的底部低频的区别）；

为进一步探索尺度序列scales的影响，选择尺度长度为128，设置不同常数项cparams来观察对时频变换的影响：


import numpy as np
 
import matplotlib.pyplot as plt
 
import pywt
 
import pandas as pd
 
 
# 连续小波变换参数
 
# 采样频率
 
sampling_rate = 1024
 
# 尺度长度
 
totalscal = 128  
 
 
wavename1 = 'cmor1.5-2'
 
fc1 = pywt.central_frequency(wavename1)
 
cparam1 = 1 * fc1 * totalscal
 
scales1 = cparam1 / np.arange(totalscal, 0, -1)
 
 
wavename2 = 'cmor1.5-2' #
 
fc2 = pywt.central_frequency(wavename2)
 
cparam2 = 2 * fc2 * totalscal
 
scales2 = cparam2 / np.arange(totalscal, 0, -1)
 
 
wavename3 = "cmor1.5-2"
 
fc3 = pywt.central_frequency(wavename3)
 
cparam3 = 3 * fc3 * totalscal  
 
scales3 = cparam3 / np.arange(totalscal, 0, -1)
 
 
wavename4 = 'cmor1.5-2'
 
fc4 = pywt.central_frequency(wavename4)
 
cparam4 = 4 * fc4 * totalscal  
 
scales4 = cparam4 / np.arange(totalscal, 0, -1)
 
 
# 进行连续小波变换
 
coefficients1, frequencies1 = pywt.cwt(data_list3, scales1, wavename1, sampling_period)
 
coefficients2, frequencies2 = pywt.cwt(data_list3, scales2, wavename2, sampling_period)
 
coefficients3, frequencies3 = pywt.cwt(data_list3, scales3, wavename3, sampling_period)
 
coefficients4, frequencies4 = pywt.cwt(data_list3, scales4, wavename4, sampling_period)
 
 
# 小波系数矩阵绝对值
 
amp1 = abs(coefficients1)
 
amp2 = abs(coefficients2)
 
amp3 = abs(coefficients3)
 
amp4 = abs(coefficients4)
 
 
# 根据采样频率 sampling_period 生成时间轴 t
 
t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)
 
进行可视化
 
plt.figure(figsize=(20,10), dpi=300)
 
plt.subplot(2,2,1)
 
plt.contourf(t, frequencies1, amp1, cmap='jet')
 
plt.title('滚珠-32')
 
plt.subplot(2,2,2)
 
plt.contourf(t, frequencies2, amp2, cmap='jet')
 
plt.title('滚珠-64')
 
plt.subplot(2,2,3)
 
plt.contourf(t, frequencies3, amp3, cmap='jet')
 
plt.title('滚珠-128')
 
plt.subplot(2,2,4)
 
plt.contourf(t, frequencies4, amp4, cmap='jet')
 
plt.title('滚珠-256')
 
plt.show()

对于不同尺度序列scales的比较，更能说明之前的结论：

高频特征：如果关注信号的高频特征，应该选择较小的尺度；
低频特征：如果关注信号的低频特征，应该选择较大的尺度；

对于滚珠故障类型数据，从时频图结果来看，应该选择2倍的cparams参数，有着较高的频率分辨率，和我们感兴趣的频率区域。

2.4.4 比较cmor小波函数不同参数 ----中心频率，带宽参数

在轴承故障诊断中，中心频率和带宽的具体设置取决于多个因素，包括轴承类型、工作条件和故障特征等。由于每个应用场景和故障类型都有所不同，没有一个通用的固定数值。以下是一些常见的参考范围和建议：

中心频率（Center Frequency）：

对于滚动轴承，常见的故障频率范围在几百赫兹到几千赫兹之间。可以选择在这个范围内的适当值作为中心频率。
根据具体的故障类型，例如滚动体故障、内圈故障或外圈故障，可能需要设置不同的中心频率。

带宽（Bandwidth）：

带宽的选择取决于信号的特征和所需的时频分辨率。
通常情况下，较小的带宽值（如1-5）适用于较短时域特征，能够提供更好的时频局部化能力。
如果需要更好的频率分辨率，可以选择较大的带宽值（如10-20），但可能会牺牲一部分时频局部化能力。

这些值仅供参考，实际应用时需要进行实验和调整，根据信号的频谱特征和故障频率范围来确定最佳的参数配置。建议通过观察生成的时频图，确保故障频率和特征得到适当的捕捉和展示。同时，根据实际的故障案例和经验，不断调整参数以提高故障诊断的准确性和可靠性。

经过大量的对比实验和观察，本文得出最后的参数结论设置：


# 尺度长度
 
totalscal = 128    
 
# 小波基函数
 
wavename = 'cmor100-1'
 
# 小波函数中心频率
 
fc = pywt.central_frequency(wavename)
 
# 常数c
 
cparam = 2 * fc * totalscal  
 
# 小波尺度序列
 
scales = cparam / np.arange(totalscal, 0, -1)