首页/文章/ 详情

信号处理基础之噪声与降噪(五) | 小波降噪及python代码实现

9月前浏览6972

本文摘要(由AI生成):

本文介绍了小波降噪方法,包括小波分解和小波降噪的基本原理,阈值的确定方法,以及软阈值和硬阈值两种处理方式。文章还提供了Python代码实现。


    接上期信号处理基础之噪声与降噪(四) | 进击的EMD族降噪及python代码实现本期为大家介绍小波降噪,并给出python代码。后续将会介绍基于DTWT和EWT的降噪应用,敬请关注。由于作者水平有限,文章的原理介绍与python代码难免有误,欢迎交流讨论!

目录

1 小波分解与小波降噪

1.1 小波分解的基本原理

1.2 小波降噪的基本原理

1.3 小波降噪的基本过程

2 阈值的确定

2.1 固定阈值原则与python实现

2.2 无偏似然估计原则与python实现

2.3 极值阈值原则与python实现

2.4 启发式阈值原则

3 阈值函数介绍

3.1 软阈值方法

3.2 硬阈值方法

4 算法测试

4.1 测试用例

4.2 降噪结果

4.3 降噪指标

5 参考文献

1 小波分解与小波降噪

1.1 小波分解的基本原理

小波降噪方法通过小波分解将噪声信号从原始信号中分离出来,信号经小波分解后可以得到一系列高频信号和低频信号,小波分解的过程如图1所示,记原始信号为S,第i次分解后的低频信号为    ,第i次分解后的高频信号为    ,信号采样频率为fs。

图1 小波分解基本原理

1.2 小波降噪的基本原理

      噪声和信号在小波域中形态的表现不同,它们的小波系数的幅值随尺度变化的趋势不同:随着尺度增加,真实信号的小波系数的幅值基本不变,而噪声的小波系数的幅值将很快衰减为零[1]。利用该机理,构造相应的降噪准则,对含噪的小波系数进行处理(减小甚至完全剔除),最终达到降噪的目的。

      注:小波变换具有线性特性,含噪信号经过小波分解后的小波系数等价于噪声的小波变换和信号的小波变换的代数和[1]

      小波降噪的关键是确定降噪准则,目前的方法较多,本文仅对阈值法降噪进行简单介绍,常用的阈值选取原则有4种[2]:固定阈值原则、无偏似然估计原则、极值阈值原则和启发式阈值原则。在阈值确定后,需采用阈值量化函数对小波系数进行处理,常用的处理方式有硬阈值函数和软阈值函数两种。

1.3 小波降噪的基本过程

      (1)确定小波分解的层数,对信号进行分解;
      (2)确定各个分解层下细节信号的阈值,对细节信号进行阈值量化处理;
      (3)利用阈值处理后的细节信号和逼近信号重构,得到降噪后的信号。
注:
      (1)第(1)步还可能涉及最优小波基的选取,这里不进行展开,具体可参考文献[1]的2.3.3章节;
      (2)文中仅给出了两种最常见的阈值函数,更多的方法此处不进行展开,具体可参考文献[3]的2.3章节。

2 阈值的确定

2.1 固定阈值原则与python实现
      固定阈值原则
    
式中,σ为噪声标准偏差,N为信号长度。
      固定阈值小波降噪的python实现如下:













import pywtimport numpy as np
def denoised_with_pywt_universal(data, wavelet='db1', mode='soft', level=1):    # 计算多级小波变换的系数    coeffs = pywt.wavedec(data, wavelet, mode='symmetric', level=level)    # 计算阈值    sigma = np.median(np.abs(coeffs[-1])) / 0.6745    threshold = sigma * np.sqrt(2 * np.log(len(data)))    # 固定阈值去噪    coeffs[1:] = (pywt.threshold(i, value=threshold, mode=mode) for i in coeffs[1:])    # 重构信号    return pywt.waverec(coeffs, wavelet, mode='symmetric')
2.2 无偏似然估计原则与python实现
      无偏似然估计方法是一种自适应阈值选择方法,其基于最小化估计误差的原理,通过估计原始信号和去噪信号之间的均方误差来自动选择最优阈值。该方法的核心思想是通过最小化一个无偏风险估计量来选择阈值,这个估计量考虑了去噪后信号与原始信号之间的误差。对于小波系数,该方法提供了一种计算阈值的方式,以期达到最小的均方误差。
      对于离散信号    ,构造新序列s
  
      定义风险集 合    ,计算各元素的风险系数:
  
设    的最小值为    ,该值对应的系数为    ,则阈值为
  
      无偏似然估计阈值小波降噪的python实现如下:























