信号处理基础之噪声与降噪(六) | DT-CWT、EWT降噪及python代码实现

7月前浏览8626

1 双树复小波降噪

1.1 双树复小波变换

双树复小波变换（Dual-Tree Complex Wavelet Transforms, DT-CWT）由实部树和虚部树的两个并行的实小波变换构成[1]。复小波表示为：

式中，h(t)和g(t)分别表示复小波的实部和虚部，两个函数均为实函数[2]。双树复小波变换包含两个平行的分解树，分别定义为树A和树B，树A给出双树复小波变换的实部，树B给出虚部，变换过程中，对树A和树B分别执行高低通滤波，即可实现双树分解，DT-CWT分解的流程如图1所示。

图1 双树复小波变换流程

DT-CWT显著改善了DWT的平移敏感性和方向选择，关于DWT的平移敏感性和方向选择的说明详见文献[3]，此处不再进行展开。

对变换中的每一棵树分别使用双正交滤波器（设计成与对应的分析滤波器具有完全重构特性的滤波器）进行反变换，最后对两棵树的输出结果进行平均，可以保证整个系统近似的平移不变性。经过双树复小波变换的实部树和虚部树的联合重构信号可以表示为：

式中，

和

为第j层分解的实部树和虚部树的联合，可分别表示为：

式中，

和

分别表示实数部小波变换的小波系数和尺度系数；

和

分别表示虚数部小波变换的小波系数和尺度系数。于是双树复小波变换分解的小波系数和尺度系数可表示为

可以看出，双树复小波变换的分解和重构实现了实部树和虚部树信息的互补，能够保持信号的完全重构性。

1.2 双树复小波降噪

双树复小波变换将原始信号分解成小波系数和尺度系数，也称为细节分量和近似分量，其中噪声成分主要集中在细节分量中（高频部分）。与小波降噪相同，随着分解层数的增加，噪声的能量逐渐衰减，双树复小波分解得到的小波系数变小，确定合理的阈值可以对噪声进行剔除。阈值确定方法详见信号处理基础之噪声与降噪(五) | 小波降噪及python代码实现。

双树复小波降噪的基本流程如下：

（1）信号分解：使用双树复小波变换对信号进行多级分解。DT-CWT 生成两棵树的小波系数，每棵树提供了信号的一半信息，这两棵树合在一起可以提供更完整的频率和方向信息；

（2） 阈值处理：对小波系数进行阈值处理以去除噪声。常用的阈值方法包括软阈值和硬阈值处理。软阈值可以更平滑地处理信号，而硬阈值可能会导致信号失真；

（3）系数重构：使用处理后的小波系数通过逆双树复小波变换重构信号。

1.3 双树复小波降噪的python实现

双树复小波降噪的python实现如下：

import pywtimport numpy as npimport matplotlib.pyplot as plt
def denoised_with_dtwt(signal, wavelet, mode='soft', level=1):    coeffs = pywt.wavedec(signal, wavelet, mode='per')    sigma = np.median(np.abs(coeffs[-1])) / 0.6745    threshold = sigma * np.sqrt(2 * np.log(len(signal)))    new_coeffs = [pywt.threshold(c, value=threshold, mode=mode) for c in coeffs]    new_coeffs[0] = coeffs[0]  # 保留近似系数不变    signal_denoised = pywt.waverec(new_coeffs, wavelet, mode='per')    return signal_denoised

2 经验小波降噪

2.1 经验小波变换

经验小波变换(Empirical Wavelet Transform, EWT)是由Gilles提出的一种基于小波框架的自适应自适应小波分析方法，其将小波变换理论与EMD方法相结合，规避了EMD计算量大和模态混叠等缺点[4]。EWT的实质是对信号频谱进行自适应分割，建立合适的小波滤波器组，加小波窗后提取紧支撑傅里叶频谱的调幅-调频成分[5]。

经验小波变换的基本步骤如下：

（1）傅里叶频谱的自适应划分：定义[0,

]的傅里叶谱将信号自适应分割成连续N段，每个频段表示为

，以每个

为中心，定义宽度为

的过渡段，如图2所示[6]。

图2 傅里叶轴的分割

（2）各分割区带通滤波器设计：依据Little wood-Paley和Meyer理论，构造经验小波函数

和经验尺度函数

，构造方法详见文献[4-6]。

（3）细节系数与近似系数确定：定义经验小波变换

，记经验小波函数

和经验尺度函数

的逆傅里叶变换分别为

和

，则细节系数可表示为经验小波函数与信号的内积：

