第三章：Python科学计算栈——神经信号的四大金刚

做神经信号处理这些年，我越来越觉得Python科学计算栈就像一套瑞士军刀。NumPy是刀刃，SciPy是精密工具，Matplotlib是显示器，Pandas是数据管家。今天咱们就把这四样东西掰开揉碎了聊。

3.1 NumPy数组操作：神经信号的数字骨架

说实话，我见过太多新手一上来就用Pandas处理原始神经信号。其实不对。原始信号本质上是多维数组，NumPy才是第一选择。

举个例子，你从采集设备拿到的脑电数据，通常是 通道数 × 时间点 的矩阵。我习惯用 np.loadtxt() 或 np.fromfile() 直接读入，然后做切片、统计、归一化。

import numpy as np

# 模拟64通道、1000Hz采样、10秒的EEG数据
eeg_data = np.random.randn(64, 10000)

# 提取第3通道前5秒的数据
ch3_5s = eeg_data[2, :5000]

# 计算每个通道的均值、方差
channel_mean = np.mean(eeg_data, axis=1)
channel_std = np.std(eeg_data, axis=1)

# 归一化：Z-score标准化
eeg_norm = (eeg_data - channel_mean[:, np.newaxis]) / channel_std[:, np.newaxis]

NumPy的索引操作也很关键。比如你想剔除某个坏通道，或者截取刺激呈现前后的时间段，用切片就能搞定。我个人最常用的是 np.where() 来定位特定事件的时间点。

# 找到信号中超过3个标准差的异常点
threshold = 3 * np.std(eeg_data[0])
bad_indices = np.where(np.abs(eeg_data[0]) > threshold)[0]
print(f"发现 {len(bad_indices)} 个异常点")

3.2 SciPy信号处理基础：滤波与变换的利器

SciPy的 signal 模块，说白了就是信号处理的百宝箱。滤波、频谱分析、小波变换，全在这。

我记得刚入行时，导师让我设计一个带通滤波器，把0.5-40Hz的脑电信号留下来。那时候我还不懂滤波器阶数、截止频率这些概念，结果滤波后的信号相位全乱了。后来才明白，要用 filtfilt 做零相位滤波。

from scipy import signal

# 设计一个4阶巴特沃斯带通滤波器，0.5-40Hz
fs = 1000  # 采样率
lowcut = 0.5
highcut = 40.0
nyquist = 0.5 * fs
low = lowcut / nyquist
high = highcut / nyquist
b, a = signal.butter(4, [low, high], btype='band')

# 零相位滤波（重要！）
eeg_filtered = signal.filtfilt(b, a, eeg_data[0])

频谱分析也是家常便饭。计算功率谱密度，看看alpha波（8-12Hz）有没有增强，或者beta波（13-30Hz）有没有异常，用 signal.welch 最方便。

# 计算功率谱密度
freqs, psd = signal.welch(eeg_filtered, fs, nperseg=1024)

# 提取alpha频段的平均功率
alpha_mask = (freqs >= 8) & (freqs <= 12)
alpha_power = np.mean(psd[alpha_mask])
print(f"Alpha波功率：{alpha_power:.3f} µV²/Hz")

3.3 Matplotlib可视化：让神经信号说话

做神经科学实验，可视化不是装饰，是诊断工具。我每次处理新数据，第一件事就是画图看看——信号有没有饱和？有没有工频干扰？基线漂移严重不严重？

Matplotlib的 pyplot 接口，配合 subplots 可以轻松画出多通道脑电波形图。我个人习惯把时间轴统一，通道按顺序排列，这样一眼就能看出哪个通道有问题。

import matplotlib.pyplot as plt

# 绘制前8个通道的波形
fig, axes = plt.subplots(8, 1, figsize=(12, 8), sharex=True)
time = np.arange(10000) / fs

for i, ax in enumerate(axes):
    ax.plot(time[:2000], eeg_data[i, :2000], lw=0.5)
    ax.set_ylabel(f'Ch{i+1}')
    ax.set_ylim(-5, 5)

axes[-1].set_xlabel('时间 (秒)')
plt.tight_layout()
plt.show()

除了波形图，频谱图、时频图、拓扑图也经常用到。比如用 plt.specgram 画时频图，能看出不同频段随时间的变化，对分析睡眠纺锤波特别有用。

3.4 Pandas时间序列数据处理：时间索引的艺术

嗯，这里要注意。Pandas虽然强大，但我不建议用它存原始神经信号——太慢了。它的真正战场是事件日志、行为数据、元数据这些带时间戳的东西。

比如你有一个实验，记录了每个试次的刺激呈现时间、反应时间、正确与否。用Pandas的 DataFrame 来管理，配合 datetime 索引，做筛选、对齐、重采样都非常方便。

import pandas as pd

# 创建事件日志
events = pd.DataFrame({
    'trial': range(1, 101),
    'stim_time': pd.date_range('2024-01-01 10:00:00', periods=100, freq='2s'),
    'reaction_time': np.random.uniform(0.2, 0.8, 100),
    'correct': np.random.choice([True, False], 100, p=[0.8, 0.2])
})

# 设置时间索引
events.set_index('stim_time', inplace=True)

# 筛选正确试次中反应时间小于0.5秒的
fast_correct = events[(events['correct']) & (events['reaction_time'] < 0.5)]
print(f"快速正确试次数：{len(fast_correct)}")

Pandas的 resample 方法也很实用。比如行为数据是100Hz采的，神经信号是1000Hz，你想把行为数据降采样到和神经信号对齐，一行代码就搞定。

# 假设行为数据是100Hz，重采样到1000Hz（线性插值）
behavior_1khz = behavior_100hz.resample('1ms').interpolate()

---