大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新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的颜色应一样亮：

在这里插入图片描述

发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/231289.html原文链接：https://javaforall.net

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

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

工作不易，转载请引用！！！

仿真信号

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

a、EMD分解+Hilbert时频图

b、EEMD分解+Hilbert时频图

c、VMD分解+Hilbert时频图

2、CWT小波时频图

相关推荐

使用Windows Live Writer 常见的几个Blog接口

ArcGIS二次开发基础教程(02)：地图导航和鹰眼「建议收藏」

delphi xe5 激活成功教程

通过主机名得到对应ip地址_如何查看电脑主机名或Ip地址

非主流文字生成_非主流文字转换器

二进制数的减法计算_二进制乘法计算器

发表回复