运动想象脑电特征分析：ERSP 值计算与保存（Python 版）

在脑电预处理完成后，接下来进入运动想象脑电特征分析的关键步骤——事件相关谱扰动（ERSP）的计算与保存。将 ERSP 值的计算独立出来，是因为该过程耗时较长，而绘图阶段往往需要反复调整参数。提前保存计算结果，可以直接调用，避免重复运算。

ERSP 值原理简介

ERSP 概念最早由 Makeig 于 1993 年提出，旨在解决听觉实验中事件相关电位（ERP）无法被稳定捕获的问题。ERP 是锁时锁相的活动，而 ERSP 则是锁时但不锁相的节律性活动。ERSP 测量相对于实验事件的宽带脑电频谱的平均振幅随时间的动态变化。

单个试次的事件相关谱估计公式为： $$ P_i = \frac{|F_i(f, t)|^2}{\mu_B(f, i)} $$ 其中 $F_i(f, t)$ 是第 i 个试次的频率 f 和时间 t 的谱估计。通常使用快速傅里叶变换（FFT），配合汉宁窗进行频谱估计。$\mu_B(f, i)$ 代表基线期的平均频谱估计。

考虑到 EEG 信号的非平稳性和高变异性，直接观察单次试次难以捕捉规律。因此，我们通过对多个试次的 ERSP 值取平均，获得更平稳的特征。平均 ERSP 计算公式如下： $$ ERSP_{mean} = 10 \log_{10}(S(f, t)) $$ 其中 $S(f, t)$ 是所有试次在频率 f 处的平均谱估计除以基线平均值。

数据准备与滤波

首先加载预处理后的数据 data，维度通常为（导联数，时间点数，试次数）。为了便于后续计算，需要将维度转换为（时间点数，导联数，试次数）。虽然预处理阶段可能已经做过滤波，但为了保证一致性，这里再次应用带通滤波器（例如 1-40Hz）。

import os
import numpy as np
from pathlib import Path
import hdf5storage
from eeg_filter import ERPs_Filter
from eeg_stft import m_tfr
from scipy.signal import butter, cheby1, filtfilt, resample

# 应用滤波函数
Data = ERPs_Filter(data, freqs=freqwindow, fs=fs, filterflag='filtfilt')
# 转置维度：(time, channels, trials)
X1 = np.transpose(Data, (1, 0, 2))

这里的 ERPs_Filter 是一个封装好的函数，负责通道选择、时间窗划分和滤波操作。其核心逻辑包括去趋势处理以及 Butterworth 滤波器的应用。

# eeg_filter.py
import numpy as np
from scipy.signal import butter, filtfilt, lfilter, detrend

def ERPs_Filter(data, freqs, channel=None, timewindow=, fs=, filterorder=, filterflag=):
    
    X = np.array(data, dtype=, copy=)
    
    
     channel     (channel) > :
         (channel) > X.shape[]:
            ()
        :
            X = X[np.array(channel)-, :, :]
            
    
     timewindow     (timewindow) > :
        t0, t1 = (timewindow[], timewindow[-])  (timewindow) >   (, timewindow[])
        t1 = (t1, X.shape[]/fs)
        s0 = ((t0*fs))
        s1 = ((t1*fs))
        X = X[:, s0:s1, :]
        
    
    low, high = freqs
    Wn = [*low/fs, *high/fs]
    b, a = butter(filterorder, Wn, btype=)
    Y = np.zeros_like(X)
    
     s  (X.shape[]):
        data1 = detrend(X[:, :, s].T, axis=)
         filterflag == :
            data2 = filtfilt(b, a, data1, axis=)
        :
            data2 = lfilter(b, a, data1, axis=)
        Y[:, :, s] = data2.T
     Y

# eeg_stft.py 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})

运动想象脑电特征分析：ERSP 值计算与保存（Python 版）