Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”

Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”  信号处理中常需要分析时域统计量、频率成分,但不平稳信号的时域波形往往复杂、无序,且傅里叶变换得到的频率成分是该时间段内的平均频率,无法分析频率随时间变化的情况。随后,短时傅里叶变换(STFT)、小波变换(WT)、希尔伯特变换(HHT)等时频分析方法相继而出。  其中,STFT受时间窗口的影响、WT则需要自己选择小波、HHT在变换时需要预先将信号分解为平稳信号。由于网上只有CWT小波时频图的python代码,笔者自

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全家桶1年46,售后保障稳定

Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”

  信号处理中常需要分析时域统计量、频率成分,但不平稳信号的时域波形往往复杂、无序,且傅里叶变换得到的频率成分是该时间段内的平均频率,无法分析频率随时间变化的情况。随后,短时傅里叶变换(STFT)、小波变换(WT)、希尔伯特变换(HHT)等时频分析方法相继而出。
  其中,STFT受时间窗口的影响、WT则需要自己选择小波、HHT在变换时需要预先将信号分解为平稳信号。由于网上只有CWT小波时频图的python代码,笔者自编了不同分解算法+Hilbert时频图的代码与其比较。

  工作不易,转载请引用!!!

仿真信号

  构造频率随时间变化的信号如下:

#构造测试信号
import nump as np
Fs=1000   #采样频率
t = np.arange(0, 1.0, 1.0 / Fs)
f1,f2,f3 = 100,200,300
signal = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
                    [lambda t: np.sin(2 * np.pi * f1 * t), lambda t: np.sin(2 * np.pi * f2 * t),
                     lambda t: np.sin(2 * np.pi * f3 * t)])                       #仿真信号1

Jetbrains全家桶1年46,售后保障稳定

  原始信号时域波形图:
在这里插入图片描述  原始信号频谱图,可以看出频率成分包括100、200、300Hz,但不能得知这些频率在何时刻出现:
在这里插入图片描述

1、EMD\EEMD\VMD分解+Hilbert时频图

  经验模态分解(EMD)由Hilbert提出,目的在于将不平稳信号分解为各平稳的IMF分量,但其“端点效应”与“模态混叠”缺点较突出。在其基础上,集成经验模态分解(EEMD)在EMD分解前加不同的高斯白噪声,一定程度上抑制了“模态混叠”,但增加了计算成本。变分模态分解(VMD)可以实现信号频域内各个分量的自适应分割,但需要指定模态个数K等参数。具体原理可自行补习。
  本文全部代码基于python 3.9,EMD\EEMD分解采用的是 PyEMD工具包(注意大小写!),VMD分解采用的是GitHub上的vmdpy代码,fftlw是笔者之前博文写的快速傅里叶变化代码,请自行下载。EMD\EEMD\VMD分解+Hilbert时频图的函数代码如下,其中,只需在调用decompose_lw()时改method即可以换不同的分解方法:

