生物医学网络科学Biomedical Network Science
Example04 / 11

相位同步与主导特征向量（LEiDA 思路）

所属专题
动态功能连接与多层脑网络
下载
Notebook · .py

对窄带信号做希尔伯特变换得到瞬时相位 θ_i(t)，逐时间点构造相位一致性矩阵 P_ij(t) = cos(θ_i − θ_j)。 P(t) 的秩不超过 2，取其主导特征向量 v1(t) 即可概括该时刻的相位结构，再对 v1(t) 做 k-means 得到状态。 模拟：8 个“脑区”，前半段 0–4 与 5–7 两组反相，后半段全部同相（模拟数据，仅用于演示方法）。 输出解读：主导特征向量的符号直接给出“谁和谁同相”；聚类标签在中点附近切换。 注意：瞬时相位只对窄带信号有意义，真实 BOLD 需先带通滤波（如 0.01–0.08 Hz 内再取较窄频带）。

Output输出图

相位同步与主导特征向量（LEiDA 思路）：输出图 1

Output文字输出

前半段主导特征向量（取中位数）符号： [-1 -1 -1 -1 -1  1  1  1]
后半段主导特征向量（取中位数）符号： [-1 -1 -1 -1 -1 -1 -1 -1]
两段各自主导标签的占比（平均）= 1.00；前后两段主导标签是否不同： True

Code代码

# 对窄带信号做希尔伯特变换得到瞬时相位 θ_i(t)，逐时间点构造相位一致性矩阵 P_ij(t) = cos(θ_i − θ_j)。
# P(t) 的秩不超过 2，取其主导特征向量 v1(t) 即可概括该时刻的相位结构，再对 v1(t) 做 k-means 得到状态。
# 模拟：8 个“脑区”，前半段 0–4 与 5–7 两组反相，后半段全部同相（模拟数据，仅用于演示方法）。
# 输出解读：主导特征向量的符号直接给出“谁和谁同相”；聚类标签在中点附近切换。
# 注意：瞬时相位只对窄带信号有意义，真实 BOLD 需先带通滤波（如 0.01–0.08 Hz 内再取较窄频带）。
import numpy as np
from scipy.signal import hilbert, butter, filtfilt
from scipy.cluster.vq import kmeans2

rng = np.random.default_rng(0)
N, T, fs = 8, 1200, 1.0                      # 1 个时间单位 = 1 个采样点
t = np.arange(T)
f0 = 0.05                                    # 共同振荡频率（每采样点 0.05 周期）
base = 2 * np.pi * f0 * t + 0.3 * np.cumsum(rng.standard_normal(T)) * 0.05
offset = np.zeros((N, T))
offset[5:, : T // 2] = np.pi                 # 前半段：5–7 号与 0–4 号反相
X = np.cos(base + offset) + 0.4 * rng.standard_normal((N, T))

b, a = butter(3, [0.03, 0.07], btype="band", fs=fs)
Xf = filtfilt(b, a, X, axis=1)
theta = np.angle(hilbert(Xf, axis=1))       # (N, T) 瞬时相位

def leading_eigvec(th):
    P = np.cos(th[:, None] - th[None, :])
    w, V = np.linalg.eigh(P)
    v = V[:, -1]
    return v if np.sum(v < 0) >= np.sum(v > 0) else -v   # 统一符号：使多数分量为负（v 与 −v 等价）

edge = 50                                    # 去掉两端受滤波边缘效应影响的点
V1 = np.array([leading_eigvec(theta[:, k]) for k in range(edge, T - edge)])
_, lab = kmeans2(V1, 2, seed=1, minit="++")
half = T // 2 - edge
print("前半段主导特征向量（取中位数）符号：", np.sign(np.median(V1[:half], axis=0)).astype(int))
print("后半段主导特征向量（取中位数）符号：", np.sign(np.median(V1[half:], axis=0)).astype(int))
acc = np.mean(lab[:half] == np.bincount(lab[:half]).argmax()) * 0.5 + \
      np.mean(lab[half:] == np.bincount(lab[half:]).argmax()) * 0.5
print(f"两段各自主导标签的占比（平均）= {acc:.2f}；前后两段主导标签是否不同：",
      np.bincount(lab[:half]).argmax() != np.bincount(lab[half:]).argmax())

import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(10, 4.2), sharex=True, gridspec_kw={"height_ratios": [3, 1]},
                         layout="constrained")
tt = np.arange(edge, T - edge)
vm = np.abs(V1).max()
im = axes[0].imshow(V1.T, aspect="auto", cmap="RdBu_r", vmin=-vm, vmax=vm, interpolation="nearest",
                    extent=[tt[0], tt[-1], N - 0.5, -0.5])
axes[0].set(ylabel="Region", title="Leading eigenvector of phase coherence: two anti-phase groups, then one")
fig.colorbar(im, ax=axes, shrink=0.8, label="v1 component")
axes[1].plot(tt, lab, color="k", lw=1)
axes[1].set(xlabel="Time (samples)", ylabel="State", yticks=[0, 1])
plt.show()

Run在浏览器中运行

首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。

相位同步与主导特征向量（LEiDA 思路）