脑电 6目录1. 为什么需要非线性分析2. 熵——信号不可预测性的度量3. 复杂度——信号中模式丰富度的量化4. 分形与自相似性5. 去趋势波动分析——长程相关性6. 非线性特征工具箱1. 为什么需要非线性分析1.1 脑电不是线性系统线性系统的典型特征输入和输出成比例不同频率成分互不干扰。但大脑是典型的非线性系统两个刺激同时呈现 ≠ 单独呈现的EEG响应之和α波和γ波之间存在跨频率耦合γ功率随α相位变化癫痫从正常→发作的转变是突然的相变非线性动力学中的分岔现象同一被试在相同状态下的EEG波形从不完全重复importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportsignal# 直观演示线性系统 vs 非线性系统tnp.linspace(0,4,500)xnp.sin(2*np.pi*2*t)0.5*np.sin(2*np.pi*8*t)# 线性系统y 2*x简单放大y_linear2*x# 非线性系统y x³输出不是输入的简单缩放y_nonlinearx**3fig,axesplt.subplots(2,2,figsize(14,7))axes[0,0].plot(t,x,b,lw1);axes[0,0].set_title(输入信号 x(t))axes[0,1].plot(t,y_linear,g,lw1);axes[0,1].set_title(线性系统输出 y2x\n(形状不变))axes[1,0].plot(t,y_nonlinear,r,lw1);axes[1,0].set_title(非线性系统输出 yx³\n(波形变形、高频增加))# 输入-输出散点图axes[1,1].scatter(x,y_linear,s2,alpha0.3,label线性(y2x))axes[1,1].scatter(x,y_nonlinear,s2,alpha0.3,label非线性(yx³))axes[1,1].set_xlabel(输入);axes[1,1].set_ylabel(输出)axes[1,1].legend();axes[1,1].set_title(输入-输出关系)plt.tight_layout();plt.show()print(线性输入×2 输出×2比例恒定)print(非线性输入和输出的关系不是直线——这就是非线性分析的动机)1.2 非线性指标能捕捉的东西线性指标非线性指标新增视角α功率样本熵α功率不变但熵降低→信号更有序了ERP峰值LZ复杂度峰值不变但复杂度降低→进入单一模式θ/β比值分形维度比值不变但D降低→信号更光滑规则频谱分布去趋势波动分析(DFA)频谱不变但长程相关性增强2. 熵——信号不可预测性的度量2.1 三种熵的对比熵类型含义冥想中的变化近似熵(ApEn)新模式出现的概率冥想加深→ApEn↓样本熵(SampEn)ApEn的改进版排除了自匹配比ApEn更稳定长度无关排列熵(PermEn)仅看大小排序模式计算极快对噪声鲁棒适合在线分析defsample_entropy(data,m2,r_factor0.2): 样本熵——衡量信号中新模式产生的概率 熵高信号不可预测混沌/警觉 熵低信号规律有序/冥想/麻醉 N,rlen(data),r_factor*np.std(data)defcount_matches(template_len):count0templatesnp.array([data[i:itemplate_len]foriinrange(N-template_len)])foriinrange(len(templates)):distsnp.max(np.abs(templates-templates[i]),axis1)countnp.sum(distsr)-1# 减1排除自身匹配returncount A,Bcount_matches(m1),count_matches(m)return-np.log(A/(B1e-10))ifB0elsefloat(inf)defpermutation_entropy(data,m3,delay1): 排列熵——只看数值大小排序模式 例窗口[1,5,3]→排序[0,2,1]→中-小-大模式 对所有窗口统计模式分布→算熵 fromitertoolsimportpermutations Nlen(data)patternslist(permutations(range(m)))pattern_count{p:0forpinpatterns}total0foriinrange(N-(m-1)*delay):windowdata[i:im*delay:delay]# 找到排序模式sorted_idxtuple(np.argsort(window))ifsorted_idxinpattern_count:pattern_count[sorted_idx]1total1# 计算熵pe0forcountinpattern_count.values():ifcount0:pcount/total pe-p*np.log(p)returnpe/np.log(len(patterns))# 归一化到[0,1]# 对比演示np.random.seed(42)noisenp.random.randn(1000)# 白噪声→高熵sinenp.sin(2*np.pi*np.arange(1000)/256*10)# 正弦→低熵mixednoise*0.6sine*0.4# 混合→中等print(信号类型对比:)forname,sigin[(白噪声(高熵),noise),(正弦波(低熵),sine),(混合(中熵),mixed)]:sesample_entropy(sig)pepermutation_entropy(sig[:500])print(f{name}: SampleEn{se:.3f}, PermEn{pe:.3f})3. 复杂度——信号中模式丰富度的量化3.1 Lempel-Ziv复杂度LZC衡量的不是信号有多乱而是信号中有多少种不同的子串模式——模式越多越复杂。deflempel_ziv_complexity(data):LZ复杂度——计算信号二值化后不同子串模式数量# 二值化大于中位数1小于0binary(datanp.median(data)).astype(int)nlen(binary)i,C,L0,1,1# C复杂度计数器, L当前窗口长度whileiLn:pattern.join(map(str,binary[i:iL]))history.join(map(str,binary[:iL-1]))ifpatterninhistory:L1# 模式出现过→扩大窗口else:C1# 新模式→复杂度1iL L1returnC/n# 归一化# 直观演示规则序列 vs 随机序列regularnp.tile([1,0,1,0],50)# 高度规则→低复杂度random_seqnp.random.randint(0,2,200)# 随机→高复杂度print(f规则序列(1010重复) LZC:{lempel_ziv_complexity(regular):.3f})print(f随机序列 LZC:{lempel_ziv_complexity(random_seq):.3f})print(→ 冥想深度↑ → 大脑活动模式趋于简单 → LZC↓)3.2 LZC在脑电中的应用麻醉深度监测LZC随麻醉加深而持续降低——比频谱指标更可靠冥想状态评估深度冥想时LZC显著低于静息状态癫痫预测发作前几分钟LZC异常降低——可用于预警意识障碍评估植物状态患者的LZC显著低于最小意识状态4. 分形与自相似性4.1 脑电是分形信号分形的核心特征自相似性——放大看局部和整体有相似的统计特性。EEG信号在毫秒到分钟的多个时间尺度上都呈现1/f频谱功率∝1/f^α这正是分形信号的特征。defhiguchi_fractal_dimension(data,kmax10): Higuchi分形维度——衡量信号的曲折程度 D范围1(光滑规则曲线) ~ 2(极其曲折填满平面) 冥想加深→D降低信号更光滑规则 麻醉加深→D降低 癫痫发作→D降低 Nlen(data)Lknp.zeros(kmax)forkinrange(1,kmax1):Lmknp.zeros(k)forminrange(k):idxsnp.arange(m,N-1,k)iflen(idxs)1:Lmk[m]np.sum(np.abs(np.diff(data[idxs])))Lmk[m]*(N-1)/((len(idxs)-1)*k)validLmk[Lmk0]Lk[k-1]np.mean(valid)iflen(valid)0else0# 最小二乘拟合: log(Lk) ~ log(1/k)xnp.log(1.0/np.arange(1,kmax1))ynp.log(Lk1e-10)D-np.polyfit(x,y,1)[0]returnDdefpetrosian_fractal_dimension(data):Petrosian分形维度——基于零交叉点计数计算极快# 零交叉信号穿过均值的次数zero_crossingsnp.sum(np.diff(datanp.mean(data))!0)Nlen(data)returnnp.log10(N)/(np.log10(N)np.log10(N/(N0.4*zero_crossings)))5. 去趋势波动分析——长程相关性5.1 DFA的原理DFADetrended Fluctuation Analysis衡量信号在不同时间尺度上的波动特征核心输出是一个指数αα ≈ 0.5白噪声——各时间尺度随机无长程相关α ≈ 1.01/f噪声粉红噪声——健康大脑的典型特征α 1.0存在长程正相关——前一刻的趋势延续到下一刻α 0.5反相关——涨落后倾向于被跌落跟随defdfa(data,scalesNone): 去趋势波动分析(DFA) 返回缩放指数α——衡量长程相关性强度 ifscalesisNone:scalesnp.unique(np.logspace(1,np.log10(len(data)//4),10).astype(int))# 1. 累积和积分ynp.cumsum(data-np.mean(data))# 2. 对每个尺度计算波动fluctuationsnp.zeros(len(scales))foridx,scaleinenumerate(scales):n_segmentslen(data)//scale rms_total0foriinrange(n_segments):segy[i*scale:(i1)*scale]# 用多项式拟合去趋势线性趋势DFA1xnp.arange(len(seg))trendnp.polyval(np.polyfit(x,seg,1),x)detrendedseg-trend rms_totalnp.mean(detrended**2)fluctuations[idx]np.sqrt(rms_total/n_segments)# 3. 拟合log(F) vs log(scale) → 斜率αvalid~np.isnan(fluctuations)(fluctuations0)alphanp.polyfit(np.log10(scales[valid]),np.log10(fluctuations[valid]),1)[0]returnalpha,scales,fluctuations# 演示不同α的对比np.random.seed(42)white_noisenp.random.randn(2000)# α≈0.5# 生成1/f噪声简化的Pink噪声pink_noisenp.cumsum(np.random.randn(2000))# α≈1.5布朗噪声forname,sigin[(白噪声(α≈0.5),white_noise),(布朗噪声(α≈1.5),pink_noise)]:alpha,_,_dfa(sig)print(f{name}: DFA α {alpha:.3f})print(健康静息EEG的DFA α通常在0.7-1.0之间)print(深度麻醉/昏迷时α可能下降到0.5左右)5.2 DFA在临床和研究中的应用应用DFA α的变化麻醉监测α从1.0(清醒)降到~0.6(深度麻醉)癫痫预测发作前α异常升高→长程相关性增强阿尔茨海默病α低于健康同龄人——脑信号失去长程组织冥想α略增——大脑维持稳定状态的能力增强睡眠分期NREM睡眠αREM睡眠α清醒α6.非线性特征工具箱classNonlinearEEGAnalyzer:非线性EEG特征提取工具箱def__init__(self,fs256):self.fsfsdefextract_all(self,eeg_segment):一键提取所有非线性特征features{}# 熵类features[sample_entropy]sample_entropy(eeg_segment)features[perm_entropy]permutation_entropy(eeg_segment,m3)# 复杂度features[lz_complexity]lempel_ziv_complexity(eeg_segment)# 分形features[higuchi_fd]higuchi_fractal_dimension(eeg_segment)features[petrosian_fd]petrosian_fractal_dimension(eeg_segment)# 长程相关alpha,_,_dfa(eeg_segment)features[dfa_alpha]alphareturnfeaturesdefcompare_states(self,eeg_segments,state_names):对比不同状态的非线性特征results[]forseg,nameinzip(eeg_segments,state_names):results.append(self.extract_all(seg))# 可视化feat_nameslist(results[0].keys())n_featlen(feat_names)fig,axesplt.subplots(1,n_feat,figsize(3*n_feat,4))fori,(feat,ax)inenumerate(zip(feat_names,axes)):values[r[feat]forrinresults]ax.bar(range(len(state_names)),values,color[steelblue,salmon,orange,green][:len(values)],edgecolorblack)ax.set_xticks(range(len(state_names)))ax.set_xticklabels(state_names,fontsize7)ax.set_title(feat,fontsize9);ax.grid(True,alpha0.3,axisy)plt.suptitle(不同脑状态下非线性特征对比,fontsize13)plt.tight_layout();plt.show()returnresults# 演示四种模拟状态analyzerNonlinearEEGAnalyzer(fs256)np.random.seed(42)tnp.arange(0,4,1/256)# 模拟四种状态rest20*np.sin(2*np.pi*10*t)np.random.randn(len(t))*3meditation25*np.sin(2*np.pi*10*t)np.random.randn(len(t))*1.5# 更干净规则drowsy10*np.sin(2*np.pi*5*t)np.random.randn(len(t))*5# 更多慢波噪声noise_signp.random.randn(len(t))*8# 纯噪声resultsanalyzer.compare_states([rest,meditation,drowsy,noise_sig],[静息,深度冥想,困倦,随机噪声])print(\n各状态特征汇总:)fori,(name,r)inenumerate(zip([静息,深度冥想,困倦,随机噪声],results)):print(f\n{name}:)fork,vinr.items():print(f{k}:{v:.4f})ifi1:print( → 深度冥想熵最低、复杂度最低、分形维最低最有序)
