Python中正弦信号的频域与时频域特征提取及小波包计算

用户生成的10秒正弦波代码 #

# Import libraries
import numpy as np
import matplotlib.pyplot as plt

Fs  = 1000                                              # Sampling frequency (Hz)
T   = 10                                                # Duration of time signal (s)
N   = T*Fs                                              # Number of data points
f   = 10                                                # Frequency of the sine wave signal (Hz)
t   = np.linspace(0, T, N)                              # Time axis data points
y   = 10 * np.sin(2*np.pi*f*t)                          # Amplitude of the sine wave signal

plt.figure(figsize=(8,8))                               # Figure size of the plot
plt.plot(t,y)                                           # time history plot of the sine wave
plt.xlabel('Time (s)')                                  # X-axis label
plt.ylabel('Amplitude')                                 # Y-axis label
##plt.xlim([0,2.5])                                     # X-axis limits
##plt.ylim([-8,18])                                     # Y-axis limits
plt.get_current_fig_manager().window.state('zoomed')    # Maximize the plot for full screen
plt.show()

1. 可提取的频域特征 #

对于这个纯10Hz正弦波信号，可提取的频域特征包括：

• 峰值频率：频谱中幅度最大值对应的频率，即信号本身的10Hz，是信号的主频率成分。
• 频谱峰值幅度：10Hz频率点对应的频谱幅度值（通过FFT计算时，实正弦信号的正频率分量幅度为(N/2)*10，其中N是采样点数）。
• 总功率/能量：信号的平均功率为(幅度^2)/2 = 50，频域上可通过对频谱幅度平方积分后除以采样频率得到。
• 谐波含量：该信号为纯正弦波，不含任何谐波成分，谐波含量为0。
• 频率分辨率：由采样频率和采样点数决定，计算公式为Fs/N = 0.1Hz，代表频谱中相邻两个频率点的间隔。

2. 可提取的时频域特征 #

作为平稳正弦信号，其时频域特征主要体现为频率不随时间变化，具体包括：

• 短时傅里叶变换（STFT）时频图：时频图中会呈现一条水平直线，集中在10Hz频率位置，所有时间点的能量都分布在该频率上。
• 小波尺度图（Scalogram）：小波变换的结果会显示在对应10Hz的尺度位置上，能量集中且不随时间波动。
• 瞬时频率：信号的瞬时频率始终稳定在10Hz，无任何起伏。
• 时频能量分布：能量在时间轴上均匀分布，且全部集中在10Hz频率区间内，无跨频率的能量扩散。

3. Python中计算小波包能量与小波包熵 #

需要使用PyWavelets库（需先通过pip install pywavelets安装），具体实现步骤如下：

代码示例 #

import numpy as np
import pywt

# 沿用用户生成的信号参数
Fs = 1000
T = 10
N = T * Fs
t = np.linspace(0, T, N)
y = 10 * np.sin(2*np.pi*10*t)

# 1. 初始化小波包分解
wavelet = 'db4'  # 选择Daubechies4小波基
level = 3        # 设置分解层数
wp = pywt.WaveletPacket(data=y, wavelet=wavelet, mode='symmetric', maxlevel=level)

# 2. 获取所有分解节点的路径
nodes = [node.path for node in wp.get_level(level, 'natural')]

# 3. 计算每个节点的能量
node_energies = []
total_energy = 0
for path in nodes:
    coeffs = wp[path].data
    energy = np.sum(np.square(coeffs))
    node_energies.append(energy)
    total_energy += energy

# 4. 计算小波包能量占比
energy_ratios = [e / total_energy for e in node_energies]

# 5. 计算小波包Shannon熵
epsilon = 1e-10  # 避免0值导致log计算报错
shannon_entropy = -np.sum([ratio * np.log2(ratio + epsilon) for ratio in energy_ratios])

# 输出结果
print("各小波包节点能量：", node_energies)
print("总能量：", total_energy)
print("小波包Shannon熵：", shannon_entropy)

说明 #

• 小波包能量：每个分解节点的能量为该节点系数的平方和，总能量等于原信号的能量（符合Parseval定理）。
• 小波包熵：这里采用Shannon熵，通过计算各节点能量占总能量的比例代入熵公式得到。对于纯正弦波，能量会集中在少数几个节点，因此熵值较低。

内容的提问来源于stack exchange，提问作者NN_Developer