Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”_信号emd后生成小波时频图

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
1
2
3
4
5
6
7
8

  原始信号时域波形图：
  原始信号频谱图，可以看出频率成分包括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])     #画希尔伯特谱
    
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

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