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

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

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

推荐几个参考的博客：

相关代码：

# %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])

一、EMD分解

HHT变换首先要做的是EMD也就是经验模态分解，原始信号就是分解后信号的累加和。开始的时候也是有些疑问的，后面尝试了下，发现真的是这样的啊，这可比小波、ICA这种分解好理解的多了，简单一句话就是将信号分解为不同的成分。

原始信号和分解后合成信号的对比：分解后再合成和原始信号基本保持一致。

开始我还以为计算错了，去掉了IMF0，再看，确实有效果：

二、瞬时频率

查了些资料，还是没搞明白瞬时频率到底是什么意思，在这个场景中我的理解是，其目的就是为了计算EMD分解后信号的频率。不过在一些需要时频分析的场景中应该更加有用吧。

先看一下5Hz+10Hz+25Hz+t信号的情况：EMD分解了4个组分IMF0+IMF1+IMF2+IMF3，IMF0\IMF1\IMF2其对应瞬时频率的中位数为25.24、10.18、5.01，和生成信号的频率相对应，还是比较准确的，取中位数主要是为了去掉边界的影响。

这里有一点需要注意：tftb.processing.inst_freq计算的瞬时频率是正则化后的，如果想要得到真实频率需要进行换算。Matlab中的单位频率指的是奈圭斯特频率(采样频率的一半)，在Python中试验了下这里应该是采样频率，而非其一半。

三、应用到真实数据中

HHT这种方法有一个好处，需要确定的参数比较少。和STFT相比不用设置窗长，和WP相比不用设置基函数。选取了一段竞赛数据处理看着还像那么回事，这次就处理到这，后面有机会再进行后续的相关工作吧。

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