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

大家好，欢迎来到IT知识分享网。

 信号处理中常需要分析时域统计量、频率成分，但不平稳信号的时域波形往往复杂、无序，且傅里叶变换得到的频率成分是该时间段内的平均频率，无法分析频率随时间变化的情况。随后，短时傅里叶变换（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 

  原始信号时域波形图：

  原始信号频谱图，可以看出频率成分包括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]) #画希尔伯特谱 000

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() 8

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