当前位置： 首页 > ops >正文

77、完全植入式脑机接口神经数据编码解码数据处理等问题答疑[嘿！快看，馒头老师在这里蹲着！]

ops 2025/6/29 13:39:30

       今天出个外勤，对接项目事宜。整理了一下自己学生问的一些问题，我看了一下汇总的内容挺多的，自己做侵入式的学生不多，现在只有两位学子，而且全侵入Spikes数据的处理问题网上相关资料较少，所以我把我们在复现一些侵入式论文过程中遇到的问题和各位学子分享一下，内容较多，此次分享分三次博客撰写。

                                                                                            —— 馒头   2025.6.4日于北京畅春园

研究背景：

做神经信号处理pipeline和indy数据集分析，在尝试用indy数据集里的原始数据跑通spikesorting全流程，遇到了以下问题。

问题1：学生在用spikeinterface时需要引入探针信息，如何找到或自定义与研究范式一致的探针类？（10*10的方形电极，但四个角落没有功能，96通路）

问题需求分析：

  • 探针结构：10×10 方形网格，总 100 个电极位置，但四个角落电极（坐标 (0,0), (0,9), (9,0), (9,9)）不工作，实际有效电极为 96 个。

  • 方法：用 probeinterface 生成该探针布局，并绑定到 SpikeInterface 的 Recording 对象。

解决方法：

第一步，定义探针电极：create_10x10_custom_probe.py

"""
代码步骤：
1. 使用ProbeInterface生成10x10的探针。
2. 确定四个角落电极的索引并排除。
3. 应用掩码，创建子集探针。
4. 可视化探针布局以确认正确性。
5. 将探针绑定到Recording对象。
"""
import numpy as np
from probeinterface import Probe, generate_multi_columns_probe
from probeinterface.plotting import plot_probe
import matplotlib.pyplot as pltdef create_10x10_custom_probe():"""创建一个 10x10 方形电极阵列，排除四个角落电极，共 96 个有效通道。"""# === 步骤 1：生成 10x10 探针 ==="""生成一个 10 列 × 10 行的电极阵列，电极间距为 40μm。电极按 列优先 顺序排列（第0列电极索引为 0~9，第1列为 10~19，依此类推）。"""probe = generate_multi_columns_probe(num_columns=10,                 # 10 列num_contact_per_column=10,      # 每列 10 个电极xpitch=40,                      # 电极水平间距 40μmypitch=40,                      # 电极垂直间距 40μmcontact_shapes='circle',        # 电极形状为圆形contact_shape_params={'radius': 5}  # 半径 5μm)# === 步骤 2：标记要排除的四个角落电极 ===# 假设电极按列优先排列（第0列：0~9，第1列：10~19，...，第9列：90~99）"""根据列优先排列规则，四个角落电极的索引为 [0, 9, 90, 99]。"""corner_indices = [0,   # 第0列第一个电极 (0,0)9,    # 第0列最后一个电极 (0,9)90,   # 第9列第一个电极 (9,0)99    # 第9列最后一个电极 (9,9)]# === 步骤 3：创建掩码排除角落电极 ==="""通过布尔掩码 (mask) 将这些电极标记为无效，生成子探针 probe_subset。"""n_contacts = len(probe.contact_ids)mask = np.ones(n_contacts, dtype=bool)mask[corner_indices] = False  # 将角落电极标记为无效# === 步骤 4：生成子探针（仅保留有效电极）===probe_subset = probe.get_subsets(mask)# === 步骤 5：定义探针物理形状（可选）==="""create_auto_shape(probe_type='rect') 根据电极位置自动生成探针的物理边界（矩形）。"""probe_subset.create_auto_shape(probe_type='rect', margin=20)return probe_subset# === 测试代码 ===
if __name__ == "__main__":# 创建探针custom_probe = create_10x10_custom_probe()# 打印探针信息print(f"电极数量: {len(custom_probe.contact_ids)}")  # 应输出 96# 可视化探针布局fig, ax = plt.subplots(figsize=(8, 8))plot_probe(custom_probe, ax=ax, with_channel_index=True)ax.set_title("10x10 方形探针（排除四个角落）")plt.show()#代码输出应为 电极数量: 96