import numpy as npimport pywt
def sure_threshold(coeffs):    n = len(coeffs)    risks = np.zeros(n)    for i in range(n):        threshold = coeffs[i]        # 计算硬阈值处理后的风险        risks[i] = n - 2 * np.sum(coeffs > threshold) + np.sum((coeffs - threshold)**2)    # 选择最小风险对应的阈值    min_risk_index = np.argmin(risks)    return coeffs[min_risk_index]
def denoised_with_pywt_sure(data, wavelet='db1', mode='soft', level=1):    # 计算多级小波变换的系数    coeffs = pywt.wavedec(data, wavelet, mode='symmetric', level=level)    # 对每一层系数应用SURE方法计算阈值并去噪    for i in range(1, len(coeffs)):        threshold = sure_threshold(np.abs(coeffs[i]))        coeffs[i] = pywt.threshold(coeffs[i], value=threshold, mode=mode)    # 重构信号    return pywt.waverec(coeffs, wavelet, mode='symmetric')

2.3 极值阈值原则与python实现
极值阈值原则的原理是使估计的最大风险最小化,该阈值方法可实现最大均方误差的最小化。其数学表达如下:
  
      极值阈值小波降噪的python实现如下:
















import pywtimport numpy as np
def denoised_with_pywt_minimax(data, wavelet='db1', mode='soft', level=1):    # 计算多级小波变换的系数    coeffs = pywt.wavedec(data, wavelet, mode='symmetric', level=level)    # 计算阈值    sigma = np.median(np.abs(coeffs[-1])) / 0.6745    if len(data) >= 32:        threshold = sigma * (0.3936+(0.1829*np.log(len(data))/np.log(2)))    else:        threshold = 0    # 固定阈值去噪    coeffs[1:] = (pywt.threshold(i, value=threshold, mode=mode) for i in coeffs[1:])    # 重构信号    return pywt.waverec(coeffs, wavelet, mode='symmetric')

2.4 启发式阈值原则
启发式阈值原则是一种基于经验的阈值选择方法,它旨在找到一个合适的阈值来最大程度地去除噪声,同时保留信号的重要特征。这种方法通常依赖于对特定应用领域的理解和经验判断,而不是严格的数学准则。一种常用的启发式阈值原则将固定阈值和无偏似然估计阈值结合,若信噪比较小,采用无偏似然估计原则,否则使用固定阈值原则。

3 阈值函数介绍

3.1 硬阈值方法
      硬阈值函数的数学表达为:
  
      硬阈值法只保留大于阈值的小波系数,而将其他小波系数置为0。
3.2 软阈值方法
      软阈值函数的数学表达为:
  
      软阈值将小于阈值的系数置为0,将其余系数向0收缩。
       软硬阈值函数如图2所示。
图2 软硬阈值函数

4 算法测试

4.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_data_universal = denoised_with_pywt_universal(y, wavelet='db4', mode='soft', level=8)denoised_data_sure = denoised_with_pywt_sure(y, wavelet='db4', mode='soft', level=8)denoised_data_minimax = denoised_with_pywt_minimax(y, wavelet='db4', mode='soft', level=8)value_universal = denoise_evaluate(s, r, denoised_data_universal)value_sure = denoise_evaluate(s, r, denoised_data_sure)value_minimax = denoise_evaluate(s, r, denoised_data_minimax)
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=16)
plt.subplot(5, 1, 2)plt.plot(t, y, linewidth=1.5, color='b')plt.title('noised signal', fontname='Times New Roman', fontsize=16)
plt.subplot(5, 1, 3)plt.plot(t, denoised_data_universal, linewidth=1.5, color='b')plt.title('固定阈值法', fontname='Simsun', fontsize=16)
plt.subplot(5, 1, 4)plt.plot(t, denoised_data_sure, linewidth=1.5, color='b')plt.title('无偏似然估计阈值', fontname='Simsun', fontsize=16)
plt.subplot(5, 1, 5)plt.plot(denoised_data_minimax, linewidth=1.5, color='b')plt.title('极值阈值', fontname='Simsun', fontsize=16plt.show()
4.2 降噪结果
降噪结果如图3所示。
图3 降噪结果
4.3 降噪指标
降噪指标如表1所示。
表1 降噪指标

5 参考文献

[1]胡先伟. 小波降噪算法及IP软核实现技术[D]. 陕西:西安电子科技大学,2013. DOI:10.7666/d.Y2379886.
[2] 陈侨. 小波降噪在结构损伤识别中的应用[D]. 湖北:华中科技大学,2015. DOI:10.7666/d.D736262.
[3] 刘尊民. 小波降噪和时空轨迹数据精细化理论及在采油集输监控系统中的应用[D]. 山东:青岛理工大学,2019.


来源:故障诊断与python学习
电子pythonUM理论数字孪生
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-03-19
最近编辑:9月前
故障诊断与python学习
硕士 签名征集中
获赞 72粉丝 76文章 152课程 0
点赞
收藏
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