一、脑电信号获取
脑电信号是由脑部神经元的兴奋和抑制行为产生的微弱电流信号。为了获取脑电信号，需要使用电极阵列将电极头置于头皮上。脑电信号是一种高噪声、低幅值、高时间分辨率的信号，因此需要进行信号放大和滤波。
二、脑电信号预处理
脑电信号需要进行预处理处理，以消除噪声、伪迹和其他干扰。预处理包括：滤波、空间滤波、伪迹去除。
1. 滤波
def butter_bandpass_filter(data,lowcut,highcut,fs,order=5):
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
y = filtfilt(b, a, data)
return y
使用巴特沃斯滤波器进行带通滤波，过滤掉低频和高频噪音。这个函数的输入数据data是一维的，函数返回的也是一维的数据。
2. 空间滤波
def common_average_reference(data,anchors,bads = []):
'''
common average reference
@data：(n_channels,n_samples)脑电信号数据矩阵
@anchors：(n_channels,n_channels)电极距离矩阵或者十二个电极头的距离矩阵
@bads：(n_bad_chls)默认是空通道数组
'''
n_chls,n_samps = data.shape
all_chls = np.arange(n_chls)
if len(bads)!=0:
all_chls = np.delete(all_chls,bads)
n_active_chls = len(all_chls)
loc = pycsd.estimate_locations(
anchors[all_chls,:][:,all_chls],#使用最小描述长度方法估计电极的三维坐标
ch_xyz=False)
loc = loc[:n_active_chls,:]
tmp_data = data[all_chls,:]
mu = np.mean(tmp_data,axis=0,keepdims=True)
re_tmp_data = tmp_data - mu
csr = pycsd.constant_neighborhood(loc)
csr_noise = pycsd.CSRNoiseThreshold(estimate='median', factor=1.0)
csr.apply(re_tmp_data,out=tmp_data)
csr_noise.apply(tmp_data,out=tmp_data)
csr_back = pycsd.CSRBackProjection(loc)
csr_back.apply(tmp_data,out=tmp_data)
re_data = np.zeros((n_chls,n_samps))
re_data[all_chls,:] = tmp_data + mu
return re_data公共平均参考（car）是一种将整个电极阵列上的平均电压作为参考的预处理方法。空间滤波是一种更高级的预处理方法，它通过通过前后极性间电导率或距离权重来进行精细的空间过滤。以上代码实现的是car算法，这个函数的输入数据data是二维的（通道×时间），函数返回的也是二维数据。
3. 伪迹去除
经常会发现一些漂移伪迹，这些伪迹可以通过比较普通平均假设和初始假设之间的网格方差来检测。伪迹去除需要进行去除或比较。
def drift_removal(data,fs,design='exponential'):
'''
@design：(str)scipy.signal.butter中的滤波函数的类型，包括‘butter’,’chebyshev’,’cheb2’,’ellip’,’bessel’
'''
drift_removal = MFDFA.rolling_band_pass_filter(exponent=0.25,
low_frequency_cutoff=0.01,
high_frequency_cutoff=1,
num_frequencies=50,
sampling_frequency=fs,
filter_type=design,
roll_perc=0.02,
enforce_constant_signal_length=True)
data_removal = drift_removal(np.copy(data))
return data_removal三、信号分析
对脑电文本采用时间序列分析相关方法进行信号分析
1、频谱分析
频谱分析可以帮助我们寻找脑波的周期性。Python可以通过使用pyeeg中的fband函数来实现求功率谱密度，以下是代码示例：
import numpy as np
import pyeeg
import matplotlib.pyplot as plt
def psd_plot(data,fs):
powers, freqs = pyeeg.bin_power(
data,
band=[0.5,4,8,15,30],
fs=fs)
plt.plot(freqs,powers)
plt.show()
2、时频分析
时频分析是指在时间和频率两个维度上分析信号。它是一种用于描述时间变化信号中频率变化的信号分析方法。Python中可以使用tftb中的mmsp函数来实现时频分析，以下是代码示例：
import numpy as np
from tftb.processing import SmoothedPseudoWignerVilleDistribution
def spwvd_transform(data):
n_chls,n_samps = np.shape(data)
spwvd_transform = np.zeros((n_chls,n_samps,n_samps))
for i in range(n_chls):
spwvd_transform[i] = SmoothedPseudoWignerVilleDistribution(data[i])
return spwvd_transform3、小波变换
小波变换是一种基于时间和频率的信号分析方法，它将信号分解为时间和频率两个维度上的不同成分。Python中可以使用pywt来实现小波分析，以下是代码示例：
import pywt
def wavelet_transform(data,wavelet='db3'):
coeffs = []
n_chls,n_samps = np.shape(data)
for i in range(n_chls):
coeffs_chl = pywt.wavedec(data[i],wavelet)
coeffs.append(coeffs_chl)
return coeffs四、特征提取
在时间序列分析中，特征提取是一种用于从信号中提取更有意义信息的方法。使用脑电信号特征提取可以有效地分析脑电数据的变化趋势和特征。下面介绍几种常用的特征提取方法：
1、能量谱密度
脑电信号的功率谱密度（PSD）对于特征提取非常有用。PSD提供了对频率功率分布的深入了解，并可以帮助识别与不同事件相关的特定频率变化。
import numpy as np
from scipy import signal
def psd_feat(data,fs,freqs):
psd_feature = []
window = signal.windows.hann(len(data))
for chl in data:
freqs, psd = signal.welch(chl, fs, window=window,nperseg=len(chl))
psd_feature.append(psd)
return psd_feature2、时域特征
时域特征是对脑电波形的统计特征进行建模。可以使用PyEEG的函数来计算各种时域特征，例如Hurst指数，样品熵等。
import numpy as np
import pyeeg
def time_feat(data):
time_feature = []
for chl in data:
pe = pyeeg.spectral_entropy(chl,6, EEG = True)
time_feature.append(pe)
return psd_feature3、小波特征
小波变换可以实现频率的细分,有助于区分不同类型的脑电信号.通过对各小波系数的统计分析可得到比较鲜明的特征。例如，可以使用Hjorth参数，TKEO等方法计算一些小波形的特征。
import numpy as np
import pywt
def wavelet_feat(data):
feature = []
for chl in data:
cA, cD = pywt.dwt(chl, 'db1')
feature.append([np.mean(cD), np.std(cD), np.mean(abs(cD)), np.median(abs(cD)), np.max(cD)-np.min(cD)])
return feature