第二步，绑定到 SpikeInterface 的 Recording 对象：也就是把自定义探针关联到数据

import spikeinterface.extractors as se# 假设已加载原始数据 recording
recording = se.read_xxx(...)  # 替换为实际数据加载代码# 绑定探针
recording_with_probe = recording.set_probe(custom_probe, in_place=False)# 验证探针信息
print(recording_with_probe.get_probe())

馒头老师总结：

上述方法可精确生成一个 10×10 方形电极阵列（排除四个角落） 的探针布局，并集成到 SpikeInterface 流程中。此方法直接解决了探针几何信息与实验范式不匹配的问题。此外，还可能遇到[电极索引错误|探针形状不匹配|数据通道数不匹配]的问题，对此遇到了后面也很好解决，在此先不提。

问题2 ：学生如何选择适合数据的Sorting算法和滤波器以及超参数

问题需求分析：

【学生需要学习如何选择适合的 Spike Sorting 算法 和 滤波器参数，以及如何设置超参数。下面我按照这3个需求进行解答】

2.1.选择 Sorting 算法

推荐算法

根据 Indy 数据集 的特点（高密度探针，如 Neuropixels），以下是我推荐的排序算法：

    算法                           适用场景                           优点
Kilosort 2/3高密度探针（如 Neuropixels），需 GPU 加速速度快，适合大规模数据，自动聚类效果好
Mountainsort 4自动化流程，无需 GPU，适合中低密度数据无需 GPU，参数简单，适合快速测试
SpyKING                  CIRCUS多电极阵列（MEA）数据，基于模板匹配对噪声鲁棒，适合复杂背景下的尖峰检测
Tridesclous在线排序或自动化流程，内置参数优化自动化程度高，适合无需手动干预的场景

馒头老师建议:

  • 首选：Kilosort 3（需 GPU 支持,也是做神经科学必备算法）。

  • 备选：Mountainsort 4（无需 GPU，适合快速测试）。

2.2.滤波器选择与配置

[不同滤波器的使用和参数设置具体需要以数据为例来分析，下面先简单介绍一些基础设置，然后对于学生写的代码的问题为例，进行进一步的解析。]

1.滤波器

尖峰信号通常位于 300 Hz 至 6000 Hz 频段，以下滤波器基本配置：

1.带通滤波器：提取尖峰信号。频率范围：300 Hz 至 6000 Hz。滤波器类型：Butterworth（二阶）。代码示例：recording_filtered = spre.bandpass_filter(recording, freq_min=300, freq_max=6000)

2.去工频干扰（可选做）：去除 50/60 Hz 电源噪声。频率：50 Hz 或 60 Hz（根据实验环境）。滤波器类型：Notch 滤波器。代码示例：recording_filtered = spre.notch_filter(recording_filtered, freq=50)

3.重参考：减少公共噪声。方法：全局中值参考（Global Median Reference）。代码：recording_filtered = spre.common_reference(recording_filtered, reference='global', operator='median')

2.3.超参数设置

Kilosort 3 推荐参数配置：

参数推荐值说明
freq_min300带通滤波下限频率（Hz）
freq_max6000带通滤波上限频率（Hz）
detect_threshold5尖峰检测阈值（噪声标准差的倍数）
carTRUE是否启用公共平均参考（Common Average Reference）
nblocks5分块处理的数量，影响内存占用和速度
nt061尖峰波形的时间点数（默认值）
whitening_range32白化邻域范围（通道数）
min_spikes100每个单元的最小尖峰数（用于聚类）

我们复现了处理LFP+SUA+MUA+ESA信号的代码，下面以学生写的代码的问题为例，指出问题所在和修改原因和建议。

1. LFP（局部场电位）处理

  • 代码中的设置：

    • 使用 1-300 Hz 带通滤波。

  • 图片中的设置：

    • 使用低通滤波（LPF），截止频率为 300 Hz。

  • 问题：

    • 代码中使用的是带通滤波（1-300 Hz），而图片中建议的是低通滤波（LPF）到 300 Hz。

  • 修改建议：

    • 将 preprocess_lfp 函数中的带通滤波改为低通滤波。

  • 修改原因：

    • LFP 信号主要是低频信号，通常只需要低通滤波来去除高频噪声，带通滤波可能会引入不必要的低频成分。