近似系数可以表示为经验尺度函数与信号内积：

（4）原始信号重构：原始信号可表示为

经验模态可表示为：

2.2 经验小波降噪

在经验小波变换的基础上，结合峭度法和相关系数法等方法（详见信号处理基础之噪声与降噪(三) | EMD降噪与VMD降噪及python代码实现）可对信号中的噪声进行去除，基本步骤如下：

（1）通过EWT对信号进行分解，得到若干经验模态函数；

（2）计算各经验模态函数的峭度值或相关系数，筛选有用信息；

（3）信号重构。

2.3 经验小波降噪的python实现

经验小波降噪的python实现如下：

import numpy as npimport matplotlib.pyplot as pltimport ewtpy
def denoised_with_ewt(signal):    ewt, mfb, boundaries = ewtpy.EWT1D(signal)    ewt = ewt.T    ewt_denoised = [component for component in ewt if np.corrcoef(signal, component)[0, 1] >= 0.4]    signal_denoised = np.sum(ewt_denoised, axis=0)    return signal_denoised

3 算法测试

3.1 测试用例

注：用例中denoised_evaluate()函数见”信号处理基础之噪声与降噪(二) | 时域降噪方法(平滑降噪、SVD降噪)python代码实现“

import matplotlib.pyplot as pltimport matplotlib# 创建一个信号n = 500  # 生成500个点的信号t = np.linspace(0, 30*np.pi, n, endpoint=False)s = np.cos(0.1*np.pi*t) + np.sin(0.3*np.pi*t) + np.cos(0.5*np.pi*t) + np.sin(0.7*np.pi*t)  # 原始信号r = 0.5 * np.random.randn(n)y = s + r  # 加噪声
denoised_signal_dtwt = denoised_with_dtwt(y, 'db4', 'soft', level=8)denoised_signal_ewt = denoised_with_ewt(y)value_dtwt = denoise_evaluate(s, r, denoised_signal_dtwt)value_ewt = denoise_evaluate(s, r, denoised_signal_ewt)
print(value_dtwt)print(value_ewt)
fig = plt.figure(figsize=(16, 12))plt.subplots_adjust(hspace=0.5)plt.subplots_adjust(left=0.05, right=0.98, top=0.95, bottom=0.05)font = {'family': 'Times New Roman', 'size': 16, 'weight': 'normal',}matplotlib.rc('font', **font)
plt.subplot(5, 1, 1)plt.plot(t, s, linewidth=1.5, color='b')plt.title('pure signal', fontname='Times New Roman', fontsize=20)
plt.subplot(5, 1, 2)plt.plot(t, y, linewidth=1.5, color='b')plt.title('noised signal', fontname='Times New Roman', fontsize=20)
plt.subplot(5, 1, 3)plt.plot(t, denoised_signal_dtwt, linewidth=1.5, color='b')plt.title('DTWT', fontname='Times New Roman', fontsize=20)
plt.subplot(5, 1, 4)plt.plot(t, denoised_signal_ewt, linewidth=1.5, color='b')plt.title('EWT', fontname='Times New Roman', fontsize=20)
plt.subplot(5, 1, 5)plt.plot(t, y-denoised_signal_ewt, linewidth=1.5, color='b')plt.title('y-EWT', fontname='Times New Roman', fontsize=20)plt.show()

3.2 降噪结果

降噪结果如图3所示。

图3 降噪结果

3.3 降噪指标

降噪指标如表1所示。

表1 降噪指标

4 参考文献

[1] 刘嘉辉,秦仙蓉,王玉龙,等. 基于双树复小波变换与样本熵的自适应降噪法[J]. 振动、测试与诊断,2022,42(2):285-291.

[2] 刘蕾. 基于双树复小波变换的图像去噪[D]. 北京:北京化工大学,2010.

[3] 石宏理,胡波. 双树复小波变换及其应用综述[J]. 信息与电子工程,2007,5(3):229-234.

[4] 吕跃刚,何洋洋. EWT和ICA联合降噪在轴承故障诊断中的应用[J]. 振动与冲击,2019,38(16):42-48,70.

[5] 骆春林,刘其洪,李伟光. EWT-SSA联合降噪及其在滚刀振动信号分析中的应用[J]. 中国测试,2022,48(8):109-116.

[6] 李志农,朱明,褚福磊,等. 基于经验小波变换的机械故障诊断方法研究[J]. 仪器仪表学报,2014(11):2423-2432.

来源：故障诊断与python学习

HUST Bearing公开数据集(含不同转速、复合故障)

