滑动窗口相关

# 生成两组两通道信号，各 600 个时间点：A（平稳）全程真实相关 0.3；
# B（切换）前 300 点真实相关 0.7、后 300 点真实相关 0。sliding_corr 返回每个窗口的相关系数。
# 输出解读：B 在中点附近明显下降；A 的真实相关恒定，窗口相关却同样明显起伏——
# 有限窗口的采样变异本身就会产生“波动”，起伏不能直接说明连接在变化。
import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.stride_tricks import sliding_window_view

def ar1_smooth(e, phi=0.8):
    """一阶自回归平滑，模仿 BOLD 的自相关（模拟数据，仅用于演示方法）。"""
    out = np.zeros_like(e)
    for t in range(1, e.shape[-1]):
        out[..., t] = phi * out[..., t - 1] + e[..., t]
    return out

def correlated_pair(rho, T, rng):
    z = rng.standard_normal((2, T))
    return ar1_smooth(np.vstack([z[0], rho * z[0] + np.sqrt(1 - rho ** 2) * z[1]]))

def sliding_corr(x, y, w, step=1):
    """窗长 w、步长 step 的滑动窗口皮尔逊相关（矩形窗，向量化实现）。"""
    xw = sliding_window_view(x, w)[::step]
    yw = sliding_window_view(y, w)[::step]
    xc = xw - xw.mean(axis=1, keepdims=True)
    yc = yw - yw.mean(axis=1, keepdims=True)
    return (xc * yc).sum(1) / np.sqrt((xc ** 2).sum(1) * (yc ** 2).sum(1))

rng = np.random.default_rng(1)
A = correlated_pair(0.3, 600, rng)
B = np.hstack([correlated_pair(0.7, 300, rng), correlated_pair(0.0, 300, rng)])

w = 40
rA = sliding_corr(A[0], A[1], w)
rB = sliding_corr(B[0], B[1], w)
print(f"A（平稳 0.3）：窗口相关范围 {rA.min():.2f} ~ {rA.max():.2f}，标准差 {rA.std():.2f}")
print(f"B（切换）：前 200 个窗口均值 {rB[:200].mean():.2f}，后 200 个窗口均值 {rB[-200:].mean():.2f}")

plt.plot(rA, label="A: true r = 0.3 (stationary)")
plt.plot(rB, label="B: true r 0.7 -> 0")
plt.axhline(0.3, ls="--", c="gray")
plt.xlabel("window start")
plt.ylabel("windowed correlation")
plt.legend()
plt.title(f"Sliding-window correlation (w = {w}, simulated data)")
plt.show()