第三讲：时域与频域特征提取——从神经信号里“挖”出有用信息

大家好，欢迎来到实战第三讲。

上一讲我们聊了怎么把原始脑电数据读进来、做预处理。但说实话，预处理只是“打扫干净屋子”。真正要干活，还得从信号里提取出有意义的特征。

这一讲，我们就来聊聊怎么从神经信号中提取时域和频域特征。我会用 Scipy 和 MNE 这两个库，带大家一步步实现。

1. 时域特征：信号最直观的“长相”

时域特征，说白了就是直接看信号波形本身。你想想看，一段脑电信号，它高高低低、起起伏伏，这里面其实藏着很多信息。

均值（Mean）

均值反映的是信号的直流分量。在脑电里，它通常代表电极与头皮之间的接触电位。我个人习惯在预处理阶段就把均值归零（去均值），这样后续分析更干净。

import numpy as np

# 假设 data 是形状为 (channels, timepoints) 的数组
mean_val = np.mean(data, axis=1)
print(f"各通道均值: {mean_val}")

方差（Variance）

方差衡量信号的波动程度。脑电信号活跃的时候，方差会变大；安静或者睡眠时，方差会变小。这个特征在检测事件相关电位（ERP）时特别有用。

var_val = np.var(data, axis=1)
print(f"各通道方差: {var_val}")

过零率（Zero Crossing Rate, ZCR）

过零率就是信号穿过零点的次数。高频信号过零率大，低频信号过零率小。嗯，这里要注意：如果信号已经去均值了，过零率能很好地反映频率成分的分布。

from scipy.signal import zero_crossings

def compute_zcr(signal):
    # zero_crossings 返回布尔数组，True 表示过零点
    crossings = zero_crossings(signal)
    return np.sum(crossings) / len(signal)

zcr_val = compute_zcr(data[0, :])  # 计算第一个通道的过零率
print(f"过零率: {zcr_val:.4f}")

2. 频域特征：换个角度看信号

时域特征虽然直观，但有些信息在时域里看不清楚。比如，alpha波（8-12Hz）在时域里可能只是一堆小波浪，但转到频域，它就是一个明显的峰值。

功率谱密度（Power Spectral Density, PSD）

PSD 告诉我们信号在不同频率上的能量分布。计算 PSD 最常用的方法是 Welch 法——说白了就是把信号分段、加窗、做 FFT，然后平均。

from scipy.signal import welch
import matplotlib.pyplot as plt

fs = 250  # 采样率，假设250Hz
freqs, psd = welch(data[0, :], fs=fs, nperseg=fs*2)  # 每段2秒

plt.semilogy(freqs, psd)
plt.xlabel('频率 (Hz)')
plt.ylabel('功率谱密度 (dB)')
plt.title('单通道 PSD')
plt.show()

用 MNE 做这件事更简单，它封装好了很多细节：

import mne

# 假设 raw 是 MNE 的 Raw 对象
raw.plot_psd(fmin=0.5, fmax=45, tmin=0, tmax=60)

频带能量（Band Power）

频带能量就是把特定频段内的 PSD 积分（求和）。常用的频带包括：

频带	频率范围	生理意义
Delta	0.5 - 4 Hz	深度睡眠
Theta	4 - 8 Hz	冥想、困倦
Alpha	8 - 12 Hz	放松、闭眼
Beta	12 - 30 Hz	专注、活跃
Gamma	30 - 45 Hz	高级认知

def band_power(psd, freqs, band):
    idx = np.where((freqs >= band[0]) & (freqs <= band[1]))[0]
    return np.trapz(psd[idx], freqs[idx])  # 用梯形法积分

alpha_power = band_power(psd, freqs, [8, 12])
beta_power = band_power(psd, freqs, [12, 30])
print(f"Alpha 能量: {alpha_power:.3f}, Beta 能量: {beta_power:.3f}")

3. 完整示例：从 MNE Raw 对象提取特征

最后，我们串起来走一遍完整流程。假设你已经有了一个预处理好的 Raw 对象：

import mne
import numpy as np
from scipy.signal import welch

# 加载数据（示例用MNE内置数据）
sample_data_folder = mne.datasets.sample.data_path()
raw = mne.io.read_raw_fif(sample_data_folder + '/MEG/sample/sample_audvis_raw.fif', preload=True)
raw.pick_types(meg=False, eeg=True)  # 只保留EEG通道
raw.filter(0.5, 45)  # 带通滤波

# 提取数据
data = raw.get_data()
fs = raw.info['sfreq']
n_channels = data.shape[0]

# 时域特征
mean_vals = np.mean(data, axis=1)
var_vals = np.var(data, axis=1)
zcr_vals = np.array([compute_zcr(data[ch, :]) for ch in range(n_channels)])

# 频域特征
freqs, psd = welch(data, fs=fs, nperseg=int(fs*2), axis=1)
alpha_power = np.array([band_power(psd[ch, :], freqs, [8, 12]) for ch in range(n_channels)])
beta_power = np.array([band_power(psd[ch, :], freqs, [12, 30]) for ch in range(n_channels)])

# 整理成特征矩阵
features = np.column_stack([mean_vals, var_vals, zcr_vals, alpha_power, beta_power])
print(f"特征矩阵形状: {features.shape}")  # (通道数, 特征数)

好了，这一讲的内容就到这里。时域和频域特征是最基础、也是最常用的两类特征。我个人觉得，初学者先把这两个玩熟，后面再接触时频分析、非线性特征就会轻松很多。

记住：特征提取没有银弹。不同任务、不同数据，适合的特征组合也不一样。多试、多对比，慢慢就会有自己的心得。

---

公众号:蓝海资料掘金营，微信deep3321