# 希爾伯特黃轉換(Hilbert Huang Transform)
這邊需要用到 `PyEMD` 套件，以及翻譯自黃院士的 matlab libary 的 `nspplota`
首先來設定一下執行整套 HHT 需要的計算資源
```python
import mne
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from mne.time_frequency import tfr_stockwell
# HHT tools
from PyEMD import EEMD as py_eemd
from scipy.signal import hilbert, medfilt
#from eemd_faphi import *
from nspplota import nspplota
from scipy.ndimage import gaussian_filter, uniform_filter
EEMD_class = py_eemd(trials = 40, noise_width = 0.4, max_imf = 6)
EEMD_class.EMD.FIXE = 10 # number of sifting
# temp_: 2D MEG/EEG data, n_trial x times
# ntrl: integer, number of trial
# nimf: integer, number of IMFs
# ntime: integer, number of time frames
def DoEEMD(temp_, ntrl, nimf, ntime, EEMD_class):
print('Do EEMD....')
data_imfs = np.zeros((ntime, nimf + 1, ntrl))
for trl in range(ntrl):
imf_ = EEMD_class(temp_[trl,:], max_imf = nimf).T
data_imfs[:,:,trl] = imf_
return data_imfs
def calc_inst_fa_hilbert(sig, dt):
""""Extracts instantaneous frequency through the Hilbert Transform."""
analytic_signal = hilbert(sig) # Apply Hilbert transform to each row
inst_a = np.abs(analytic_signal)
inst_phase = np.unwrap(np.angle(analytic_signal)) # Compute angle between img and real
inst_f1 = np.diff(inst_phase)/(2*np.pi*dt)
inst_f = np.hstack((inst_f1, np.tile(inst_f1[:, [-1]], 1))) # Duplicate last points to make dimensions equal
# medien filter
for n in range(inst_f.shape[0]):
inst_f[n] = medfilt(inst_f[n], 5)
pass
inst_a = np.abs(analytic_signal)
return inst_f, inst_a, inst_phase
data_path = '/Users/kevinhsu/Documents/D/000_experiment/project_cospro/PG_meg_data/001'
sid = 1
epochs_file = data_path + '/%.3d_raw_epo.fif' % sid
epochs = mne.read_epochs(epochs_file)
epochs
```
------------
## Part 1.
這裡面用到 `mne` 的 `tfr_stockwell`，純粹是製作一個儲存 time-frequency 資料的容器：
```python
freqs = np.arange(3., 60., 2.)
fmin, fmax = freqs[[0, -1]]
width = 0.9
power1kHz = tfr_stockwell(epochs['Tone1kHz'], fmin=fmin, fmax=fmax, width=width)
powerAM40hz = tfr_stockwell(epochs['AM40Hz'], fmin=fmin, fmax=fmax, width=width)
print(type(power1kHz))
```
<class 'mne.time_frequency.tfr.AverageTFR'>
這邊建立了兩個 `AverageTFR` class 的物件，`power1kHz` 以及`powerAM40hz`。和 `raw`, `epochs`, `evoked` 這三個 class 一樣，`AverageTFR` 底下有一個 numpy ndarray 的資料。
以下的步驟是先做基準線校正，以事件發生之前 (-0.1~0s) 的區間的 power 平均值當作分母，去對比每一個時間點的 power, 取得相對值之後，在做 log 轉換。
最後結果如果是正值，表示事件發生之後的訊號大於baseline的訊號 (event-related synchronization ERS)
如果是負值，表示事件發生之後的訊號小於baseline的訊號 (event-related desynchronization ERD)
```python
power1kHz.apply_baseline((-0.1, 0), mode = 'logratio')
powerAM40hz.apply_baseline((-0.1, 0), mode = 'logratio')
```
然後從兩個物件裡面取出資料並且相減，以凸顯 40Hz 活動在兩組情境之下的差異：
```python
powerDiff_amp = powerAM40hz.data - power1kHz.data
```
這時候 `powerDiff_amp` 只是兩個 numpy ndarray 相減的結果，本身也是 numpy ndarray, 接著把`powerDiff_amp` 變成 `AverageTFR` class
實際上的做法是，利用`mne.time_frequency.AverageTFR`這個指令，產生一個 `AverageTFR` class 的物件，並且把 `powerDiff_amp` 放到這個新的物件裡面。根據網頁說明，一個 `AverageTFR` 至少需要下面這些資訊：
?mne.time_frequency.AverageTFR
```
info : The measurement info.
data : ndarray, shape (n_channels, n_freqs, n_times)
--The data.
times : ndarray, shape (n_times,)
--The time values in seconds.
freqs : ndarray, shape (n_freqs,)
--The frequencies in Hz.
nave : int
--The number of averaged TFRs.
```
基本上除了第二項要用`powerDiff_amp`，其他四項可以直接繼承自 `power1kHz` 或 `powerAM40hz`
```python
powerDiff = mne.time_frequency.AverageTFR(powerAM40hz.info,
powerDiff_amp, powerAM40hz.times,
powerAM40hz.freqs, powerAM40hz.nave)
```
然後就能繪圖：
```python
powerDiff.plot_topo()
```
------------
## Part 2.
到此為止，沒有用到任何 HHT 的概念。我們接下來要做的是：
1. 從`epochs` 取出 data array，維度是 (193, 157, 700), events x channels x times
2. 針對每一段訊號 (event by event, channel by channel)逐一執行 HHT 的第一步驟，EEMD。假設每一個訊號分解成 7 個 IMFs，資料的維度變成 (193, 157, 700, 7)
3. 針對每一組 IMFs (總共 193 x 157 x 700) 組，取出前面 6 項 IMF, 然後計算瞬時頻率、瞬時震幅。這兩個資料的大小和 IMFs 的量一樣，都是 (193, 157, 700, 7)
4. 每一組 IMFs 會有一對 瞬時頻率、瞬時震幅。將這兩項資訊投射到 time-frequency 空間。頻率至個維度取 50 個 bins。然後將 events 這個維度平均再一起。最後的結果是 (157, 50, 700) 的三維矩陣。
### 2.1. **EEMD**
第一個步驟要用的就是 `epochs.get_data()`，這邊先只挑一種事件來作HHT。這邊要注意一個重點，mne裡面的資料是用沒有轉換小數點的格式，MEG是T (tesla), EEG是V (voltage)
由於 MEG sensor 的解析度是到fT，所以一開始用`get_data()`取得的數值會很小，取出來的資料需要轉換成fT
```python
MEG_data = epochs['AM40Hz'].get_data() * 1e+15
MEG_data = MEG_data.transpose((1,0,2))
MEG_data.shape
```
(157, 95, 700)
為了快速示範，這邊只挑一個 channel (MEG 103)
```python
MEG_data_select = MEG_data[102, :, :] # T >> fT
MEG_data_select.shape
```
(95, 700)
現在取得95段event的訊號，每一段的長度是700個取樣點。以下針對這95段訊號執行EEMD
```python
ntrl, ntimes = MEG_data_select.shape # ntrl: number of trials; ntimes: number of time frames
nimf = 6 # number of IMFs
#data_imfs = DoEEMD( MEG_data_select, 193, 6, 700, EEMD_class)
data_imfs = DoEEMD( MEG_data_select, ntrl, nimf, ntimes, EEMD_class)
```
這邊要注意一下，`nimf` 是設定IMF的數量。PyEMD會再額外加上一個IMF不能解釋的殘差當作做後一項資料，所以做後輸出的IMF的維度是 `nimf + 1`：
```python
data_imfs.shape
```
(700, 7, 95) >> times x imf x events
### 2.2. **Hilbert transform**
接著，使用 IMF 計算瞬時頻率、瞬時震幅。這邊先放棄瞬時相位...Orz
等之後我補上院士在2002年使用的 `direct quadrature` 再來玩 phase
```python
dt = 1.0 / epochs.info['sfreq']
# data_imfs[:,0:-1,0].T.shape >> (6, 700)
inst_f, inst_a, inst_phase = calc_inst_fa_hilbert(data_imfs[:,0:-1,0].T, dt)
```
可以用這個指令看一下瞬時震幅與IMF之間的關聯：
```python
fig, axes = plt.subplots(6,1)
temp_imf = data_imfs[:,0:-1,0]
for num in range(6):
ax = axes[num]
ax.plot(epochs.times, temp_imf[:, num], 'b')
ax.plot(epochs.times, inst_a[num, :], 'r')
#ax.plot(epochs.times, inst_f[num, :], 'b')
ax.set_ylabel("IMF " + str(num+1))
```
接下來就是有一點噁心的步驟了。上面的 `calc_inst_fa_hilbert` 只有抓一段 event 訊號的 6 個 IMFs。下面要做的是，把95段event各自的6個IMF取出來執行`calc_inst_fa_hilbert`...
```python
inst_f_full = np.zeros((data_imfs.shape[2], data_imfs.shape[1]-1, data_imfs.shape[0]))
inst_a_full = np.zeros((data_imfs.shape[2], data_imfs.shape[1]-1, data_imfs.shape[0]))
for num in range(data_imfs.shape[2]):
print(num)
inst_f, inst_a, inst_phase = calc_inst_fa_hilbert(data_imfs[:,0:-1,num].T, dt)
inst_f_full[num] = inst_f
inst_a_full[num] = inst_a
pass
```
這邊的`inst_f_full` 和 `inst_a_full` 都是 events x imf x times 的矩陣。把event 放在第一個維度主要是方便後續的分析步驟
### 2.1. **make a beautiful picture**
到此為止，前面說的1-4個步驟裡，我們用一個 channel 完成了1~3
接下來要做第四步。先用一段MEG訊號的instantaneous frequency、instantaneous amplitude來示範怎麼使用`nspplota`:
```python
# 定義nspplota的參數
t0 = epochs.times[0] # 時間序列訊號的起始時間，等於epochs物件底下的 times 的第一個數值 (s)
t1 = epochs.times[-1] # 時間序列訊號的起始時間，用epochs物件底下的 times 的最後一個數值 (s)
fres = 50 # 頻率維度的 bin 數量
tres = 100 # 時間維度的 bin 數量
fw0 = 0 # 頻率維度的最小值 (Hz)
fw1 = 100 # 頻率維度的最大值 (Hz)
tw0 = t0 # 同 t0
tw1 = t1 # 同 t1
# 用第一段event的instantaneous frequency、instantaneous amplitude。請注意輸入變數的順序要放對
nt, fscale, tscale, nt_imf_ms = nspplota(inst_f_full[0], inst_a_full[0], dt, t0, t1, fres, tres, fw0, fw1, tw0, tw1)
# nt 的大小等於 fres x tres, 在這邊的結果是 50 x 100
# 第一個維度對應 frequency, 第二個維度則是 times
# 依照院士在 1998 年的報告中，HHT spectra 需要做 兩次 2D
# smoothing 才能凸現 time-frequency 的特徵
nt2 = gaussian_filter(nt, sigma = 0.6)
nt2 = uniform_filter(nt2, size = 5)
nt2 = np.sqrt(nt2)
# 左圖是沒有 smoothing 的 HHT spectrum, 右圖是
fig, axes = plt.subplots(1,2)
ax = axes[0]
ax.imshow(nt, extent=[tscale[0],tscale[-1],fscale[0],fscale[-1]],
interpolation='none', aspect='auto', origin = 'lower')
ax = axes[1]
ax.imshow(nt2, extent=[tscale[0],tscale[-1],fscale[0],fscale[-1]],
interpolation='none', aspect='auto', origin = 'lower')
```
好了，那再來跑一遍令人頭暈的分析，把同一個MEG channel的95段event 的HHT spectra計算出來，然後平均在一起，得到一個平均之後的 spectra。
```python
t0 = epochs.times[0]
t1 = epochs.times[-1]
fres = 50
tres = 100
fw0 = 0
fw1 = 100
tw0 = t0
tw1 = t1
nt_full = np.zeros((inst_f_full.shape[0], fres, tres))
nt_imf_ms_full = np.zeros((inst_f_full.shape[0], nimf, fres))
for num in range(inst_f_full.shape[0]):
print(num)
nt, fscale, tscale, nt_imf_ms = nspplota(inst_f_full[num], inst_a_full[num], dt, t0, t1, fres, tres, fw0, fw1, tw0, tw1)
nt2 = gaussian_filter(nt, sigma = 0.6)
nt2 = uniform_filter(nt2, size = 5)
nt2 = np.sqrt(nt2)
nt_full[num] = np.nan_to_num(nt2)
nt_imf_ms_full[num] = nt_imf_ms
pass
nt_mean = np.mean(nt_full, axis = 0)
plt.imshow(nt_mean, extent=[tscale[0],tscale[-1],fscale[0],fscale[-1]],
interpolation='none', aspect='auto', origin = 'lower', )
```
到這邊為止，我們把 `MEG 103` 號 channel 上的資料做了下面的處理：
1. 95 段event 的訊號 (95x700), 轉換成 IMF (700x7x95)
2. 每一段event 的 IMF (700x7) 變成 instantaneous frequency、instantaneous amplitude。所以這個步驟產生兩個同樣是 95x7x700 的矩陣
3. 每一段event 的instantaneous frequency、instantaneous amplitude投射到由 fres (50) 和 tres (100) 定義的 HHT spectra空間上。這邊會形成一個 95x50x100 的矩陣
4. 最後，把前一個步驟的矩陣沿著第一個維度 (event) 平均在一起，產生一個 50x100 的 time-frequency spectra
恭喜你！完成人生中的第一個 event-related HHT spectra!!!!
那麼....
# My wish list:
1. 我們總共有 157 個 channels，請把上面的步驟變成一個大大的 loop，最後產生一個 157 x freq x tres 的 numpy ndarray
3. 然後請用 part 1 的方式，把這些 time-frequency spectra 放入 `AverageTFR` class 裡面 (人生中的第一個 全頭式 MEG 的 ER-HHT spectra!!!!)。
放入 `AverageTFR` class 裡面的步驟大致上會像這樣：
```python
dummy_power_ndarray = np.ones((157, fres, tres))
dummy_power_ndarray[102] = nt_mean
dummy_power = mne.time_frequency.AverageTFR(epochs.info, dummy_power_ndarray, tscale, fscale, data_imfs.shape[2])
```
Sign in with Wallet
Connect another wallet