HHT方法的python代码完美实现——结构时变模态参数识别方法

针对非线性、非平稳的数据，基函数的形式无法提前确定，需要自适应基函数的方法。

21世纪初，黄鳄提出了基于经验的HHT方法，在工程信号分析领域有广泛的应用。

HHT主要包括两个步骤：

1. 经验模态分解（EMD）
2. 希尔伯特变换（Hilbert）

• 首先合成示例信号

# 生成0-1时间序列，共2048个点
N = 2048
t = np.linspace(0, 1, N)
# 生成信号
signal = (2+np.cos(8*np.pi*t)) * np.cos(40*np.pi*(t+1)**2) + np.cos(20*np.pi*t + 5*np.sin(2*np.pi*t))

• 画出该合成信号的时域图和频谱

•  进行HHT分析，画出每一个可能IMF的时域图和频谱

分析结果为前两阶为待分析信号的两个IMF。

后两阶与前两阶的幅值差异较大，且前两阶IMF信号的频谱图组合基本上就是原始信号的频谱图。

•  选择前2个IMF，画出其时域图和时频图

对前2个IMF进行希尔伯特变换得到的时频图。

•  验证前n个IMF能否合成原始信号

对前2阶IMF信号求和，画出他们的合成信号与原始信号进行对比，发现基本吻合。

下面给出的代码完美实现这一过程，用户只需要修改采样点数和信号函数即可方便使用。

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from pyhht import EMD
from scipy.signal import hilbert
import tftb.processing
matplotlib.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
matplotlib.rcParams['axes.unicode_minus']=False       #显示负号

def picture(x, y, N):
    '''
    画信号的时域图和频谱
    输入：
    x: 0-1时间序列
    y: 信号
    N: 1s内采样点数
    输出：
    信号的时域图和频谱
    '''
    ax1 = plt.subplot(2, 1, 1) 
    plt.plot(x, y)  
    plt.xlabel('时间/s') 
    plt.ylabel('幅值')
    plt.title('合成信号时域曲线')
    yf = np.fft.fft(y)
    xf = np.linspace(0.0, N/2, N//2)
    ax2 = plt.subplot(2, 1, 2) 
    plt.plot(xf, 2.0/N * np.abs(yf[:N//2])) #频谱幅值归一化，需要*2/N
    plt.xlabel('频率/Hz') 
    plt.ylabel('幅值')
    plt.title('合成信号频谱')
    plt.show() 

def HHTAnalysis(t, signal, N):
    '''
    进行HHT分析，画出每一个IMF的时域图和频谱
    输入：
    t: 0-1时间序列
    signal: 信号
    N: 1s内采样点数
    输出：
    信号每一个IMF的时域图和频谱
    且返回IMFs，维度为(IMFs个数，信号点数N)
    '''
    # 进行EMD分解
    decomposer = EMD(signal)
    # 获取EMD分解后的IMF成分
    imfs = decomposer.decompose()
    # 分解后的IMF个数
    n_components = imfs.shape[0]
    # 画出每一个IMF的时域图和频谱
    fig1, axes = plt.subplots(n_components, 2, figsize=(10, 7), sharex='col', sharey=False)
    for i in range(n_components):
        #画时域图
        axes[i][0].plot(t, imfs[i])
        axes[i][0].set_title('IMF{}'.format(i + 1))
        # 画fft图
        yf = np.fft.fft(imfs[i])   
        xf = np.linspace(0.0, N/2, N//2) 
        axes[i][1].plot(xf, 2.0/N * np.abs(yf[:N//2])) #频谱幅值归一化，需要*2/N
        axes[i][1].set_title('IMF{}'.format(i + 1))
    plt.show()   
    return imfs

def HHTPicture(t, imfs, N, n):
    '''
    画出指定个数的IMFs的时域图和时频图
    输入：
    t: 0-1时间序列
    imfs: IMFs成分
    N: 1s内采样点数
    n: 指定画前几个IMFs成分
    输出：
    前n个IMFs的时域图和时频图
    '''
    fig2, axes = plt.subplots(n, 2, figsize=(10, 7), sharex='col', sharey=False)
    # 计算并绘制各个组分
    for iter in range(n):
        # 绘制分解后的IMF时域图
        axes[iter][0].plot(t, imfs[iter])
        axes[iter][0].set_xlabel('时间/s')
        axes[iter][0].set_ylabel('幅值')
        # 计算各组分的Hilbert变换
        imfsHT = hilbert(imfs[iter])
        # 计算各组分Hilbert变换后的瞬时频率
        instf, timestamps = tftb.processing.inst_freq(imfsHT)
        # 绘制瞬时频率，这里乘以fs是正则化频率到真实频率的转换
        fs = N
        axes[iter][1].plot(timestamps/fs, instf * fs)
        axes[iter][1].set_xlabel('时间/s')
        axes[iter][1].set_ylabel('频率/Hz')
        # 计算瞬时频率的均值和中位数
        axes[iter][1].set_title('Freq_Mean{:.2f}----Freq_Median{:.2f}'.format(np.mean(instf * fs), np.median(instf * fs)))
        
def HHTFilter(signal, componentsRetain):
    '''
    定义HHT的滤波函数，提取部分EMD组分
    输入：
    signol: 信号
    componentsRetain: IMF的列表 []
    输出：
    一幅图，同时包含原始信号和合成信号
    '''
    # 进行EMD分解
    decomposer = EMD(signal)
    # 获取EMD分解后的IMF成分
    imfs = decomposer.decompose()
    # 选取需要保留的EMD组分，并且将其合成信号
    signalRetain = np.sum(imfs[componentsRetain], axis=0)
    # 绘图
    plt.figure(figsize=(10, 7))
    # 绘制原始数据
    plt.plot(signal, label='RawData')
    # 绘制保留组分合成的数据
    plt.plot(signalRetain, label='HHTData')
    # 绘制标题
    plt.title('RawData-----HHTData')
    # 绘制图例
    plt.legend()
    plt.show()
    return signalRetain


# 生成0-1时间序列，共2048个点
N = 2048
t = np.linspace(0, 1, N)
# 生成信号
signal = (2+np.cos(8*np.pi*t)) * np.cos(40*np.pi*(t+1)**2) + np.cos(20*np.pi*t + 5*np.sin(2*np.pi*t))
# 画出原始信号的时域图和频谱
picture(t, signal, N)
# 进行HHT分析，画出所有IMFs的时域图和频谱
imfs = HHTAnalysis(t, signal, N)
# 画出前2个IMFs的时域图和时频图
HHTPicture(t, imfs, N, 2)
# 进行验证，判断与原始信号的差异
signalRetain = HHTFilter(signal, [0, 1])

第一次发文章，如果对你有用，可以请您点个赞吗