修改后的代码：

def preprocess_lfp(recording, probe):logging.info("Preprocessing LFP: lowpass filter (300 Hz) and common reference")recording_filtered = spre.lowpass_filter(recording, freq_max=300)  # 改为低通滤波recording_cmr = spre.common_reference(recording_filtered, reference='global', operator='median')recording_cmr = recording_cmr.set_probe(probe, in_place=False)return recording_cmr

2. SUA（单元活动）处理

  • 代码中的设置：

    • 使用 Mountainsort4 进行 spike sorting，带通滤波范围为 250-5000 Hz。

  • 图片中的设置：

    • 带通滤波范围为 500-5000 Hz（对于数据集 I）或 250-5000 Hz（对于数据集 II）。

  • 问题：

    • 代码中的带通滤波范围是 250-5000 Hz，而图片中对于数据集 I 的建议是 500-5000 Hz。

  • 修改建议：

    • 根据数据集类型调整带通滤波范围。如果使用的是数据集 I，应将带通滤波范围改为 500-5000 Hz。

  • 修改原因：

    • 不同的数据集可能需要不同的滤波范围来优化 spike sorting 的效果。

修改后的代码：

def run_sua_sorting(recording_saved):logging.info("==== SUA Sorting Process Start ====")try:import spikeinterface.sorters as sis# 设置临时文件目录（根据实际情况调整）os.environ['TEMPDIR'] = 'E:/BaiduNetdiskDownload/bi_data/temp'logging.info("Running Mountainsort4 sorter...")sorting = sis.run_sorter(sorter_name="mountainsort4",recording=recording_saved,output_folder="mountainsort4_output",verbose=False,docker_image=False,freq_min=500,  # 改为 500 Hz（适用于数据集 I）freq_max=5000,num_workers=8,clip_size=50,detect_threshold=3,detect_interval=10)# 提取 spike train 数据unit_ids = sorting.get_unit_ids()spike_data = []for unit in unit_ids:spike_times = sorting.get_unit_spike_train(unit)unit_labels = np.full(len(spike_times), unit)spikes = np.column_stack((spike_times, unit_labels))spike_data.append(spikes)if spike_data:spike_array = np.vstack(spike_data)else:spike_array = np.empty((0, 2))# 保存结果np.save("sorting_results.npy", spike_array)np.savetxt("sorting_results.txt", spike_array, fmt="%d", header="spike_time unit_id")logging.info("SUA spike sorting completed. %d units detected. Results saved to sorting_results.txt", len(unit_ids))return sortingexcept Exception as e:logging.error("SUA sorting failed: %s", str(e))return None

3. MUA（多单元活动）处理

  • 代码中的设置：

    • 使用 500-5000 Hz 带通滤波，负向阈值检测尖峰事件。

  • 图片中的设置：

    • 带通滤波范围为 500-5000 Hz（对于数据集 I）或 250-5000 Hz（对于数据集 II）。

  • 问题：

    • 代码中的带通滤波范围是 500-5000 Hz，与图片中的设置一致。但如果使用的是数据集 II，需要调整为 250-5000 Hz。

  • 修改建议：

    • 根据数据集类型调整带通滤波范围。如果使用的是数据集 II，应将带通滤波范围改为 250-5000 Hz。

  • 修改原因：

    • 不同的数据集可能需要不同的滤波范围来优化 MUA 检测的效果。

修改后的代码：

def detect_mua(recording, n_channels):logging.info("**** MUA Detection ****")recording_spike = spre.bandpass_filter(recording, freq_min=500, freq_max=5000)  # 根据数据集类型调整 freq_min# 如果使用数据集 II，改为 freq_min=250def detect_spikes(trace, threshold):above_threshold = trace < -thresholdspike_indices = np.where(np.diff(above_threshold.astype(int)) == 1)[0] + 1return spike_indicesmua = {}for ch in range(n_channels):ch_trace = recording_spike.get_traces(channel_ids=[ch]).flatten()thr = 5 * np.std(ch_trace)spike_idx = detect_spikes(ch_trace, thr)mua[ch] = spike_idxlogging.debug("Channel %d: detected %d spikes", ch, len(spike_idx))logging.info("MUA detection completed.")return mua

