Python环境下一种基于改进小波变换的信号时频分析方法

发布时间:2024年01月19日

小波,不严格的说即有限持续时间的波形,其平均值为零。许多我们感兴趣的信号和图像表现出瞬变行为。例如语音信号的特点是辅音短脉冲编码,然后元音稳态振荡;自然图像边缘突变;金融时间序列表现出瞬态行为,经济状况的快速上升和下降。与傅里叶基不同,小波基比较擅长稀疏表示分段规则信号和图像,其中包括众多瞬态行为。将小波与正弦波进行比较,正弦波是傅里叶分析的基础,正弦曲线的持续时间没有限制—负无穷延伸到正无穷。正弦曲线是平滑的,小波往往是不规则和不对称的。

针对小波变换存在的依赖母函数选择问题,受多分辨率思想的启发,对小波进行了修正,以一系列小波集代替传统小波母函数,对小波进行阶数选择及循环数缩放,得到能够在时频面有着较高分辨率的小波变换结果。

需要提前安装JAX模块,可以参考:

https://blog.csdn.net/weixin_50008473/article/details/126589113

程序部分代码如下:

import sys
sys.path.insert(0, '..')

import jax.numpy as jnp
import matplotlib.pyplot as plt
from improvewavelets import wavelet_transform, adaptive_wavelet_transform

fs = 1024
burst_freqs = [20, 40, 60]
f_shift = 10
n_cycles = 11
n_neighb_cycles = 12

ys = []

# create a 0.1s blank signal to start
ys.append(jnp.zeros(int(fs*0.1)))

for f in burst_freqs:
  # frequency contaminated signal
  t = 1/f * n_cycles
  x = jnp.linspace(0, t, int(t * fs))
  y = jnp.sin(2*jnp.pi*f*x) + jnp.sin(2*jnp.pi*(f+f_shift)*x - jnp.pi/1.5)
  ys.append(y)

  # time contaminated signal, 2 cycles later
  ys.append(jnp.zeros(int(fs*(1/f)*2)))
  
  t2 = 1/f * n_neighb_cycles
  x = jnp.linspace(0, t2, int(t2 * fs))
  y = jnp.sin(2*jnp.pi*f*x)
  ys.append(y)

  # space between bursts of 0.1s
  ys.append(jnp.zeros(int(fs*0.1)))

signal = jnp.concatenate(ys)


fig, ax = plt.subplots(figsize=(15, 2), dpi=300)
ax.set_xlabel("Time (ms)")
ax.plot(jnp.linspace(0, len(signal)/fs, len(signal)), signal)

freqs = jnp.linspace(10, 80, 141)

fig, ax = plt.subplots(ncols=3, figsize=(15, 5), dpi=300)

for i, c in enumerate([3, 16, 33]):
    scalogram = wavelet_transform(signal, freqs, c, fs)
    ax[i].imshow(jnp.abs(scalogram)**2, aspect=1/40, cmap="jet", interpolation="none", origin="lower", extent=[0, len(signal)/fs, freqs[0], freqs[-1]])
    ax[i].set_title(f"Cycles: {c}")
    ax[i].set_xlabel("Time (s)")
    ax[i].set_ylabel("Frequency (Hz)")

fig, ax = plt.subplots(ncols=3, figsize=(15, 5), dpi=300)

for (i, (base_cycle, min_order, max_order)) in enumerate(zip([3, 5, 1], [1, 1, 5], [30, 30, 40])):
    scalogram = adaptive_wavelet_transform(signal, freqs, sampling_freq=fs, 
                                        base_cycle=base_cycle, min_order=min_order, max_order=max_order, mode="add")
    ax[i].imshow(jnp.abs(scalogram)**2, aspect=1/40, cmap="jet", interpolation="none", origin="lower", extent=[0, len(signal)/fs, freqs[0], freqs[-1]])
    ax[i].set_title(f"Base cycles: {base_cycle}, Orders: {min_order}-{max_order}")
    ax[i].set_xlabel("Time (s)")
    ax[i].set_ylabel("Frequency (Hz)")

出图如下:

工学博士,担任《Mechanical System and Signal Processing》审稿专家,担任
《中国电机工程学报》优秀审稿专家,《控制与决策》,《系统工程与电子技术》,《电力系统保护与控制》,《宇航学报》等EI期刊审稿专家,担任《计算机科学》,《电子器件》 , 《现代制造过程》 ,《电源学报》,《船舶工程》 ,《轴承》 ,《工矿自动化》 ,《重庆理工大学学报》 ,《噪声与振动控制》 ,《机械传动》 ,《机械强度》 ,《机械科学与技术》 ,《机床与液压》,《声学技术》,《应用声学》,《石油机械》,《西安工业大学学报》等中文核心审稿专家。
擅长领域:现代信号处理,机器学习,深度学习,数字孪生,时间序列分析,设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。

文章来源:https://blog.csdn.net/weixin_39402231/article/details/135694864
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。