【代码学习】使用DTCWT进行信号分解并计算分量的能量熵值

最新推荐文章于 2024-09-04 07:14:47 发布

王摇摆

最新推荐文章于 2024-09-04 07:14:47 发布

阅读量279

点赞数 3

文章标签：学习

本文链接：https://blog.csdn.net/weixin_44943389/article/details/140223382

版权

"""
# 使用双数复小波变换进行信号分解，得到分量，将分量进行可视化，同时计算每个分量的能量熵值，作为能量特征来进行模型输入
# -*- coding: utf-8 -*-
# @Time    : 2024/7/6 8:20
# @Author  : 王摇摆
# @FileName: test01.py
# @Software: PyCharm
# @Blog    ：https://blog.csdn.net/weixin_44943389?type=blog
"""
import numpy as np
import matplotlib.pyplot as plt
import dtcwt

def compute_energy_entropy(coeffs):
    energy = np.array([np.sum(np.abs(c)**2) for c in coeffs])
    total_energy = np.sum(energy)
    normalized_energy = energy / total_energy
    entropy = -np.sum(normalized_energy * np.log(normalized_energy + np.finfo(float).eps))
    return entropy, energy, normalized_energy

# 生成带噪声的信号
np.random.seed(0)
t = np.linspace(0, 1, 1024)
signal = np.sin(2 * np.pi * 10 * t) + np.sin(2 * np.pi * 20 * t) + 0.5 * np.random.normal(size=t.shape)

# 执行DTCWT
transform = dtcwt.Transform1d()
coeffs = transform.forward(signal, nlevels=3)

# 计算每个子带的能量熵
entropies = []
for i in range(len(coeffs.highpasses)):
    entropy, energy, normalized_energy = compute_energy_entropy(coeffs.highpasses[i])
    entropies.append(entropy)
    print(f"Level {i+1} Energy Entropy: {entropy}")

# 可视化原始信号
plt.figure(figsize=(12, 12))
plt.subplot(4, 2, 1)
plt.plot(t, signal)
plt.title('Original Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')

# 可视化每个分量
subplot_idx = 2
for i in range(3):  # 假设分解3层
    num_directions = coeffs.highpasses[i].shape[1]  # 获取每层的方向数量
    for j in range(num_directions):
        plt.subplot(4, 2, subplot_idx)
        plt.plot(np.real(coeffs.highpasses[i][:, j]))
        plt.title(f'Level {i+1} - Direction {j+1}')
        plt.xlabel('Time')
        plt.ylabel('Amplitude')
        subplot_idx += 1

plt.tight_layout()
plt.show()

# 可视化每个子带的能量熵
plt.figure(figsize=(10, 4))
plt.bar(range(len(entropies)), entropies)
plt.title('Energy Entropy of Each Subband')
plt.xlabel('Subband')
plt.ylabel('Energy Entropy')
plt.show()