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

13天前浏览932

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 pltdef 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 ewtpydef 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()