《python希尔伯特变换-Python中HHT（希尔伯特-黄变换）以及其在EEG数据处理中》由会员分享，可在线阅读，更多相关《python希尔伯特变换-Python中HHT（希尔伯特-黄变换）以及其在EEG数据处理中（4页珍藏版）》请在金锄头文库上搜索。
1、python希尔伯特变换_Python中HHT（希尔伯特-黄变换）以及其在EEG数据处理中。在脑电信号的处理过程中去除伪迹是很关键的个步骤，常的有ICA和波等法。不过这些法多是针对多通道脑电数据进处理的，单通道的脑电数据如何去除伪迹呢？推荐篇章单通道脑电信号眼电伪迹去除算法研究，在章中提到了种WT-EEMD-ICA法，该法是波-集合经验模态分解-独成分分析的结合。具体内容感兴趣的可以精读下这篇章，在对应的下载附件中有这篇章。上说的和本篇的内容关系不，我就是在看了章后对提到的HHT法感兴趣，就研究了下。下主要说的是HHT的实现以及如何准确计算瞬时频率。推荐个参考的博客：相关代码：# %matpl
2、otlib inlineimport matplotlib.pyplot as pltimport numpy as npfrom pyhht import EMDfrom scipy.signal import hilbertimport tftb.processingimport mne# 定义HHT的计算分析函数def HHTAnalysis(eegRaw, fs):# 进EMD分解decomposer = EMD(eegRaw)# 获取EMD分解后的IMF成分imfs = decomposer.decompose()# 分解后的组分数n_components = imfs.shape0
3、# 定义绘图，包括原始数据以及各组分数据fig, axes = plt.subplots(n_components + 1, 2, figsize=(10, 7), sharex=True, sharey=False)# 绘制原始数据axes00.plot(eegRaw)# 原始数据的Hilbert变换eegRawHT = hilbert(eegRaw)# 绘制原始数据Hilbert变换的结果axes00.plot(abs(eegRawHT)# 设置绘图标题axes00.set_title(Raw Data)# 计算Hilbert变换后的瞬时频率instf, timestamps = tftb
4、.processing.inst_freq(eegRawHT)# 绘制瞬时频率，这乘以fs是正则化频率到真实频率的转换axes01.plot(timestamps, instf * fs)# 计算瞬时频率的均值和中位数axes01.set_title(Freq_Mean:.2f-Freq_Median:.2f.format(np.mean(instf * fs), np.median(instf * fs)# 计算并绘制各个组分for iter in range(n_components):# 绘制分解后的IMF组分axesiter + 10.plot(imfsiter)# 计算各组分的Hil
5、bert变换imfsHT = hilbert(imfsiter)# 绘制各组分的Hilber变换axesiter + 10.plot(abs(imfsHT)# 设置图名axesiter + 10.set_title(IMF.format(iter)# 计算各组分Hilbert变换后的瞬时频率instf, timestamps = tftb.processing.inst_freq(imfsHT)# 绘制瞬时频率，这乘以fs是正则化频率到真实频率的转换axesiter + 11.plot(timestamps, instf * fs)# 计算瞬时频率的均值和中位数axesiter + 11.se
6、t_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(imfscomponentsRetain,
7、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 eegRetainif _name_ = _main_:#例数据分析# 成0-1时间序列，共100个点t = np.linspace(0, 1, 100)# 成频率为5Hz、10Hz、25Hz的正弦信号累加modes =
8、 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(rF:BaiduNe
9、tdiskDownloadBCICompetitionBCICIV_2a_gdfTrainFifA02T_epo.fif)# 获取采样频率sfreq = epochs.infosfreq# 想要分析的数据eegData = epochs.get_data()00HHTAnalysis(eegData, sfreq)HHTFilter(eegData, 0, 1, 2, 3, 4, 5, 6, 7, 8)、EMD分解HHT变换先要做的是EMD也就是经验模态分解，原始信号就是分解后信号的累加和。开始的时候也是有些疑问的，后尝试了下，发现真的是这样的啊，这可波、ICA这种分解好理解的多了，简单句话就
10、是将信号分解为不同的成分。原始信号和分解后合成信号的对：分解后再合成和原始信号基本保持致。开始我还以为计算错了，去掉了IMF0，再看，确实有效果：、瞬时频率查了些资料，还是没搞明瞬时频率到底是什么意思，在这个场景中我的理解是，其的就是为了计算EMD分解后信号的频率。不过在些需要时频分析的场景中应该更加有吧。先看下5Hz+10Hz+25Hz+t信号的情况：EMD分解了4个组分IMF0+IMF1+IMF2+IMF3，IMF0IMF1IMF2其对应瞬时频率的中位数为25.24、10.18、5.01，和成信号的频率相对应，还是较准确的，取中位数主要是为了去掉边界的影响。这有点需要注意：tftb.processing.inst_freq计算的瞬时频率是正则化后的，如果想要得到真实频率需要进换算。Matlab中的单位频率指的是奈圭斯特频率(采样频率的半)，在Python中试验了下这应该是采样频率，其半。三、应到真实数据中HHT这种法有个好处，需要确定的参数较少。和STFT相不设置窗长，和WP相不设置基函数。选取了段竞赛数据处理看着还像那么回事，这次就处理到这，后有机会再进后续的相关作吧。
