脑电信号4

脑电信号4目录1. 盲源分离——从鸡尾酒会问题说起2. ICA的数学原理3. ICA在脑电中的实战应用4. EEG微状态——大脑的基本字母5. 微状态分析流程6. Python完整实战1. 盲源分离——从鸡尾酒会问题说起BSS盲源分离就是仅从混合录音中分离各自的声音不知道混合比例也不知道源信号长什么样。脑电跟BSS一样头皮64个电极录到的是脑内多个源的混合。BSS帮我们把它们拆开。importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportsignalfromsklearn.decompositionimportFastICA# 直观演示两个源→混合→ICA还原np.random.seed(42)tnp.linspace(0,2,500)s1np.sin(2*np.pi*2*t)# 源1慢正弦模拟脑电s2signal.sawtooth(2*np.pi*8*t)# 源2快锯齿模拟眼电x10.8*s10.5*s2# 电极1混合x20.3*s10.9*s2# 电极2混合Xnp.c_[x1,x2]icaFastICA(n_components2,random_state42,max_iter1000)S_hatica.fit_transform(X)fig,axesplt.subplots(2,2,figsize(14,6))axes[0,0].plot(t,s1,b,lw1.5);axes[0,0].plot(t,s2,r,lw1.5)axes[0,0].set_title(真实信号源);axes[0,0].legend([脑电,眼电])axes[0,1].plot(t,x1,gray,lw0.8);axes[0,1].plot(t,x2,k,lw0.8)axes[0,1].set_title(头皮电极录到的混合信号);axes[0,1].legend([电极1,电极2])axes[1,0].plot(t,S_hat[:,0],b,lw1.5);axes[1,0].plot(t,S_hat[:,1],r,lw1.5)axes[1,0].set_title(ICA还原的独立成分)axes[1,1].imshow(ica.mixing_,cmapRdBu_r,aspectauto)axes[1,1].set_title(混合矩阵A);axes[1,1].set_xticks([0,1]);axes[1,1].set_yticks([0,1])plt.suptitle(ICA盲源分离仅从头皮混合信号还原出独立源,fontsize13)plt.tight_layout();plt.show()ICA vs PCA——很多人混淆这两个PCAICA目标找方差最大的方向找统计独立的方向假设信号服从高斯分布信号非高斯输出成分之间不相关正交成分之间统计独立更强2. ICA的数学原理ICA模型X A·S。已知X头皮EEG求W≈A⁻¹得 Ŝ W·X。核心假设源信号是非高斯的。为什么中心极限定理——多个独立非高斯信号的混合比单个信号更像高斯分布。ICA做的就是找到让输出最不像高斯分布的分离方向——“非高斯性≈独立性”。用负熵量化非高斯性J(y) H(y_gauss) − H(y)。高斯分布熵最大负熵越大→越不像高斯→越可能是独立源。FastICA通过最大化负熵迭代求解W。defdemo_ica_geometry():直观感受ICA的几何意义np.random.seed(0)s1np.random.uniform(-1,1,1000)# 两个独立均匀源s2np.random.uniform(-1,1,1000)Xnp.array([[0.7,0.5],[0.3,0.8]]) np.vstack([s1,s2])S_hatFastICA(n_components2,random_state42).fit_transform(X.T).T fig,axesplt.subplots(1,3,figsize(15,4.5))axes[0].scatter(s1,s2,s3,alpha0.5)axes[0].set_title(源信号正方形独立);axes[0].set_aspect(equal)axes[1].scatter(X[0],X[1],s3,alpha0.5,cred)axes[1].set_title(混合后菱形相关);axes[1].set_aspect(equal)axes[2].scatter(S_hat[0],S_hat[1],s3,alpha0.5,cgreen)axes[2].set_title(ICA还原恢复正方形);axes[2].set_aspect(equal)plt.suptitle(ICA几何直觉混合旋转拉伸→ICA解旋还原,fontsize13)plt.tight_layout();plt.show()demo_ica_geometry()3. ICA在脑电中的实战应用3.1 三大去伪迹场景伪迹类型ICA识别依据关键特征眼电(EOG)前额通道权重高 低频功率大波形呈尖峰、峰度5肌电(EMG)高频功率大(30Hz)颞区权重高、不规则心电(ECG)规律尖峰(~1Hz)空间分布广泛3.2 自动分类成分defclassify_ica_components(ica,raw_eeg,fs256):自动分类ICA成分眼电/肌电/心电/脑电n_compica.mixing_.shape[1]sourcesica.transform(raw_eeg.T).T mixingica.mixing_ labels[]foriinrange(n_comp):frontal_wnp.abs(mixing[:3,i]).mean()/(np.abs(mixing[:,i]).mean()1e-10)f,psdsignal.welch(sources[i],fs,npersegfs*2)totalnp.trapz(psd,f)1e-10low_rnp.trapz(psd[(f0.5)(f4)],f[(f0.5)(f4)])/total high_rnp.trapz(psd[(f30)(f50)],f[(f30)(f50)])/total kurtnp.mean((sources[i]-sources[i].mean())**4)/(np.std(sources[i])**41e-10)iffrontal_w1.5andlow_r0.3andkurt5:labels.append(眼电)elifhigh_r0.3:labels.append(肌电)elif0.05low_r0.4andkurt4:labels.append(心电)else:labels.append(脑电)returnlabels使用ICA的注意事项先做1Hz高通滤波再做ICA——低频漂移会让分解不稳定成分数≤通道数数据越长ICA越稳定自动标记后建议人工复查——误删脑电成分比漏删伪迹更糟糕去坏段后再ICA——运动伪迹会污染所有成分4. EEG微状态——大脑的基本字母4.1 什么是微状态观察EEG的多通道地形图你会发现头皮电位不是随机变化的——它在几种固定空间模式之间快速切换每种模式持续约60-120ms。这些准稳态的地形图就是微状态Microstate。打个比方微状态像字母大脑用有限的字母约4种的不同排列组合来表达不同认知状态。4.2 经典四类微状态微状态地形图特征功能关联A类左枕→右额梯度听觉/语言网络B类右枕→左额梯度视觉网络C类额-枕对称中线突显网络——自我参照、冥想相关D类中央-额叶分布注意网络——任务正相关4.3 四个关键时间参数参数含义冥想中的变化持续时间平均持续多少ms冥想相关微状态持续时间↑覆盖率占总时间百分比冥想相关微状态占比↑出现频率每秒出现次数切换频率↓→大脑更稳定转移概率从A切换到B的概率某些转移路径概率改变有研究发现资深冥想者的C类微状态覆盖率显著高于普通人且切换频率更低——大脑更倾向于稳定地保持在自我参照状态。5. 微状态分析流程5.1 四步走预处理EEG → 计算GFP → 取GFP峰值处地形图 → K-means聚类(4类) → 回标到原始EEG取GFP峰值处的原因GFP峰值处信噪比最高、地形图最稳定。谷值处地形图不可靠。deffind_gfp_peaks(eeg_multi,fs256,min_dist_ms50):找到GFP峰值时刻gfpnp.std(eeg_multi,axis0)peakssignal.find_peaks(gfp,distanceint(min_dist_ms/1000*fs))[0]returnpeaks,gfpdefextract_microstates(eeg_data,n_states4,n_peaks500,fs256):从EEG提取微状态GFP峰值→聚类→回标# 1. 找GFP峰值处的地形图peaks,gfpfind_gfp_peaks(eeg_data,fs)top_idxpeaks[np.argsort(gfp[peaks])[-n_peaks:]]toposeeg_data[:,top_idx].T# 2. K-means聚类fromsklearn.clusterimportKMeans kmeansKMeans(n_clustersn_states,random_state42,n_init10)kmeans.fit(topos)# 3. 回标每个时间点分配给最近聚类中心eeg_teeg_data.T distsnp.array([np.sum((eeg_t-c)**2,axis1)forcinkmeans.cluster_centers_])seqnp.argmin(dists,axis0)returnseq,kmeans.cluster_centers_defmicrostate_params(seq,fs256,min_dur_ms30):计算微状态时间参数min_sint(min_dur_ms/1000*fs)params{}forsinnp.unique(seq):is_s(seqs).astype(int)edgesnp.diff(np.concatenate([[0],is_s,[0]]))startsnp.where(edges1)[0]endsnp.where(edges-1)[0]durs(ends-starts)/fs*1000validdursmin_s params[fMS{s}]{平均持续(ms):np.mean(durs[valid]),覆盖率(%):is_s.sum()/len(seq)*100,频率(次/s):len(starts[valid])/(len(seq)/fs),}returnparams6. Python完整实战classEEGBSSAnalyzer:ICA去伪迹 微状态分析def__init__(self,fs256):self.fsfs;self.icaNonedefrun_ica(self,eeg):n_cheeg.shape[0]self.icaFastICA(n_componentsn_ch,random_state42,max_iter2000)sourcesself.ica.fit_transform(eeg.T).Treturnsources,self.ica.mixing_defremove_artifacts(self,eeg,artifact_idx):sourcesself.ica.transform(eeg.T).T sources[artifact_idx]0returnself.ica.mixing_ sources# 完整演示 fs,dur256,10tnp.arange(0,dur,1/fs);nlen(t);n_ch8np.random.seed(42)# 生成模拟多通道EEGeegnp.array([(10-ch*0.8)*np.sin(2*np.pi*10*tch*0.3)(5ch*0.3)*np.sin(2*np.pi*20*tch*0.5)np.random.randn(n)*3forchinrange(n_ch)])# 加眼电前额通道强forbtin[2.5,5.0,7.5]:blink100*np.exp(-((t-bt)/0.1)**2)forchinrange(n_ch):eeg[ch]blink*(1.5-0.18*ch)# 1. ICA去眼电analyzerEEGBSSAnalyzer(fs)sources,mixinganalyzer.run_ica(eeg)artifact_idxnp.argmax(np.abs(mixing[:2]).mean(axis0))# 前额权重最高的成分eeg_cleananalyzer.remove_artifacts(eeg,[artifact_idx])# 2. 微状态分析seq,centersextract_microstates(eeg_clean,n_states4,fsfs)paramsmicrostate_params(seq,fs)# 可视化fig,axesplt.subplots(2,2,figsize(14,8))# ICA前后对比axes[0,0].plot(t,eeg[0],gray,lw0.5,alpha0.7,label含眼电)axes[0,0].plot(t,eeg_clean[0],blue,lw1,label去眼电后)axes[0,0].set_title(f通道0(Fp1) ICA前后对比);axes[0,0].legend();axes[0,0].grid(True,alpha0.3)# 混合矩阵imaxes[0,1].imshow(mixing,aspectauto,cmapRdBu_r)axes[0,1].set_title(f混合矩阵(成分{artifact_idx}眼电))axes[0,1].axvline(artifact_idx,colorred,lw2)plt.colorbar(im,axaxes[0,1])# 微状态覆盖率nameslist(params.keys())axes[1,0].bar(names,[params[k][覆盖率(%)]forkinnames],color[red,blue,green,orange],edgecolorblack)axes[1,0].set_title(四类微状态覆盖率);axes[1,0].set_ylabel(%)# 微状态时序segseq[500:900]colors[red,blue,green,orange]axes[1,1].fill_between(range(len(seg)),seg-0.4,seg0.4,color[colors[int(s)]forsinseg],alpha0.5)axes[1,1].plot(seg,black,lw0.5)axes[1,1].set_title(微状态时序~1.6秒片段);axes[1,1].set_ylabel(状态)axes[1,1].set_ylim(-0.5,3.5);axes[1,1].set_yticks([0,1,2,3])plt.suptitle(ICA去眼电 微状态聚类 完整演示,fontsize14)plt.tight_layout();plt.show()print(\n微状态时间参数:)forkinsorted(params.keys()):pparams[k]print(f{k}: 持续{p[平均持续(ms)]:.0f}ms, 覆盖{p[覆盖率(%)]:.1f}%, 频率{p[频率(次/s)]:.2f}次/s)关键发现ICA能有效分离眼电成分前额高权重低频尖峰特征明显EEG微状态通常聚类为4类其中C类与自我参照/冥想状态密切相关微状态时间参数比传统频谱指标更稳定——是个体脑指纹的候选去伪迹后脑电成分保留完好重构信号波形自然无失真