4. ESA（整体尖峰活动包络）处理

无需修改

总结：

  • LFP：应将带通滤波改为低通滤波，截止频率为 300 Hz。

  • SUA：根据数据集类型调整带通滤波范围，数据集 I 为 500-5000 Hz，数据集 II 为 250-5000 Hz。

  • MUA：根据数据集类型调整带通滤波范围，数据集 I 为 500-5000 Hz，数据集 II 为 250-5000 Hz。

  • ESA：无需修改，代码中的设置与文献图片中的建议一致。

问题3 ：学生主要问题是spikeinterface有这么多种算法，滤波器也有很多种类，他们所用的超参数组合更是无法遍历，我们希望能听一下您的经验，一般使用什么样的滤波器和sorting算法组合

馒头前言：

        我们在使用 SpikeInterface 进行 spike sorting 时，选择合适的滤波器和 sorting 算法组合是非常重要的。不同的侵入式数据（如 LFP、SUA、MUA 等）和任务（如单单元活动检测、多单元活动分析等）可能需要不同的处理流程和参数设置。以下是我的一些实操经验和建议，希望帮助你们在处理 Indy 数据集时选择合适的滤波器和 sorting 算法组合。

1. 数据预处理：滤波器的选择

滤波器的选择取决于你要提取的信号类型（如 LFP、SUA、MUA 、ESA）。以下是一些常见的滤波器设置：

（1）LFP（局部场电位）

  • 任务：提取低频信号，通常用于分析局部神经元群体的电活动。

  • 滤波器：低通滤波器（LPF）

  • 参数：

    • 截止频率：300 Hz

    • 滤波器类型：Butterworth 或 Chebyshev

（2）SUA（单单元活动）

  • 任务：提取单个神经元的尖峰信号，通常用于 spike sorting。

  • 滤波器：带通滤波器（BPF）

  • 参数：

    • 截止频率：300 Hz - 5000 Hz（适用于大多数侵入式数据）

    • 滤波器类型：Butterworth 或 Chebyshev

（3）MUA（多单元活动）

  • 任务：提取多个神经元的尖峰信号，通常用于分析电极周围的总体活动。

  • 滤波器：带通滤波器（BPF）

  • 参数：

    • 截止频率：500 Hz - 5000 Hz（适用于大多数侵入式数据）

    • 滤波器类型：Butterworth 或 Chebyshev

4）ESA（整体尖峰活动包络）

  • 任务：提取整体尖峰活动的包络信号，通常用于分析神经活动的总体趋势。

  • 滤波器：高通滤波器（HPF） + 低通滤波器（LPF）

  • 参数：

    • 高通截止频率：300 Hz

    • 低通截止频率：12 Hz

    • 滤波器类型：Butterworth 或 Chebyshev

2. Spike Sorting 算法的选择

SpikeInterface 提供了多种 spike sorting 算法，每种算法适用于不同的数据类型和任务。以下是一些常见的算法及其适用场景：

（1）Mountainsort4

  • 适用场景：适用于大多数侵入式数据，尤其是高质量的 SUA 数据。

  • 优点：

    • 准确性高，适合处理单单元活动。

    • 支持多线程，计算速度较快。

  • 参数建议：

    • freq_min: 300 Hz

    • freq_max: 5000 Hz

    • detect_threshold: 3.0（根据数据噪声水平调整）

    • clip_size: 50（尖峰波形长度）

（2）Kilosort2

  • 适用场景：适用于高通道数的侵入式数据（如 Neuropixels 数据）。

  • 优点：

    • 能够处理高通道数的数据，适合大规模记录。

    • 准确性高，适合复杂的数据集。

  • 参数建议：

    • freq_min: 300 Hz

    • freq_max: 5000 Hz

    • detect_threshold: 3.0（根据数据噪声水平调整）

