信号处理基础（2）小波阈值去噪方法-CSDN博客

本文链接：https://blog.csdn.net/m0_63478459/article/details/147922180

上篇文章主要介绍了小波阈值去噪的原理：信号处理基础（1）小波阈值去噪原理详解-CSDN博客接下来讲下如何使用python进行实现

一、小波阈值去噪的关键步骤

小波分解：选择合适的小波基和分解层数，将信号分解为近似系数和细节系数。
阈值处理：确定阈值并选择阈值函数对细节系数进行处理。
小波重构：利用处理后的系数重构信号。

二、Python 实现小波阈值去噪

主要分为两部分：小波分解与重构、阈值函数实现

1.小波分解与重构

import numpy as np
import pywt
import matplotlib.pyplot as plt

def wavelet_decomposition(signal, wavelet='db4', level=3):
    """
    小波分解函数
    Parameters:
        signal (array): 输入信号
        wavelet (str): 小波基，默认为'db4'
        level (int): 分解层数
    Returns:
        coeffs (list): 分解后的系数列表
        sigma (float): 噪声标准差估计
    """
    # 小波分解
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    
    # 计算噪声估计 (使用最高频细节系数)
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745  # 噪声标准差估计
    
    return coeffs, sigma

def wavelet_reconstruction(processed_coeffs, wavelet='db4', original_length=None):
    """
    小波重构函数
    Parameters:
        processed_coeffs (list): 处理后的系数列表
        wavelet (str): 小波基，这里使用'db4'
        original_length (int): 原始信号长度，用于处理边界效应
    Returns:
        denoised_signal (array): 重构后的信号
    """
    # 小波重构
    denoised_signal = pywt.waverec(processed_coeffs, wavelet)
    
    # 处理边界效应（可能长度不一致）
    if original_length is not None:
        return denoised_signal[:original_length]
    return denoised_signal

2.阈值函数实现部分：硬阈值、软阈值和 Garrote 阈值

def hard_threshold_denoising(coeffs, sigma, threshold_mode='universal'):
    """
    硬阈值去噪函数
    Parameters:
        coeffs (list): 分解后的系数列表
        sigma (float): 噪声标准差估计
        threshold_mode (str): 阈值计算方式 ('universal', 'manual')
    Returns:
        processed_coeffs (list): 处理后的系数列表
    """
    # 计算通用阈值
    if threshold_mode == 'universal':
        N = len(np.concatenate(coeffs))
        threshold = sigma * np.sqrt(2 * np.log(N))
    else:
        threshold = sigma  # 可自定义阈值
    
    # 应用硬阈值处理（不处理近似系数）
    processed_coeffs = [coeffs[0]]  # 保留近似系数
    for i in range(1, len(coeffs)):
        processed_coeffs.append(coeffs[i] * (np.abs(coeffs[i]) >= threshold))
    
    return processed_coeffs

def soft_threshold_denoising(coeffs, sigma, threshold_mode='universal'):
    """
    软阈值去噪函数
    Parameters:
        coeffs (list): 分解后的系数列表
        sigma (float): 噪声标准差估计
        threshold_mode (str): 阈值计算方式 ('universal', 'manual')
    Returns:
        processed_coeffs (list): 处理后的系数列表
    """
    # 计算通用阈值
    if threshold_mode == 'universal':
        N = len(np.concatenate(coeffs))
        threshold = sigma * np.sqrt(2 * np.log(N))
    else:
        threshold = sigma  # 可自定义阈值
    
    # 应用软阈值处理（不处理近似系数）
    processed_coeffs = [coeffs[0]]  # 保留近似系数
    for i in range(1, len(coeffs)):
        processed_coeffs.append(np.sign(coeffs[i]) * np.maximum(np.abs(coeffs[i]) - threshold, 0))
    
    return processed_coeffs

def garrote_threshold_denoising(coeffs, sigma, threshold_mode='universal'):
    """
    Garrote阈值去噪函数
    Parameters:
        coeffs (list): 分解后的系数列表
        sigma (float): 噪声标准差估计
        threshold_mode (str): 阈值计算方式 ('universal', 'manual')
    Returns:
        processed_coeffs (list): 处理后的系数列表
    """
    # 计算通用阈值
    if threshold_mode == 'universal':
        N = len(np.concatenate(coeffs))
        threshold = sigma * np.sqrt(2 * np.log(N))
    else:
        threshold = sigma  # 可自定义阈值
    
    # 应用Garrote阈值处理（不处理近似系数）
    processed_coeffs = [coeffs[0]]  # 保留近似系数
    for i in range(1, len(coeffs)):
        abs_coeffs = np.abs(coeffs[i])
        processed_coeffs.append(np.where(abs_coeffs >= threshold, 
                                        (1 - (threshold**2)/(abs_coeffs**2)) * coeffs[i], 0))
    
    return processed_coeffs

三、参数的选择与分析

小波基选择：代码中使用了 Daubechies 4 (db4) 小波基，这是一种常用的紧支集小波基。不同的小波基对去噪效果有一定影响，实际应用中需要根据信号特性选择合适的小波基。
分解层数确定：分解层数level决定了信号被分解的精细程度。层数越多，细节信息越丰富，但计算复杂度也会增加。一般根据信号长度和噪声特性来选择。
阈值计算方法：
- Universal 阈值：基于 Donoho 提出的通用阈值公式 λ=σ√(2logN)
- 手动阈值：根据实际情况自定义阈值
阈值函数比较：
- 硬阈值：直接将小于阈值的系数置零，能较好保留信号的边缘和细节，但可能产生振荡。
- 软阈值：将大于阈值的系数进行收缩，去噪后的信号更加平滑，但可能导致边缘模糊。
- Garrote 阈值：介于硬阈值和软阈值之间的折衷方法，既保留了信号的细节，又避免了硬阈值的振荡问题。