活动介绍
file-type

精神分裂患者静息态脑磁信号的IAF-EMD-样本熵特征提取研究

PDF文件

599KB | 更新于2024-09-01 | 132 浏览量 | 3 评论 | 2 下载量 举报 收藏
download 立即下载
本文主要探讨了基于个体化α峰频(Individual Alpha Frequency, IAF)的静息态脑磁信号特征提取方法在精神分裂患者中的应用。针对精神分裂症患者静息状态下脑磁信号的分析,研究人员提出了一个结合独立成分分析(Independent Component Analysis, ICA)预处理、IAF波段划分(包括快α、慢α1和慢α2)以及经验模式分解(Empirical Mode Decomposition, EMD)和样本熵分析的策略。 首先,对静息态MEG(Magnetoencephalography, 磁共振脑电图)数据进行预处理,利用ICA来减少混叠信号的影响,提高信号质量。ICA通过分解信号为多个独立成分,能够有效地分离出脑磁信号中的潜在源,有助于后续特征分析。 接着,根据IAF理论,将α波段进一步细分为慢α1(IAF的60%到80%)、慢α2(IAF的80%到100%)和快α(IAF到120%)亚频。这是因为不同亚频对应不同的神经功能,如注意力、抑制等功能,对于精神分裂症的认知异常有重要意义。 在特征提取过程中,采用了EMD算法来分析信号的非线性和非稳态特性。EMD将复杂信号分解成一系列独立的IMF成分,每个IMF代表信号的一个特定频率成分,这有助于揭示脑磁信号的内在结构。 样本熵作为一种复杂度度量工具,用于评估信号的时间序列变化规律的不确定性。研究结果显示,精神分裂患者各波段的样本熵普遍高于正常人,尤其在慢α1波段,大脑左半球的额叶、枕叶和颞叶区域表现出显著的异常。这可能反映出精神分裂症患者大脑在静息状态下存在某些认知功能的异常,样本熵的增加可能是认知处理过程中的不稳定性的表现。 本文的工作提供了新的视角来理解和诊断精神分裂症,通过结合IAF、EMD和样本熵,研究者们能够在个体化的基础上更深入地探究脑磁信号在疾病状态下的特征,为精神分裂症的早期识别和治疗提供了可能的技术支持。

相关推荐

filetype