继前期推荐的两个公开数据集：哈工大航空发动机轴承数据集(HIT bearing dataset)北交风力发电机行星齿轮箱数据集(WT-planetary-gearbox-dataset)今天给大家推荐一个轴承故障公开数据集，该数据集有4种转速工况，9种健康状态，考虑了复合故障，诊断起来更具有挑战性。该数据集是华中科技大学沈卫明老师团队在2024年公开，因此目前基于该数据集的论文不是很多，小伙伴们赶紧用起来吧！对于研究轴承故障诊断的小伙伴们，再也不用担心写论文找不到数据啦！数据集基本信息该数据集包括轴承在4种不同工况条件下的9种不同健康状态的振动信号。这些数据集是公开的，任何人都可以使用它们来验证滚动轴承的诊断算法。使用HUSTbearing数据集的出版物请引用以下论文：参考论文：Domain Generalization for Cross-Domain Fault Diagnosis: an Application-oriented Perspective and a Benchmark Study论文期刊：Reliability Engineering and System SafetyDoi：https://doi.org/10.1016/j.ress.2024.109964作者：Chao Zhao(a c), Enrico Zio(b c), Weiming Shen(a)机构：a State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science & Technology, Wuhan 430074, Chinab MINES Paris PSL University, CRC, Sophia Antipolis, Francec Energy Department, Politecnico di Milano, Milan, Italy数据集下载：https://github.com/CHAOZHAO-1/HUSTbearing-dataset公开日期：2024年作者简介：沈卫明，(Fellow, IEEE) 1983年和1986年分别在北京交通大学获得学士和硕士学位，1996年在法国贡比涅工业大学获得博士学位。沈老师目前是华中科技大学的教授，以及加拿大西安大略大学的兼职教授。主要研究方向为智能软件代理、无线传感器网络、物联网、大数据及其在工业中的应用。沈教授是加拿大工程院和加拿大工程学院的院士，同时也是加拿大安大略省注册工程师。上述链接若打不开，文末有百度网盘链接。目录1 实验简介1.1 轴承故障实验台1.2 测试轴承参数1.3 传感器及相关设置1.4 工况条件1.5 采样设置2 数据集细节3 数据读取与展示3.1 数据读取3.2 绘制时域图3.3 绘制频域图3.4 绘制stft时频域图4 参考文献1 实验简介1.1 轴承故障实验台轴承故障试验使用Spectra-Quest机械故障实验台进行，如图1所示。图1 HUST轴承数据集实验台实验台上从左到右依次为①速度控制器、②电机、③轴、④加速度传感器、⑤轴承和⑥数据采集卡。9种健康状态的轴承如图2所示，分别表示(1) 正常，(2) 内圈中度故障，(3) 内圈重度故障，(4) 外圈中度故障，(5) 外圈重度故障，(6) 滚动体中度故障，(7) 滚动体重度故障，(8) 外圈内圈复合中度故障，(9) 外圈内圈复合重度故障。需要注意的是，复合故障表示内圈和外圈都存在故障。所有故障都是人为预设的。图2 故障轴承照片1.2 测试轴承参数被测轴承类型为ER-16K，详细参数见表1。表1 测试轴承参数[1]ER-16K深沟球轴承的轴承故障特征频率如表2所示，为转速，单位r/min。表2 轴承故障特征频率[2]1.3 传感器及相关设置用于数据采集的加速度传感器如图3所示，传感器模型如图4所示。具体信号采集设置如图5所示。图3 三向加速度传感器的照片图4 传感器模型图5 信号采集设置1.4 工况条件实验共设置4种不同的速度工况条件，工况条件包括:1) 65Hz (3900rpm)2) 70Hz (4200rpm)3) 75Hz (4500rpm)4) 80Hz (4800rpm)1.5 采样设置采样频率设置为25.6 kHz。如图6所示，每次采样共记录262144个数据点(即10.2s)。图6 文件描述2 数据集细节原始数据文件包括36个文件(9个健康状态乘以4个工作条件)，每个文件为Excel格式。例如，文件名“0.5X_B_65Hz”表示在65Hz工作条件下发生滚动体中度故障，其中0.5X表示中度故障。故障状态用以下代码表示：H：健康I：内圈故障O：外圈故障B：球故障C：组合故障。例如，“O_80Hz”表示在80Hz工作条件下，外圈重度故障。图7 部分文件示例3 数据读取与展示以0.5X_B_65Hz.xls数据为例，展示其时域图、频域图、stft时频域图。3.1 数据读取首先是先了解如何用python读取数据。它的格式是xls，用pd.read_csv()函数读取，但还需要再处理一下。下面定义了1个data_read()函数。## 导入包from matplotlib import pyplot as pltfrom matplotlib import rcParamsimport numpy as npimport pandas as pdimport osconfig = { "font.family": 'serif', # 衬线字体 "font.size": 14, # 相当于小四大小 "font.serif": ['SimSun'], # 宋体 "mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大 'axes.unicode_minus': False # 处理负号，即-号}rcParams.update(config)def data_read(file_path): """ :fun: 读取xls数据 :param file_path: 文件路径 :return df: Dataframe，五列分别为'Time', 'Speed', 'Acc_x', 'Acc_y', 'Acc_z' """ df = pd.read_csv(file_path)[17:] # 第17行开始为加速度数据 df.columns = ['Column'] # 使用 split 方法分割每行的单元格，expand=True 会将分割后的每个元素作为单独的列返回 df_split = df['Column'].str.split('\t', expand=True) # 将分割后的所有列转换为 float 类型 df_split = df_split.astype(float) # 重置索引，以便将分割后的列与原始 DataFrame 的索引对齐 df_split.reset_index(drop=True, inplace=True) # 现在 df_split 包含了分割后的五列数据 # 为这些新列设置列名 df_split.columns = ['Time', 'Speed', 'Acc_x', 'Acc_y', 'Acc_z'] return df_split file_path = r'E:\03-公开数据集\HUST-bearing-dataset\Raw data (原始数据)/0.5X_B_65Hz.xls'df = data_read(file_path)acc_z_arr= df['Acc_z'] # 选择z轴数据df图8 一个xls文件内具体内容可知，该数据为5列，分别为'Time', 'Speed', 'Acc_x', 'Acc_y', 'Acc_z'。3.2 绘制时域图选择z轴的加速度传感器数据进行展示。##========绘制时域信号图========##def plt_time_domain(arr, fs=1600, ylabel='Amp($m/s^2$)', title='原始数据时域图', img_save_path=None, x_vline=None, y_hline=None): """ :fun: 绘制时域图模板 :param arr: 输入一维数组数据 :param fs: 采样频率 :param ylabel: y轴标签 :param title: 图标题 :return: None """ import matplotlib.pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 plt.rcParams['axes.unicode_minus'] = False # 显示负号 font = {'family': 'Times New Roman', 'size': '20', 'color': '0.5', 'weight': 'bold'} plt.figure(figsize=(12,4)) length = len(arr) t = np.linspace(0, length/fs, length) plt.plot(t, arr, c='g') plt.xlabel('t(s)') plt.ylabel(ylabel) plt.title(title) if x_vline: plt.vlines(x=x_vline, ymin=np.min(arr), ymax=np.max(arr), linestyle='--', colors='r') if y_hline: plt.hlines(y=0.2, xmin=np.min(t), xmax=np.max(t), linestyle=':', colors='y') #===保存图片====# if img_save_path: plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight') plt.show()fs = 25600file_name = r'0.5X_B_65Hz.xls'##=====绘制时域数据====##time_img_save_path = file_path.replace('.xls', '时域_.png')plt_time_domain(acc_z_arr, fs=fs, title=file_name, img_save_path=time_img_save_path)共10.2s的数据，单位。3.3 绘制频域图##========绘制频域信号图========##def plt_fft_img(arr, fs, ylabel='Amp(mg)', title='频域图', img_save_path=None, vline=None, hline=None, xlim=None): """ :fun: 绘制频域图模板 :param arr: 输入一维时域数组数据 :param fs: 采样频率 :param ylabel: y轴标签 :param title: 图标题 :return: None """ # 计算频域幅值 length = len(arr) t = np.linspace(0, length/fs, length) fft_result = np.fft.fft(arr) fft_freq= np.fft.fftfreq(len(arr), d=t[1]-t[0]) # FFT频率 fft_amp= 2*np.abs(fft_result)/len(t) # FFT幅值 # 绘制频域图 plt.figure(figsize=(12,4)) plt.title(title) plt.plot(fft_freq[0: int(len(t)/2)], fft_amp[0: int(len(t)/2)], label='Frequency Spectrum', color='b') plt.xlabel('频率 (Hz)') plt.ylabel('幅值') plt.legend() if vline: plt.vlines(x=vline, ymin=np.min(fft_amp), ymax=np.max(fft_amp), linestyle='--', colors='r') if hline: plt.hlines(y=hline, xmin=np.min(fft_freq), xmax=np.max(fft_freq), linestyle=':', colors='y') #===保存图片====# if img_save_path: plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight') if xlim: # 图片横坐标是否设置xlim plt.xlim(0, xlim) plt.tight_layout() plt.show()fs = 25600file_name = r'0.5X_B_65Hz.xls'##=====绘制频域数据====##fft_img_save_path = file_path.replace('.xls', '频域_.png')plt_fft_img(acc_z_arr, fs=fs, title=file_name, img_save_path=fft_img_save_path)可见明显的主频及倍频。3.4 绘制stft时频域图def plt_stft_img(arr, fs, ylabel='Amp(mg)', title='频域图', img_save_path=None, vline=None, hline=None, xlim=None): """ :fun: 绘制stft时频域图模板 :param arr: 输入一维时域数组数据 :param fs: 采样频率 :param ylabel: y轴标签 :param title: 图标题 :return: None """ import scipy.signal as signal import numpy as np import matplotlib.pyplot as plt f, t, nd = signal.stft(arr, fs=fs, window='hann', nperseg=128, noverlap=64,nfft=None, detrend=False, return_onesided=True, boundary='odd', padded=False, axis=-1) # fs:时间序列的采样频率, nperseg:每个段的长度，默认为256(2^n) noverlap:段之间重叠的点数。如果没有则noverlap=nperseg/2 #window ：字符串或元组或数组，可选需要使用的窗。 # #如果window是一个字符串或元组，则传递给它window是数组类型，直接以其为窗，其长度必须是nperseg。 # 常用的窗函数有boxcar，triang，hamming， hann等，默认为Hann窗。 #nfft ：int，可选。如果需要零填充FFT，则为使用FFT的长度。如果为 None，则FFT长度为nperseg。默认为无 # detrend ：str或function或False，可选 # 指定如何去除每个段的趋势。如果类型参数传递给False，则不进行去除趋势。默认为False。 # return_onesided ：bool，可选 # 如果为True，则返回实际数据的单侧频谱。如果 False返回双侧频谱。默认为 True。请注意，对于复杂数据，始终返回双侧频谱。 # boundary ：str或None，可选 # 指定输入信号是否在两端扩展，以及如何生成新值，以使第一个窗口段在第一个输入点上居中。 # 这具有当所采用的窗函数从零开始时能够重建第一输入点的益处。 # 有效选项是['even', 'odd', 'constant', 'zeros', None]. # 默认为‘zeros’,对于补零操作[1, 2, 3, 4]变成[0, 1, 2, 3, 4, 0] 当nperseg=3. # padded：bool，可选 # 指定输入信号在末尾是否填充零以使信号精确地拟合为整数个窗口段，以便所有信号都包含在输出中。默认为True。 # 填充发生在边界扩展之后，如果边界不是None，则填充为True，默认情况下也是如此。 # axis ：int，可选 # 绘制STFT时频域图 plt.figure(figsize=(12,4)) plt.pcolormesh(t, f, np.abs(nd), vmin = np.min(np.abs(nd)), vmax = np.max(np.abs(nd))) plt.title(title) plt.xlabel('时间（t）') plt.ylabel('频率 (Hz)') if vline: plt.vlines(x=vline, ymin=np.min(fft_amp), ymax=np.max(fft_amp), linestyle='--', colors='r') if hline: plt.hlines(y=hline, xmin=np.min(fft_freq), xmax=np.max(fft_freq), linestyle=':', colors='y') #===保存图片====# if img_save_path: plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight') if xlim: # 图片横坐标是否设置xlim plt.xlim(0, xlim) plt.tight_layout() plt.show()fs = 25600file_name = r'0.5X_B_65Hz.xls'##=====绘制STFT时频域数据====##stft_img_save_path = file_path.replace('.xls', '时频域_.png') plt_stft_img(acc_z_arr, fs=fs, img_save_path=stft_img_save_path)与fft频域图一致，频率主要集中在低频(0-1000Hz)。4 参考文献[1] Luo W, Yan C, Yang J, et al. Vibration response of defect-ball-defect of rolling bearing with compound defects on both inner and outer races[C]//IOP Conference Series: Materials Science and Engineering. IOP Publishing, 2021, 1207(1): 012006.[2] Mishra C, Samantaray A K, Chakraborty G. Ball bearing defect models: A study of simulated and experimental fault signatures[J]. Journal of Sound and Vibration, 2017, 400: 86-112.来源：故障诊断与python学习