实践工作
I. 概述
本工作的目标是：(i) 从电生理信号的噪声混合物中提取感兴趣的信号，在耐药性癫痫的背景下；(ii) 使用模拟的类癫痫脑内脑电图（iEEG）信号研究脑连接性。假设使用 N_c 个传感器在观察时间段 T 内观察到三个电生理信号的线性瞬时混合物，这些传感器放置在癫痫患者的头皮上。获取的数据随后以时空矩阵的形式给出，用 \mathbf{X} (N_c \times T) 表示，使得
其中 \mathbf{A}(N_c \times P) 是混合矩阵，建模从源子空间向观察空间的传播介质，\mathbf{S}(P \times T) 是一个矩阵，其第 p 行表示第 p 个源的时间轮廓。请注意，本研究中考虑了仪器无噪声观测模型。在整个工作中，考虑两个传感器。
在第II和第III节中，涉及三个生物物理信号：
- s_1(n)：棘间波，代表大脑中一个神经群体的癫痫活动
- s_2(n)：眼动，代表伪影信号
- s_3(n)：肌肉活动。
请注意，在本工作中，眼动和肌肉活动都被视为干扰感兴趣信号（即棘间波）的伪影。如前所述，本实践工作的一个目标是将感兴趣的信号与伪影信号（眼动和肌肉活动）分离。为了处理不同的伪影去除技术，首先将使用自适应滤波方法从时空混合物中去除眼动活动。其次，将使用独立成分分析（ICA）技术去除剩余的肌肉活动。最后，基于模拟的去噪类癫痫脑内脑电图（iEEG）信号，将讨论神经群体之间的连接性问题。
II. 第一步：使用带参考通道的自适应噪声降低去除眼动伪影
给定两个传感器，本节的目标是从传感器输出中去除眼动活动。为此，需要考虑两个伪影去除系统（每个通道一个）。为了清晰起见，这里仅描述一个对应于一个通道（即一个传感器）的伪影去除系统，请记住，在考虑另一个通道（即另一个传感器）时需要使用相同的步骤。
我们想要设计一个噪声降低系统，其中两个观测值 y_1(n) 和 y_2(n) 可用，使得
其中 s(n) 是感兴趣的信号，b_1(n) 代表要从观测信号 y_1(n) 中去除的噪声/伪影。参考信号 b_2(n) 是必须与目标噪声/伪影 b_1(n) 相关的信号。因此，b_2(n) 可以被视为 b_1(n) 的滤波版本。滤波器 h 是自适应的，我们希望最小化 E[(\hat{s})^2]，其中 \hat{s} = s + b_1 - h(b_2)。这等价于最小化均方误差 E[(b_1 - h(b_2))^2]。在实践中，滤波器 h 是一个具有 M 个系数（M = 20）的横向滤波器，其目标是通过计算噪声/伪影 b_1 的最佳估计 \hat{b}_1 = h(b_2)，从可用的参考信号 b_2 中提供信号 s 的最佳估计 \hat{s}。描述LMS算法的主要方程如下：
图1. 噪声降低系统
1. 在测试真实癫痫信号之前，我们建议在以下全模拟噪声数据上测试自适应算法：信号 s(n) 是振幅 A = 2 的正弦信号，在100 Hz振荡并以1 kHz采样，b_2(n) 是零均值单位方差的白高斯噪声，b_1(n) 是使用传递函数等于 1/(1-\alpha z^{-1}) 的滤波器从 b_2(n) 生成的。编写相应的程序并针对 \alpha = 0.9，M = 20 和不同的 \mu 值（例如：\mu = 0.01, \mu = 0.0005）进行测试。
2. 对于 \mu = 0.01, \mu = 0.005, \mu = 0.001 和 \mu = 0.0005，绘制从时刻1到时刻100以及从时刻5000到时刻5120的输出信号。
3. 对于 \mu = 0.01, \mu = 0.005, \mu = 0.001 和 \mu = 0.0005，绘制所有时刻的项 (s - \hat{s})^2。
4. 讨论步长 \mu 在稳态和收敛速度上的影响。
5. 通过对最后500个点取平均来计算稳态误差。
6. 在这个问题中，感兴趣的信号是癫痫棘间波序列。
a) 使用"signals.mat"加载三个信号，使用Matlab命令"load"，其中信号"s1"是棘间波，"s2"表示眼动，信号"s3"代表肌肉活动。
b) 按以下方式构建时空混合矩阵 \mathbf{X}(2 \times T)：
其中 \mathbf{S}(i,:) 是 \mathbf{S} = [\mathbf{s}_1^T, \mathbf{s}_2^T, \mathbf{s}_3^T]^T 的第 i 行，\mathbf{A} = [\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3] = \begin{bmatrix} 1 & 0.3 & 0.5 \\ -0.7 & 0.6 & 0.5 \end{bmatrix} 且 \sigma_1, \sigma_2 是一些标量，用于调整信噪比（SNR）。为简单起见，我们在这里考虑 \sigma_1 = \sigma_2 = \sigma，其中 SNR = 20\log(\sigma_s/\sigma)，\sigma_s 是感兴趣信号的标准差。我们设置 \sigma_s = 1 且 SNR = 5dB。
c) 绘制信号 s_1(n)，s_2(n)，s_3(n)，x_1(n)（即 \mathbf{X}(1,:)）和 x_2(n)（即 \mathbf{X}(2,:)）。
7. 根据上面定义的混合模型，混合矩阵 \mathbf{X} 的每一行，这里用 x_i(n) i \in \{1,2\} 表示，代表在头皮表面记录的棘间癫痫波与肌肉活动的噪声混合物，如前所述。在这一部分中，对于参考噪声（眼动活动），b_2(n) 是与目标眼动活动 b_1(n)（例如纯记录的眼动活动）相关的信号，并在**"Ref_Ocl.mat"** Matlab文件中提供。编写程序以实现自适应LMS算法，从棘间波和肌肉伪影混合物中去除眼动活动的贡献。执行上述步骤3和4。
III. 第二步：使用ICA去除肌肉运动
在使用第一步去除眼动活动后，我们现在对从结果混合物 \hat{\mathbf{X}}_{(2\times N)} = [\hat{x}_1^T, \hat{x}_2^T]^T 中去除剩余伪影（即肌肉运动）感兴趣。由于 \hat{X} 现在是来自两个生理学上不同源（即iEEG棘间波和肌肉运动）的两个源（N：样本数）的混合物，这些源在统计上是独立的，ICA技术非常适合处理这个盲源分离（BSS）问题。回顾一下，源的统计独立性是ICA [1] 技术用于解决BSS问题所依赖的主要假设。众所周知的SOBI（二阶盲辨识）方法[2]将用于从其伴随伪影（即肌肉运动）中提取感兴趣的信号（即iEEG信号）。基于预白化数据，SOBI方法通过联合对角化从预白化数据构建并在不同时间滞后处获取的一组协方差矩阵来解决BSS问题。因此，我们将编辑一个Matlab脚本来：
1. 使用PCA预白化可用观测值 \hat{\mathbf{X}}_{(2\times N)}：
a. 计算零均值观测值 \hat{X}_c 的协方差矩阵
b. 计算 \mathbf{R}_{\hat{x}} 的特征值分解（EVD）：\mathbf{R}_{\hat{x}} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^T = \tilde{\mathbf{A}}\mathbf{R}_{\tilde{s}}\tilde{\mathbf{A}}^{\dagger}，其中 \mathbf{R}_{\tilde{s}} 是新源向量 \tilde{\mathbf{s}}(n) = [s_1(n), s_2(n)]^T 的协方差矩阵。将 \mathbf{R}_{\tilde{s}} 的特征值按降序排序，并相应地重新排列 \mathbf{U} 的列以得到 \tilde{\mathbf{U}}。
c. 设置预白化算子 \mathbf{W} = \mathbf{Q}^{-1}\mathbf{R}_{\tilde{s}}^{-1/2}\tilde{\mathbf{A}}^{-1} = \mathbf{\Lambda}^{-1/2}\tilde{\mathbf{U}}^{\dagger}，其中 \dagger 是伪逆算子（Matlab命令"pinv"用于执行此类操作），\mathbf{Q} 是酉矩阵。
d. 将原始数据矩阵 \hat{X}_c 投影到新子空间上：\mathbf{Z} = \mathbf{W}\hat{X}_c = \mathbf{Q}\mathbf{S}，其中 \mathbf{Z} 表示预白化数据矩阵。
2. 肌肉活动去除和iEEG提取使用SOBI方法：
a. 估计一组基于样本二阶（SO）矩的空间协方差矩阵：\mathbf{R}_{\mathbf{Z}}(\tau_k) = \frac{1}{N}\sum_{n=1}^{N}\mathbf{z}[n]\mathbf{z}[n+\tau_k]^T，\tau_k = 1:10
b. 通过联合对角化计算的协方差矩阵集 \mathbf{R}_{\mathbf{Z}}(\tau_k) 来估计新的混合矩阵 \mathbf{Q}，使用提供的"jad.m"方法。
c. 估计包含iEEG活动和眼动的源向量，通过：\hat{\mathbf{S}} = \hat{\mathbf{Q}}^{\dagger}\mathbf{Z}，其中 \hat{\mathbf{Q}} 是 \mathbf{Q} 的估计。
d. 绘制提取的两个信号（即iEEG和肌肉运动）。
e. 使用Matlab命令"corr"，计算真实iEEG信号 s_1(n) 与提取的信号 \hat{s}_1(n) 之间的相关系数。验证您的结果。
IV. 第三步：脑连接性
本研究的主要目标是基于各自的iEEG信号来量化两个神经集合之间的功能连接性。这种连接性使用相干函数来估计。
初步问题
在第一步中，您必须验证，给定两个在二阶平稳的中心化随机过程 u(n) 和 v(n)，我们有 0 \leq |\rho_{uv}(f)|^2 \leq 1，其中 \rho_{uv}(f) 是相干函数，定义为 \rho_{uv}(f) = \gamma_{uv}(f)/\sqrt{\gamma_{uu}(f)\gamma_{vv}(f)}（\gamma_{uv}(f)，\gamma_{uu}(f)，\gamma_{vv}(f) 分别是两个信号的互功率谱密度（psd）和功率谱psd）。为此，让我们设计 \lambda 为复数任意数，并计算过程 z(n) 的psd，定义为：z(n) = u(n) + \lambda v(n)。然后选择 \lambda = r\gamma_{uv}(f)/|\gamma_{uv}(f)| 并给出证明。
实验
现在，让我们考虑对应于以下向量线性自回归模型的两个iEEG信号（见图2）：
(1)
其中 w_j，j = 1,2 是独立白噪声 w_j 的实现，具有零均值和单位方差。在方程(1)中，y_1(n) 和 y_2(n) 可以假定为iEEG信号。
[图2：显示节点1指向节点2的连接图]
图2
a. 在第一步中，我们仅考虑 y_1(n) 作为方程(1)中给出的AR模型，y_2(n) 是 y_1(n) 的延迟版本，使得 y_2(n) = 0.5y_1(n-1)。计算并绘制两个信号在4个1024点块上的频谱，每个块有50%的重叠（使用"pwelch"命令）。计算并绘制这两个信号之间在不同块上的幅度平方相干（MSC），并得出结论（使用"mscohere"命令）。
b. 在第二步中，我们考虑方程(1)中给出的两个信号 y_1(n) 和 y_2(n)。计算两个信号在100个1024点块上的频谱，每个块有50%的重叠，以及MSC。绘制在这100个块上计算的平均MSC并评论您的结果。
c. 然后，我们引入从生理学模型获得的两个信号（在"ScenarioK1500second1000a.mat" Matlab文件中提供并使用LFP s），我们估计这些信号之间随时间变化的MSC。显示信号的频谱以及MSC。
V. 参考文献
[1] A. Hyvarinen, J. Karhunen and E. Oja, "Independent component analysis", Adaptive and learning systems for signal processing, communications and control, John Wiley, New York, Chichester, 2001.
[2] A. Belouchrani, K. Abed-Mariam, J. F. Cardoso and E. Moulines, "A Blind Source Separation Technique Using Second-Order Statistics", IEEE Trans. On Sig. Proc. Vol. 45, No. 2, 1997.