srate = 250
freqwindow = [1, 40]
TaskDuration = 4
RestDuration = 2
tlimits = [-1000*RestDuration, TaskDuration*1000 - 1]
for class_choose in [1, 2, 3, 4]:
X1_temp = X1[:, :, label == class_choose]
TFR1, t_idx, f_hz = m_tfr(X1_temp, fs=fs, N=256, win_name='gauss',
win_len=251, out_freq_bins=81, out_time_len=frames)
power1 = 10 * np.log10(np.maximum(np.abs(TFR1) ** 2, np.finfo(float).eps))
times = (np.arange(frames) / fs * 1000.0 + tlimits[0])
freqs = f_hz
bmask = (times >= -1000 * RestDuration) & (times <= 0)
base1 = power1[:, bmask, :, :].mean(axis=1, keepdims=True)
ERSP_1 = power1 - base1
ERSP_1 = ERSP_1.astype(np.float32)
hdf5storage.savemat(os.path.join(save_root_ERSP, f"sub{s}_class{class_choose}.mat"), {'ERSP_1': [ERSP_1]})
hdf5storage.savemat(os.path.join(save_root_ERSP, f"times+freqs_class{class_choose}.mat"),{'times': times, 'freqs': freqs})
from numpy.fft import fft
import numpy as np
def tftb_window(L, win_name='hamming'):
"""生成各种窗函数"""
L = int(L)
if L <= 0: raise ValueError("Window length must be positive")
if L % 2 == 0: L = L + 1
win_name = win_name.lower()
if win_name in ['rect', 'rectangular']: w = np.ones(L)
elif win_name in ['hamming', 'hamming']:
n = np.arange(L); w = 0.54 - 0.46 * np.cos(2 * np.pi * n / (L - 1))
elif win_name in ['hanning', 'hann']:
n = np.arange(L); w = 0.5 * (1 - np.cos(2 * np.pi * n / (L - 1)))
elif win_name == 'bartlett':
n = np.arange(L); w = 1 - 2 * np.abs(n - (L - 1) / 2) / (L - 1)
elif win_name == 'triang':
n = np.arange(L); w = 1 - 2 * np.abs(n - (L - 1) / 2) / L
elif win_name == 'blackman':
n = np.arange(L)
w = (0.42 - 0.5 * np.cos(2 * np.pi * n / (L - 1)) + 0.08 * np.cos(4 * np.pi * n / (L - 1)))
elif win_name in ['gauss', 'gaussian']:
n = np.arange(L); alpha = 2.5
w = np.exp(-0.5 * (alpha * (n - (L - 1) / 2) / ((L - 1) / 2)) ** 2)
elif win_name == 'flattop':
n = np.arange(L)
a0, a1, a2, a3, a4 = 0.21557895, 0.41663158, 0.277263158, 0.083578947, 0.006947368
w = (a0 - a1 * np.cos(2 * np.pi * n / (L - 1)) + a2 * np.cos(4 * np.pi * n / (L - 1)) - a3 * np.cos(6 * np.pi * n / (L - 1)) + a4 * np.cos(8 * np.pi * n / (L - 1)))
else: raise ValueError(f"Unknown window type: {win_name}")
w = w.astype(np.float64)
w = w / np.linalg.norm(w)
return w
def tfrstft(x, N=None, h=None):
"""等价 MATLAB: [tfr,t,f] = tfrstft(x, t, N, h)"""
x = np.asarray(x).reshape(-1)
T = x.shape[0]
if N is None: N = T
if h is None:
hlen = int(np.floor(N/4.0))
if hlen % 2 == 0: hlen += 1
h = tftb_window(hlen, 'hamming')
h = np.asarray(h).reshape(-1, 1)
if h.shape[0] % 2 == 0: raise ValueError("Window length must be odd")
h = h / np.linalg.norm(h)
Lh = (h.shape[0]-1)//2
tfr = np.zeros((N, T), dtype=complex)
halfN = int(np.floor(N/2.0))
for ti in range(T):
tau_min = -min(halfN-1, Lh, ti)
tau_max = min(halfN-1, Lh, T-1-ti)
tau = np.arange(tau_min, tau_max+1)
idx = ( (N + tau) % N )
tfr[idx, ti] = x[ti+tau] * np.conjugate(h[Lh+tau, 0])
tfr = fft(tfr, axis=0)
if N % 2 == 0:
f = np.concatenate([np.arange(0, N//2), np.arange(-N//2, 0)]) / N
else:
f = np.concatenate([np.arange(0, (N-1)//2 + 1), np.arange(-(N-1)//2, 0)]) / N
return tfr, f
def m_tfr(data, fs, N=256, win_name='gauss', win_len=287, out_freq_bins=53, out_time_len=1500):
"""输入 data 维度可为 (time, channels, trials) 或 (channels, time, trials)"""
if data.shape[0] in (1500, 1800, 2000): time_major = True
else: time_major = False
if time_major:
T, C, Tr = data.shape; X = data
else:
C, T, Tr = data.shape; X = np.transpose(data, (1, 0, 2))
h = tftb_window(win_len, win_name)
ttfr = np.zeros((out_freq_bins, out_time_len, Tr, C), dtype=complex)
for ci in range(C):
for tj in range(Tr):
tfr, f_norm = tfrstft(X[:, ci, tj], N=N, h=h)
tfr_pos = tfr[:out_freq_bins, :out_time_len]
ttfr[:, :, tj, ci] = tfr_pos
df = fs / float(N)
f_hz = np.arange(out_freq_bins) * df
t_idx = np.arange(1, out_time_len+1)
return ttfr, t_idx, f_hz
import os
import hdf5storage
import numpy as np
from pathlib import Path
topo_band = [(8,13),(13, 30),(13,20),(20,30)]
topo_timewin = [(0, 4),(-2,0),(0,0.5),(0.5,1),(1,1.5),(1.5,2),(2,2.5),(2.5,3),(3,3.5),(3.5,4)]
CHANNELS = [ 'FP1','FPZ','FP2','AF3','AF4','F7','F5','F3','F1','FZ','F2','F4','F6','F8',
'FT7','FC5','FC3','FC1','FCZ','FC2','FC4','FC6','FT8','T7','C5','C3','C1',
'CZ','C2','C4','C6','T8','TP7','CP5','CP3','CP1','CPZ','CP2','CP4','CP6',
'TP8','P7','P5','P3','P1','PZ','P2','P4','P6','P8','PO7','PO5','PO3','POZ',
'PO4','PO6','PO8','O1','OZ','O2' ]
save_dir = r'.\ERSP 数据'
save_root_TOPO = r'.\脑地形图数据'
Path(save_root_TOPO).mkdir(parents=True, exist_ok=True)
save_root_TF = r'.\时频图数据'
Path(save_root_TF).mkdir(parents=True, exist_ok=True)
subjects = list(range(1, 37 + 1))
CLASS_CHOOSE_LIST = [1, 2]
Class_name = ('Left', 'Right')
for class_choose in CLASS_CHOOSE_LIST:
topo, TOPO = {'s1': {}, 's2': {}}, {'s1': {}, 's2': {}}
tf, TF = {'s1': [], 's2': []}, {'s1': None, 's2': None}
tf_path = os.path.join(save_dir, f"times+freqs_class{class_choose}.mat")
times, freqs = hdf5storage.loadmat(tf_path)['times'], hdf5storage.loadmat(tf_path)['freqs']
Flen, Tlen = len(times), len(freqs)
for band in topo_band:
for time in topo_timewin:
topo['freq' + str(band) + 'time' + str(time)] = np.empty((len(CHANNELS), 0))
for s in subjects:
p1 = os.path.join(save_dir, f"sub{s}_class{class_choose}.mat")
ersp_s1 = hdf5storage.loadmat(p1)['ERSP_1'][0]
for band in topo_band:
for time in topo_timewin:
topo_temp = ersp_s1[(freqs >= band[0]).__and__(freqs < band[1]), :, :, :]
topo_temp1 = topo_temp[:, (times >= time[0] * 1000).__and__(times < time[1] * 1000), :, :]
topo_mean = np.expand_dims( np.mean(np.mean(np.mean(topo_temp1, axis=0), axis=0), axis=0), axis=1)
topo['freq' + str(band) + 'time' + str(time)] = np.concatenate(
(topo['freq' + str(band) + 'time' + str(time)], topo_mean), axis=1)
tf.append(np.mean(ersp_s1, axis=2))
TF = np.array(tf)
p = os.path.join(save_root_TF, f"TF_class{class_choose}.mat")
hdf5storage.savemat(p, {'tf': TF})
p = os.path.join(save_root_TOPO, f"TOPO_class{class_choose}.mat")
hdf5storage.savemat(p, {'topo': topo})