import scipy.io import numpy as np import matplotlib.pyplot as plt from scipy import signal, stats import matplotlib.dates as mdates from datetime import datetime, timedelta import pandas as pd import seaborn as sns import os from pathlib import Path def process_alpha_eeg_data(mat_file_path): """ 处理脑电数据,支持中文路径 参数: mat_file_path: .mat文件路径(支持中文) 返回: fig: 生成的图表对象 session_df: 包含所有会话分析结果的DataFrame """ # 1. 安全加载.mat文件数据(支持中文路径) def safe_loadmat(file_path): """安全加载MAT文件,支持中文路径""" try: # 直接尝试加载 return scipy.io.loadmat(file_path) except UnicodeDecodeError: # 如果失败,使用Path对象转换路径 path_obj = Path(file_path) return scipy.io.loadmat(str(path_obj)) except Exception as e: # 其他错误处理 raise ValueError(f"无法加载MAT文件: {str(e)}") mat_data = safe_loadmat(mat_file_path) # 2. 提取数据矩阵 if 'data' not in mat_data: data_key = [key for key in mat_data.keys() if not key.startswith('__')] if len(data_key) == 1: data_matrix = mat_data[data_key[0]] else: # 尝试搜索可能的矩阵 for key in data_key: if isinstance(mat_data[key], np.ndarray) and mat_data[key].size > 0: data_matrix = mat_data[key] print(f"使用自动检测到的数据矩阵: {key}") break else: raise ValueError("无法确定数据矩阵,请检查.mat文件结构") else: data_matrix = mat_data['data'] # 3. 解析数据矩阵结构 # 假设时间戳在第一列 timestamps = data_matrix[:, 0] # 时间戳 # 检查数据列数 num_columns = data_matrix.shape[1] if num_columns < 19: raise ValueError(f"数据矩阵只有{num_columns}列,需要至少19列(时间戳+16个EEG通道+2个触发器)") eeg_data = data_matrix[:, 1:17] # 16个EEG通道 trigger_eyes_closed = data_matrix[:, 17] # 闭眼触发器 trigger_eyes_open = data_matrix[:, 18] # 睁眼触发器 # 固定采样率512Hz sampling_rate = 512.0 # 4. 高级预处理 def preprocess(data): # 带通滤波提取Alpha波 (8-12Hz) nyquist = 0.5 * sampling_rate low = 8 / nyquist high = 12 / nyquist b, a = signal.butter(4, [low, high], btype='bandpass') alpha_data = signal.filtfilt(b, a, data) # 陷波滤波去除50/60Hz工频干扰 notch_freq = 50.0 notch_width = 2.0 freq = notch_freq / nyquist q = freq / (notch_width / nyquist) b, a = signal.iirnotch(freq, q) return signal.filtfilt(b, a, alpha_data) eeg_data_filtered = np.apply_along_axis(preprocess, 0, eeg_data) # 5. 信号质量检测 def signal_quality(signal_chunk): """检测信号质量""" pp_amplitude = np.max(signal_chunk) - np.min(signal_chunk) zero_crossings = np.sum(np.diff(np.sign(signal_chunk)) != 0) signal_power = np.mean(signal_chunk**2) noise_power = np.var(signal_chunk - signal.savgol_filter(signal_chunk, 51, 3)) snr = 10 * np.log10(signal_power / noise_power) if noise_power > 0 else 100 return pp_amplitude, zero_crossings, snr # 6. 识别所有会话 def find_sessions(trigger): """识别所有会话的开始和结束""" # 找到所有上升沿(会话开始) trigger_diff = np.diff(trigger) session_starts = np.where(trigger_diff == 1)[0] + 1 # 找到所有下降沿(会话结束) session_ends = np.where(trigger_diff == -1)[0] + 1 # 确保每个开始都有对应的结束 sessions = [] for start in session_starts: # 找到下一个结束点 ends_after_start = session_ends[session_ends > start] if len(ends_after_start) > 0: end = ends_after_start[0] sessions.append((start, end)) return sessions # 获取所有会话(闭眼和睁眼) closed_eye_sessions = find_sessions(trigger_eyes_closed) open_eye_sessions = find_sessions(trigger_eyes_open) # 7. 处理每个会话 session_results = [] condition_names = {1: "闭眼", 2: "睁眼"} # 处理闭眼会话 for session_idx, (start_idx, end_idx) in enumerate(closed_eye_sessions): # 提取会话数据 session_eeg = eeg_data_filtered[start_idx:end_idx, :] session_duration = (end_idx - start_idx) / sampling_rate # 计算信号质量指标 quality_metrics = np.apply_along_axis(signal_quality, 0, session_eeg) avg_pp_amplitude = np.mean(quality_metrics[0]) avg_zero_cross = np.mean(quality_metrics[1]) avg_snr = np.mean(quality_metrics[2]) # 计算注意力水平 (使用整个会话数据) # 计算Alpha波能量 (RMS) alpha_energy = np.sqrt(np.mean(session_eeg**2, axis=0)) # 计算注意力水平 (与Alpha能量负相关) channel_attention = 1 / (1 + alpha_energy) channel_attention = (channel_attention - np.min(channel_attention)) / np.ptp(channel_attention) * 100 session_avg_attention = np.mean(channel_attention) # 存储会话结果 session_results.append({ 'session_id': f"C{session_idx+1}", 'condition': 1, 'condition_name': "闭眼", 'start_time': timestamps[start_idx], 'duration': session_duration, 'avg_attention': session_avg_attention, 'pp_amplitude': avg_pp_amplitude, 'zero_crossings': avg_zero_cross, 'snr': avg_snr }) # 处理睁眼会话 for session_idx, (start_idx, end_idx) in enumerate(open_eye_sessions): # 提取会话数据 session_eeg = eeg_data_filtered[start_idx:end_idx, :] session_duration = (end_idx - start_idx) / sampling_rate # 计算信号质量指标 quality_metrics = np.apply_along_axis(signal_quality, 0, session_eeg) avg_pp_amplitude = np.mean(quality_metrics[0]) avg_zero_cross = np.mean(quality_metrics[1]) avg_snr = np.mean(quality_metrics[2]) # 计算注意力水平 (使用整个会话数据) # 计算Alpha波能量 (RMS) alpha_energy = np.sqrt(np.mean(session_eeg**2, axis=0)) # 计算注意力水平 (与Alpha能量负相关) channel_attention = 1 / (1 + alpha_energy) channel_attention = (channel_attention - np.min(channel_attention)) / np.ptp(channel_attention) * 100 session_avg_attention = np.mean(channel_attention) # 存储会话结果 session_results.append({ 'session_id': f"O{session_idx+1}", 'condition': 2, 'condition_name': "睁眼", 'start_time': timestamps[start_idx], 'duration': session_duration, 'avg_attention': session_avg_attention, 'pp_amplitude': avg_pp_amplitude, 'zero_crossings': avg_zero_cross, 'snr': avg_snr }) # 8. 创建会话结果DataFrame session_df = pd.DataFrame(session_results) # 9. 高级可视化 fig = plt.figure(figsize=(18, 20), constrained_layout=True) gs = fig.add_gridspec(4, 2) # 图1: 会话持续时间分布 ax1 = fig.add_subplot(gs[0, 0]) sns.histplot(data=session_df, x='duration', hue='condition_name', bins=10, kde=True, palette=['blue', 'orange'], ax=ax1) ax1.set_title('会话持续时间分布', fontsize=16) ax1.set_xlabel('持续时间 (秒)', fontsize=14) ax1.set_ylabel('频数', fontsize=14) ax1.grid(True, linestyle='--', alpha=0.3) # 图2: 平均注意力比较 ax2 = fig.add_subplot(gs[0, 1]) sns.boxplot(x='condition_name', y='avg_attention', data=session_df, ax=ax2, palette=['blue', 'orange']) ax2.set_title('不同条件下平均注意力比较', fontsize=16) ax2.set_xlabel('条件', fontsize=14) ax2.set_ylabel('平均注意力 (%)', fontsize=14) ax2.grid(True, linestyle='--', alpha=0.3) # 图3: 注意力与持续时间关系 ax3 = fig.add_subplot(gs[1, 0]) sns.scatterplot(x='duration', y='avg_attention', hue='condition_name', data=session_df, ax=ax3, palette=['blue', 'orange'], s=100) ax3.set_title('注意力与会话持续时间关系', fontsize=16) ax3.set_xlabel('持续时间 (秒)', fontsize=14) ax3.set_ylabel('平均注意力 (%)', fontsize=14) ax3.grid(True, linestyle='--', alpha=0.3) # 图4: 信号质量指标 ax4 = fig.add_subplot(gs[1, 1]) sns.scatterplot(x='pp_amplitude', y='snr', hue='condition_name', data=session_df, ax=ax4, palette=['blue', 'orange'], s=100) ax4.set_title('信号峰峰值幅度与信噪比关系', fontsize=16) ax4.set_xlabel('峰峰值幅度 (μV)', fontsize=14) ax4.set_ylabel('信噪比 (dB)', fontsize=14) ax4.grid(True, linestyle='--', alpha=0.3) # 图5: 零交叉率比较 ax5 = fig.add_subplot(gs[2, 0]) sns.boxplot(x='condition_name', y='zero_crossings', data=session_df, ax=ax5, palette=['blue', 'orange']) ax5.set_title('不同条件下零交叉率比较', fontsize=16) ax5.set_xlabel('条件', fontsize=14) ax5.set_ylabel('零交叉率 (次/秒)', fontsize=14) ax5.grid(True, linestyle='--', alpha=0.3) # 图6: 时间-频率分析 (示例闭眼会话) ax6 = fig.add_subplot(gs[2, 1]) if len(closed_eye_sessions) > 0: start_idx, end_idx = closed_eye_sessions[0] sample_session = eeg_data_filtered[start_idx:end_idx, 0] # 取第一个通道 f, t, Sxx = signal.spectrogram( sample_session, fs=sampling_rate, window='hann', nperseg=256, # 0.5秒窗口 noverlap=128 # 50%重叠 ) # 聚焦Alpha波段 (8-12Hz) alpha_mask = (f >= 8) & (f <= 12) # 绘制时频图 im = ax6.pcolormesh(t, f[alpha_mask], 10 * np.log10(Sxx[alpha_mask, :]), shading='gouraud', cmap='viridis') fig.colorbar(im, ax=ax6, label='强度 (dB)') ax6.set_title('闭眼会话Alpha波时频分析', fontsize=16) ax6.set_ylabel('频率 (Hz)', fontsize=14) ax6.set_xlabel('时间 (秒)', fontsize=14) # 图7: 时间-频率分析 (示例睁眼会话) ax7 = fig.add_subplot(gs[3, :]) if len(open_eye_sessions) > 0: start_idx, end_idx = open_eye_sessions[0] sample_session = eeg_data_filtered[start_idx:end_idx, 0] # 取第一个通道 f, t, Sxx = signal.spectrogram( sample_session, fs=sampling_rate, window='hann', nperseg=256, # 0.5秒窗口 noverlap=128 # 50%重叠 ) # 聚焦Alpha波段 (8-12Hz) alpha_mask = (f >= 8) & (f <= 12) # 绘制时频图 im = ax7.pcolormesh(t, f[alpha_mask], 10 * np.log10(Sxx[alpha_mask, :]), shading='gouraud', cmap='viridis') fig.colorbar(im, ax=ax7, label='强度 (dB)') ax7.set_title('睁眼会话Alpha波时频分析', fontsize=16) ax7.set_ylabel('频率 (Hz)', fontsize=14) ax7.set_xlabel('时间 (秒)', fontsize=14) plt.tight_layout() # 10. 保存结果(支持中文路径) base_name = os.path.splitext(os.path.basename(mat_file_path))[0] # 保存图像 fig.savefig(f"{base_name}_analysis.png", dpi=300, bbox_inches='tight') # 保存结果到CSV session_df.to_csv(f"{base_name}_results.csv", index=False, encoding='utf-8-sig') return fig, session_df # 使用示例 if __name__ == "__main__": # 中文路径示例 mat_file = "F:/Grade2/attention/2348892/subject_02.mat" try: # 处理数据并生成可视化 fig, results_df = process_alpha_eeg_data(mat_file) # 显示图像 plt.show() # 打印会话结果 print("EEG分析结果:") print(results_df[['session_id', 'condition_name', 'duration', 'avg_attention']]) # 计算统计显著性 closed_eyes = results_df[results_df['condition'] == 1] open_eyes = results_df[results_df['condition'] == 2] if not closed_eyes.empty and not open_eyes.empty: t_stat, p_value = stats.ttest_ind(closed_eyes['avg_attention'], open_eyes['avg_attention']) print(f"\n闭眼与睁眼条件间注意力差异显著性检验: t = {t_stat:.2f}, p = {p_value:.4f}") # 计算条件间差异百分比 closed_avg = closed_eyes['avg_attention'].mean() open_avg = open_eyes['avg_attention'].mean() diff_percent = ((open_avg - closed_avg) / closed_avg) * 100 print(f"睁眼条件相对于闭眼条件的注意力提升: {diff_percent:.1f}%") else: print("\n警告:缺少闭眼或睁眼会话数据,无法进行统计检验") except Exception as e: print(f"处理过程中发生错误: {str(e)}") # 可以添加详细的错误日志记录 import traceback traceback.print_exc() 分析上述代码

资源评论
用户头像
士多霹雳酱
2025.08.16
本研究深入分析了精神分裂症患者的脑磁信号特征,有效区分了正常与异常脑波段差异。
用户头像
XiZi
2025.04.21
该方法能准确捕捉精神分裂症患者在静息态下的脑活动异常,具有临床应用潜力。
用户头像
方2郭
2025.04.14
使用EMD和样本熵方法提取脑磁信号特征,为精神疾病诊断提供新视角。🌈
weixin_38738977
  • 粉丝: 6
上传资源 快速赚钱