# -*- coding: utf-8 -*-
""" Created on Fri Dec 17 21:18:48 2021 @author: lw """
import matplotlib.pyplot as plt
import numpy as np
from PyEMD import EEMD,EMD,Visualisation
from scipy.signal import hilbert
# from fftlw import fftlw
from vmdpy import VMD
#分解方法(emd、eemd、vmd)
def decompose_lw(signal,t,method='eemd',K=5,draw=1):
names=['emd','eemd','vmd']
idx=names.index(method)
#emd分解
if idx==0:
emd = EMD()
IMFs= emd.emd(signal)
#vmd分解
elif idx==2:
alpha = 2000       # moderate bandwidth constraint
tau = 0.            # noise-tolerance (no strict fidelity enforcement)
DC = 0             # no DC part imposed
init = 1           # initialize omegas uniformly
tol = 1e-7
# Run actual VMD code
IMFs, _, _ = VMD(signal, alpha, tau, K, DC, init, tol)
#eemd分解
else:
eemd = EEMD()
emd = eemd.EMD
emd.extrema_detection="parabol"
IMFs= eemd.eemd(signal,t)
#可视化
if draw==1:
plt.figure()
for i in range(len(IMFs)):
plt.subplot(len(IMFs),1,i+1)
plt.plot(t,IMFs[i])
if i==0:
plt.rcParams['font.sans-serif']='Times New Roman'
plt.title('Decomposition Signal',fontsize=14)
elif i==len(IMFs)-1:
plt.rcParams['font.sans-serif']='Times New Roman'
plt.xlabel('Time/s')
# plt.tight_layout()
return IMFs       
#希尔波特变换及画时频谱
def hhtlw(IMFs,t,f_range=[0,500],t_range=[0,1],ft_size=[128,128],draw=1):
fmin,fmax=f_range[0],f_range[1]         #时频图所展示的频率范围
tmin,tmax=t_range[0],t_range[1]         #时间范围
fdim,tdim=ft_size[0],ft_size[1]         #时频图的尺寸(分辨率)
dt=(tmax-tmin)/(tdim-1)
df=(fmax-fmin)/(fdim-1)
vis = Visualisation()
#希尔伯特变化
c_matrix=np.zeros((fdim,tdim))
for imf in IMFs:
imf=np.array([imf])
#求瞬时频率
freqs = abs(vis._calc_inst_freq(imf, t, order=False, alpha=None))
#求瞬时幅值
amp= abs(hilbert(imf))
#去掉为1的维度
freqs=np.squeeze(freqs)
amp=np.squeeze(amp)
#转换成矩阵
temp_matrix=np.zeros((fdim,tdim))
n_matrix=np.zeros((fdim,tdim))
for i,j,k in zip(t,freqs,amp):
if i>=tmin and i<=tmax and j>=fmin and j<=fmax:
temp_matrix[round((j-fmin)/df)][round((i-tmin)/dt)]+=k
n_matrix[round((j-fmin)/df)][round((i-tmin)/dt)]+=1
n_matrix=n_matrix.reshape(-1)
idx=np.where(n_matrix==0)[0]
n_matrix[idx]=1
n_matrix=n_matrix.reshape(fdim,tdim)
temp_matrix=temp_matrix/n_matrix
c_matrix+=temp_matrix
t=np.linspace(tmin,tmax,tdim)
f=np.linspace(fmin,fmax,fdim)
#可视化
if draw==1:
fig,axes=plt.subplots()
plt.rcParams['font.sans-serif']='Times New Roman'
plt.contourf(t, f, c_matrix,cmap="jet")
plt.xlabel('Time/s',fontsize=16)
plt.ylabel('Frequency/Hz',fontsize=16)
plt.title('Hilbert spectrum',fontsize=20)
x_labels=axes.get_xticklabels()
[label.set_fontname('Times New Roman') for label in x_labels]
y_labels=axes.get_yticklabels()
[label.set_fontname('Times New Roman') for label in y_labels]
# plt.show()
return t,f,c_matrix
#%%测试函数
if __name__=='__main__':
#构造测试信号
Fs=1000   #采样频率
t = np.arange(0, 1.0, 1.0 / Fs)
f1,f2,f3 = 100,200,300
signal = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
[lambda t: np.sin(2 * np.pi * f1 * t), lambda t: np.sin(2 * np.pi * f2 * t),
lambda t: np.sin(2 * np.pi * f3 * t)])                       #仿真信号1
# signal = 3*np.sin(2*np.pi*f1*t)+6*np.sin(2*np.pi*f2*t)+5*np.sin(2*np.pi*f3*t) #仿真信号2
# signal = 3*t*np.sin(2*np.pi*f1*t) #仿真信号3
#画时域图
plt.figure()
plt.plot(t,signal)
plt.rcParams['font.sans-serif']='Times New Roman'
plt.xlabel('Time/s',fontsize=16)
plt.title('Original Signal',fontsize=20) 
plt.show()
#画仿真信号频谱图
# _,_=fftlw(Fs,signal,1)
IMFs=decompose_lw(signal,t,method='vmd',K=10)                                    #分解信号
tt,ff,c_matrix=hhtlw(IMFs,t,f_range=[0,500],t_range=[0,1],ft_size=[128,128])     #画希尔伯特谱

a、EMD分解+Hilbert时频图

  EMD分解所得IMF分量,可知分量1存在模态混叠现象:
在这里插入图片描述  时频图:
在这里插入图片描述

b、EEMD分解+Hilbert时频图

  EEMD分解所得IMF分量,可知分量1仍然存在模态混叠现象:
在这里插入图片描述  时频图,相比EMD,300Hz更加集中:
在这里插入图片描述

c、VMD分解+Hilbert时频图

  VMD分解所得IMF分量,默认模态个数设为10,可知基本能准确分出100、200、300Hz分量,但还存在端点效应:

在这里插入图片描述  时频图,频率成分更加集中,效果更好:
在这里插入图片描述

2、CWT小波时频图

  连续小波时频图是转载自知乎文章

  连续小波变换(CWT)时频图绘制 python实现

# -*- coding: utf-8 -*-
""" Created on Fri Dec 17 20:17:42 2021 @author: lw """
import matplotlib.pyplot as plt
import numpy as np
import pywt
# from matplotlib.font_manager import FontProperties
# chinese_font = FontProperties(fname='/usr/share/fonts/truetype/wqy/wqy-microhei.ttc')
sampling_rate = 1000
t = np.arange(0, 1.0, 1.0 / sampling_rate)
f1 = 100
f2 = 200
f3 = 300
data = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
[lambda t: np.sin(2 * np.pi * f1 * t), lambda t: np.sin(2 * np.pi * f2 * t),
lambda t: np.sin(2 * np.pi * f3 * t)])
wavename = 'cgau8'
totalscal = 256
fc = pywt.central_frequency(wavename)
cparam = 2 * fc * totalscal
scales = cparam / np.arange(totalscal, 1, -1)
[cwtmatr, frequencies] = pywt.cwt(data, scales, wavename, 1.0 / sampling_rate)
plt.figure(figsize=(8, 8))
plt.subplot(211)
plt.plot(t, data)
# plt.xlabel(u"时间(秒)", fontproperties=chinese_font)
# plt.title(u"300Hz和200Hz和100Hz的分段波形和时频谱", fontproperties=chinese_font, fontsize=20)
plt.subplot(212)
plt.contourf(t, frequencies, abs(cwtmatr))
# plt.ylabel(u"频率(Hz)", fontproperties=chinese_font)
# plt.xlabel(u"时间(秒)", fontproperties=chinese_font)
plt.subplots_adjust(hspace=0.4)
plt.tight_layout()
plt.show()

  原始信号与时频图如下,可知小波分解的频率分量比较分散,且幅值(颜色)有点对应不上,理论上100、200、300Hz的颜色应一样亮:

在这里插入图片描述

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/231289.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • Navicat 15 for MySQL激活码【2021最新】

    (Navicat 15 for MySQL激活码)JetBrains旗下有多款编译器工具(如:IntelliJ、WebStorm、PyCharm等)在各编程领域几乎都占据了垄断地位。建立在开源IntelliJ平台之上,过去15年以来,JetBrains一直在不断发展和完善这个平台。这个平台可以针对您的开发工作流进行微调并且能够提供…

  • C语言括号匹配(栈括号匹配c语言)

    给定一串字符,不超过100个字符,可能包括括号、数字、字母、标点符号、空格,编程检查这一串字符中的(),[],{}是否匹配。输入格式:输入在一行中给出一行字符串,不超过100个字符,可能包括括号、数字、字母、标点符号、空格。输出格式:如果括号配对,输出yes,否则输出no。输入样例1:sin(10+20)输出样例1:yes输入样例2:{[}]输出样例2:no思路:题目输入一些字符串,我们就先保留括号之类的,判断是否匹配。如果遇到左括号,就入栈,如果遇到一个右括号,就与栈顶元

  • paceMaker_pacemaker怎么读

    paceMaker_pacemaker怎么读1. 简介 Pacemaker是一个集群资源管理者。他用资源级别的监测和恢复来保证集群服务(aka.资源)的最大可用性。它可以用你所擅长的基础组件(Corosync或者是Heartbeat)来实现通信和关系管理。​2. 特性 Pacemaker包含以下的关键特性:  监测并恢复节点和服务级别的故障​  存储无关,并不需要共享存储​  资源无关,任何能用脚本控制的资源

    2022年10月25日
  • windows本地用户及组的区别「建议收藏」

    windows本地用户及组的区别「建议收藏」Administrators(超级管理员组)用来管理注册用户或者具有一定权限的其他管理员,维护网站的运行。Administrators中的用户对计算机/域有不受限制的完全访问权,分配给该组的默认权

  • vue 富文本编辑框_基于vue的富文本编辑器

    vue 富文本编辑框_基于vue的富文本编辑器1、下载插件npmiwangeditor–save插件官网地址:https://www.wangeditor.com/2、封装富文本组件<templatelang=”html”><divclass=”editor”><!–<divref=”toolbar”class=”toolbar”></div>–><divref=”editor”class=”text”></div

  • 物联网实践

    物联网实践

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号