（3）Spyking Circus

  • 适用场景：适用于中等通道数的侵入式数据，适合处理 SUA 和 MUA。

  • 优点：

    • 支持多通道数据的联合排序，适合处理复杂的尖峰活动。

  • 参数建议：

    • freq_min: 300 Hz

    • freq_max: 5000 Hz

    • detect_threshold: 3.0（根据数据噪声水平调整）

（4）Tridesclous

  • 适用场景：适用于小规模数据集，适合初学者使用。

  • 优点：

    • 参数设置简单，适合快速上手。

  • 参数建议：

    • freq_min: 300 Hz

    • freq_max: 5000 Hz

    • detect_threshold: 3.0（根据数据噪声水平调整）

3. 实操经验与建议

（1）数据质量检查

  • 在运行 spike sorting 之前，先检查数据的质量（如信噪比、噪声水平等）。

  • 可以使用 spikeinterface.toolkit 中的工具（如 compute_noise_levels）来评估数据质量。

（2）参数调优

  • 检测阈值（detect_threshold）：根据数据的噪声水平调整，通常设置为 3.0 - 5.0 倍的标准差。

  • 滤波范围：根据信号类型选择合适的滤波范围（如 SUA 使用 300-5000 Hz，MUA 使用 500-5000 Hz）。

  • 尖峰波形长度（clip_size）：通常设置为 50-100 个样本点，确保能够捕捉完整的尖峰波形。

（3）多算法对比

  • 对于同一份数据，可以尝试多种 spike sorting 算法（如 Mountainsort4、Kilosort2 等），并对比结果。

  • 使用 spikeinterface.comparison 模块进行算法结果的比较和评估。

（4）并行计算

  • 对于大规模数据，可以使用多线程或分布式计算来加速 spike sorting 过程。

  • 在 run_sorter 中设置 num_workers 参数来启用多线程。

问题4 学生：不清楚探针信息对sorting影响有多大，最好是可以指导我怎么自定义代码探针

这个问题和问题1重叠了，但这里我还是说一下【探针对排序的影响】

探针的影响：

  • 几何信息：电极位置直接影响信号的空间插值、通道邻域关系（用于聚类和去噪）。

  • 信号质量：探针布局不匹配会导致信号参考错误或通道间干扰放大。

自定义探针代码示例：

请参照问题1给出的探针代码步骤学习

http://www.xdnf.cn/news/11881.html

相关文章：

  • 详解Jenkins Pipeline 中git 命令的使用方法
  • Kubernetes 网络方案：Flannel 插件全解析
  • 常用的录音芯片型号有哪些?
  • 高并发区块链系统实战：从架构设计到性能优化
  • NFS的基本配置
  • Java中的多态
  • Java SpringBoot 调用大模型 AI 构建智能应用实战指南
  • 在树莓派上添加音频输入设备的几种方法
  • Rust学习（1）
  • 采用 Docker GPU 部署的 Ubuntu 或者 windows 桌面环境
  • Elasticsearch中的刷新（Refresh）和刷新间隔介绍
  • 【Zephyr 系列 7】BLE 数据透传系统设计与实现：双向通信、缓冲区与状态同步全解析
  • c++第6天--运算符重载
  • Linux基础开发工具——yum工具
  • Flutter快速上手，入门教程
  • stm32——UART和USART
  • URL 结构说明+路由（接口）的认识
  • VSCode 工作区配置文件通用模板（CMake + Ninja + MinGW/GCC 编译器 的 C++ 或 Qt 项目）
  • 软件工程：如何在项目中把软件做好
  • Python爬虫：trafilatura 的详细使用（高效的网页正文提取工具）
  • AI自动化任务执行工具OpenManus一键启动整合包
  • k8s热更新-subPath 不支持热更新
  • 网络安全中网络诈骗的攻防博弈
  • label-studio 标注实体识别
  • 第三篇：MCP协议深度实践——从理论到生产级AI智能体系统
  • go的工具库：github.com/expr-lang/expr
  • 杰发科技AC7840——Timer修改重装载值
  • 深入解析Java17核心新特性（增强NullPointerException、强封装 JDK 内部 API、伪随机数生成器增强）
  • LVDS与GMSL和FPD-Link之间的关系
  • CppCon 2015 学习:All Your Tests are Terrible