#?%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的計算分析函數(shù)
def?HHTAnalysis（eegRaw?fs）:
????#?進行EMD分解
????decomposer?=?EMD（eegRaw）
????#?獲取EMD分解后的IMF成分
????imfs?=?decomposer.decompose（）
????#?分解后的組分數(shù)
????n_components?=?imfs.shape[0]
????#?定義繪圖，包括原始數(shù)據(jù)以及各組分數(shù)據(jù)
????fig?axes?=?plt.subplots（n_components?+?1?2?figsize=（10?7）?sharex=True?sharey=False）
????#?繪制原始數(shù)據(jù)
????axes[0][0].plot（eegRaw）
????#?原始數(shù)據(jù)的Hilbert變換
????eegRawHT?=?hilbert（eegRaw）
????#?繪制原始數(shù)據(jù)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）
????#?計算瞬時頻率的均值和中位數(shù)
????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）
????????#?計算瞬時頻率的均值和中位數(shù)
????????axes[iter?+?1][1].set_title（‘Freq_Mean{:.2f}----Freq_Median{:.2f}‘.format（np.mean（instf?*?fs）?np.median（instf?*?fs）））
????plt.show（）


#?定義HHT的濾波函數(shù)，提取部分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））
????#?繪制原始數(shù)據(jù)
????plt.plot（eegRaw?label=‘RawData‘）
????#?繪制保留組分合成的數(shù)據(jù)
????plt.plot（eegRetain?label=‘HHTData‘）
????#?繪制標題
????plt.title（‘RawData-----HHTData‘）
????#?繪制圖例
????plt.legend（）
????plt.show（）
????return?eegRetain

if?__name__?==?‘__main__‘:
####################示例數(shù)據(jù)分析##################################
????#?生成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]）
######################真實數(shù)據(jù)分析#################################
????#?加載fif格式的數(shù)據(jù)
????epochs?=?mne.read_epochs（r‘F:\BaiduNetdiskDownload\BCICompetition\BCICIV_2a_gdf\Train\Fif\A02T_epo.fif‘）
????#?獲取采樣頻率
????sfreq?=?epochs.info[‘sfreq‘]
????#?想要分析的數(shù)據(jù)
????eegData?=?epochs

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????目錄???????????0??2020-09-07?18:13??Python中HHT（希爾伯特-黃變換）以及其在EEG數(shù)據(jù)處理中的應用\
?????文件????29231312??2020-05-09?16:39??Python中HHT（希爾伯特-黃變換）以及其在EEG數(shù)據(jù)處理中的應用\A02T_epo.fif
?????文件????????3907??2020-09-07?17:49??Python中HHT（希爾伯特-黃變換）以及其在EEG數(shù)據(jù)處理中的應用\HHT_Analysize.py
?????文件?????1238670??2020-09-03?09:36??Python中HHT（希爾伯特-黃變換）以及其在EEG數(shù)據(jù)處理中的應用\pyhht-readthedocs-io-en-latest.pdf
?????文件?????3897122??2020-09-02?18:31??Python中HHT（希爾伯特-黃變換）以及其在EEG數(shù)據(jù)處理中的應用\單通道腦電信號眼電偽跡去除算法研究.pdf