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