Python中HHT(希尔伯特-黄变换)以及其在EEG数据处理中的应用

最新推荐文章于 2024-03-07 10:43:56 发布

zhoudapeng01

最新推荐文章于 2024-03-07 10:43:56 发布

阅读量1.1w

点赞数 22

分类专栏： MNE 算法 python 文章标签： python 算法 HHT

本文链接：https://blog.csdn.net/zhoudapeng01/article/details/108448218

版权

python 同时被 3 个专栏收录

20 篇文章 8 订阅

订阅专栏

算法

19 篇文章 7 订阅

订阅专栏

MNE

15 篇文章 52 订阅

订阅专栏

在脑电信号的处理过程中去除伪迹是很关键的一个步骤，常用的有ICA和小波等方法。不过这些方法大多是针对多通道脑电数据进行处理的，单通道的脑电数据如何去除伪迹呢？推荐一篇文章《单通道脑电信号眼电伪迹去除算法研究》，在文章中提到了一种WT-EEMD-ICA方法，该方法是小波-集合经验模态分解-独立成分分析的结合。具体内容感兴趣的可以精读下这篇文章，在对应的下载附件中有这篇文章。

博客相关资源下载：https://download.csdn.net/download/zhoudapeng01/12820983

上面说的和本篇的内容关系不大，我就是在看了文章后对里面提到的HHT方法感兴趣，就研究了一下。下面主要说的是HHT的实现以及如何准确计算瞬时频率。

https://www.cnblogs.com/xingshansi/p/6498913.html

https://www.cnblogs.com/RoseVorchid/p/11985368.html

https://www.pythonf.cn/read/76270

相关代码：

# %matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from pyhht import EMD
from scipy.signal import hilbert
import tftb.processing
import mne


# 定义HHT的计算分析函数
def HHTAnalysis(eegRaw, fs):
    # 进行EMD分解
    decomposer = EMD(eegRaw)
    # 获取EMD分解后的IMF成分
    imfs = decomposer.decompose()
    # 分解后的组分数
    n_components = imfs.shape[0]
    # 定义绘图，包括原始数据以及各组分数据
    fig, axes = plt.subplots(n_components + 1, 2, figsize=(10, 7), sharex=True, sharey=False)
    # 绘制原始数据
    axes[0][0].plot(eegRaw)
    # 原始数据的Hilbert变换
    eegRawHT = hilbert(eegRaw)
    # 绘制原始数据Hilbert变换的结果
    axes[0][0].plot(abs(eegRawHT))
    # 设置绘图标题
    axes[0][0].set_title('Raw Data')
    # 计算Hilbert变换后的瞬时频率
    instf, timestamps = tftb.processing.inst_freq(eegRawHT)
    # 绘制瞬时频率，这里乘以fs是正则化频率到真实频率的转换
    axes[0][1].plot(timestamps, instf * fs)
    # 计算瞬时频率的均值和中位数
    axes[0][1].set_title('Freq_Mean{:.2f}----Freq_Median{:.2f}'.format(np.mean(instf * fs), np.median(instf * fs)))

    # 计算并绘制各个组分
    for iter in range(n_components):
        # 绘制分解后的IMF组分
        axes[iter + 1][0].plot(imfs[iter])
        # 计算各组分的Hilbert变换
        imfsHT = hilbert(imfs[iter])
        # 绘制各组分的Hilber变换
        axes[iter + 1][0].plot(abs(imfsHT))
        # 设置图名
        axes[iter + 1][0].set_title('IMF{}'.format(iter))
        # 计算各组分Hilbert变换后的瞬时频率
        instf, timestamps = tftb.processing.inst_freq(imfsHT)
        # 绘制瞬时频率，这里乘以fs是正则化频率到真实频率的转换
        axes[iter + 1][1].plot(timestamps, instf * fs)
        # 计算瞬时频率的均值和中位数
        axes[iter + 1][1].set_title('Freq_Mean{:.2f}----Freq_Median{:.2f}'.format(np.mean(instf * fs), np.median(instf * fs)))
    plt.show()


# 定义HHT的滤波函数，提取部分EMD组分
def HHTFilter(eegRaw, componentsRetain):
    # 进行EMD分解
    decomposer = EMD(eegRaw)
    # 获取EMD分解后的IMF成分
    imfs = decomposer.decompose()
    # 选取需要保留的EMD组分，并且将其合成信号
    eegRetain = np.sum(imfs[componentsRetain], axis=0)

    # 绘图
    plt.figure(figsize=(10, 7))
    # 绘制原始数据
    plt.plot(eegRaw, label='RawData')
    # 绘制保留组分合成的数据
    plt.plot(eegRetain, label='HHTData')
    # 绘制标题
    plt.title('RawData-----HHTData')
    # 绘制图例
    plt.legend()
    plt.show()
    return eegRetain

if __name__ == '__main__':
####################示例数据分析##################################
    # 生成0-1时间序列，共100个点
    t = np.linspace(0, 1, 100)
    # 生成频率为5Hz、10Hz、25Hz的正弦信号累加
    modes = np.sin(2 * np.pi * 5 * t) + np.sin(2 * np.pi * 10 * t) + np.sin(2 * np.pi * 25 * t)
    # 信号和时间累加，相当于添加噪声
    eegRaw = modes + t
    # fs为采样频率，用于正则化频率和真实频率的转换
    fs = 100
    # 进行HHT分析
    HHTAnalysis(eegRaw, fs)
    # 选取需要保留的信号进行合成，也就是相当于滤波
    eegRetain = HHTFilter(eegRaw, [0, 1, 2, 3])
######################真实数据分析#################################
    # 加载fif格式的数据
    epochs = mne.read_epochs(r'F:\BaiduNetdiskDownload\BCICompetition\BCICIV_2a_gdf\Train\Fif\A02T_epo.fif')
    # 获取采样频率
    sfreq = epochs.info['sfreq']
    # 想要分析的数据
    eegData = epochs.get_data()[0][0]
    HHTAnalysis(eegData, sfreq)
    HHTFilter(eegData, [0, 1, 2, 3, 4, 5, 6, 7, 8])