一、腦電信號獲取
腦電信號是由腦部神經元的興奮和抑制行為產生的微弱電流信號。為了獲取腦電信號,需要使用電極陣列將電極頭置於頭皮上。腦電信號是一種高噪聲、低幅值、高時間分辨率的信號,因此需要進行信號放大和濾波。
二、腦電信號預處理
腦電信號需要進行預處理處理,以消除噪聲、偽跡和其他干擾。預處理包括:濾波、空間濾波、偽跡去除。
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_transform
3、小波變換
小波變換是一種基於時間和頻率的信號分析方法,它將信號分解為時間和頻率兩個維度上的不同成分。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_feature
2、時域特徵
時域特徵是對腦電波形的統計特徵進行建模。可以使用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_feature
3、小波特徵
小波變換可以實現頻率的細分,有助於區分不同類型的腦電信號.通過對各小波係數的統計分析可得到比較鮮明的特徵。例如,可以使用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
原創文章,作者:小藍,如若轉載,請註明出處:https://www.506064.com/zh-hant/n/303396.html