小波包变换

最新推荐文章于 2025-06-26 15:10:50 发布

JoinApper

最新推荐文章于 2025-06-26 15:10:50 发布

阅读量719

点赞数 5

CC 4.0 BY-SA版权

分类专栏：故障诊断信号与系统文章标签：信号处理

本文链接：https://blog.csdn.net/qq1234534215/article/details/145899087

故障诊断同时被 2 个专栏收录

10 篇文章

订阅专栏

信号与系统

9 篇文章

订阅专栏

前言

之前曾介绍了离散小波变换，本质上是通过小波基函数构造的低通滤波器和高通滤波器，将信号分解为近似系数（低频部分）和细节系数（高频部分），经过下采样后将低频部分继续分解，直至达到指定的分解层数。但这样一来，高频部分就无法做得更精细的分析了，而小波包变换则在此基础上进行了优化，对低频和高频部分都进行分解，最终形成一棵小波包树。

小波包变换原理

下面我们就来具体看看小波包变换的原理，信号通过低通滤波器和高通滤波器，再进行下采样后就可以得到近似系数(cA)和细节系数(cD)，这步骤和离散小波变换是一致的：
在这里插入图片描述
然后将每一层的近似系数和细节系数，再通过低通滤波器和高通滤波器以及下采样，进行分解，直到满足预设的分解层数即可：

Python示例

我们还是以一个实际例子来说明小波包变换的使用，仍然使用pywt这个包。
首先生成一个测试信号：

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

plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）

fs = 1000  # 采样率 (Hz)
t = np.arange(0, 1, 1/fs)
f1, f2, f3 = 50, 150, 250  # 频率分量
signal = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t) + 0.3 * np.sin(2 * np.pi * f3 * t)

在这里插入图片描述
然后可以通过小波包进行测试信号的分解：

# 小波包分解参数
wavelet = 'db4'  # 小波基函数
level = 3       # 分解层数

# 创建小波包树
wp = pywt.WaveletPacket(signal, wavelet, mode='symmetric', maxlevel=level)

# 获取所有叶子节点路径（最后一层的所有子带）
leaf_paths = [node.path for node in wp.get_level(level, 'natural')]

# 计算每个子带的能量（系数平方和）
energies = []
for path in leaf_paths:
    node = wp[path]
    coeffs = node.data
    energy = np.sum(coeffs**2)
    energies.append(energy)

# 计算每个子带的频率范围
freq_bands = []
for i in range(len(leaf_paths)):
    band_width = fs / (2 ** level)
    start_freq = i * band_width
    end_freq = (i + 1) * band_width
    freq_bands.append((start_freq, end_freq))

# 可视化
plt.figure(figsize=(12, 10))

# 原始信号
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title('测试信号')
plt.xlabel('时间')
plt.ylabel('幅值')

# 小波包子带能量分布
plt.subplot(2, 1, 2)
plt.bar(range(len(energies)), energies, tick_label=[f"{f[0]:.0f}-{f[1]:.0f} Hz" for f in freq_bands])
plt.title('小波包能量分布(Level 3)')
plt.xlabel('频带')
plt.ylabel('能量')
plt.xticks(rotation=45)

plt.tight_layout()
plt.show()

在这里插入图片描述
为了更容易理解小波包变换的分解过程，我们也可以手动实现：

def extend_signal(x, filter_len):
    # 对称延拓
    return np.pad(x, (filter_len-1, filter_len-1), mode='symmetric')


def wavelet_packet_decomposition(x, level=3, wavelet='db4'):
    """
    手动实现小波包分解
    :param x: 输入信号(长度需为2^n)
    :param wavelet: 小波基
    :param level: 分解层数
    :return: 字典结构，键为节点路径，值为系数数组
    """
    # 初始化滤波器
    lo, hi, _, _ = pywt.Wavelet(wavelet).filter_bank

    # 初始化分解树
    wp_tree = {}
    wp_tree[''] = x  # 根节点路径为空字符串

    # 递归分解函数
    def decompose(node_path, signal, current_level):
        if current_level >= level:
            return

        # 延拓
        signal = extend_signal(signal, len(lo))
        # 执行分解
        cA = np.convolve(signal, lo, mode='valid')[1::2]  # 低通+下采样
        cD = np.convolve(signal, hi, mode='valid')[1::2]  # 高通+下采样

        # 存储系数
        wp_tree[node_path + 'a'] = cA
        wp_tree[node_path + 'd'] = cD

        # 递归分解两个子节点
        decompose(node_path + 'a', cA, current_level + 1)
        decompose(node_path + 'd', cD, current_level + 1)

    # 开始递归分解
    decompose('', x, 0)
    return wp_tree


# 执行3层小波包分解
wp_tree = wavelet_packet_decomposition(signal, level=3)

# 计算每层能量
energies_manual = []
for key, value in wp_tree.items():
    if len(key) < 3:
        continue
    energy = np.sum(value ** 2)
    energies_manual.append(energy)