医学信号处理 TP：相干函数性质证明与脑连接性实验

下面给出针对所需实验步骤的详细解答，包括数学推导、从一个公式到下一个公式的推导解释以及相应的示例 MATLAB 代码。请注意，在所有公式中，我会使用

首先，整理一下题目要求与目标：

1. 在完全模拟的数据上测试自适应滤波的噪声消减算法（LMS算法）。
2. 观察不同步长\mu时滤波输出的波形以及误差\bigl(s-\hat{s}\bigr)^2。
3. 讨论步长\mu对稳态以及收敛速度的影响，并计算稳态时的误差。
4. 利用提供的真实信号文件（或模拟的癫痫尖波信号等），建立混合模型并使用自适应滤波去除眼部活动。

以下将分成两大部分进行说明：

• 第一部分：在模拟数据上进行验证（对应题目中的1—5小问）。
• 第二部分：在真实（或更复杂）的癫痫尖波 + 眼部活动 + 肌肉活动信号上进行混合及自适应滤波（对应题目中的6—7小问）。

\textbf{第一部分：在模拟数据上测试自适应滤波算法}

我们已有如下假设和数据生成方式：

• 目标信号：
  s(n) = A\sin(2\pi f \cdot \frac{n}{F_s})
  • 其中 A = 2, f = 100 \,\text{Hz}, 采样率 F_s = 1 \,\text{kHz}。
• 参考噪声：
  b_2(n) \sim \mathcal{N}(0,1) \quad (\text{白高斯噪声})
• 要去除的噪声：
  b_1(n)
  
  通过一阶滤波器
  \frac{1}{1 - \alpha z^{-1}}
  作用于 b_2(n) 得到，其中 \alpha=0.9。

总体观测模型（对于单通道而言）为：

在自适应滤波系统中，假设所用的横向滤波器长度为 M = 20。定义：

其中

LMS算法的核心迭代公式是：

在实现中，需要做以下几个主要步骤：

1. 生成模拟信号 s(n)。
2. 生成白高斯噪声 b_2(n)。
3. 用滤波器（传递函数 1/(1-\alpha z^{-1})）得到 b_1(n)。
4. 分别令 \mu 取不同的值（题目中给了 \mu=0.01,0.005,0.001,0.0005，也可测试更值），记录并分析收敛过程与误差。

下面给出一个简要的 MATLAB 参考实现示例（仅作示例，变量命名请根据你自己的习惯替换）：

示例代码（仅文字说明，无代码块分隔）：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 参数设置
fs = 1000;           % 采样率
N = 6000;            % 信号长度(足够长，保证收敛过程可观测)
n = 0:N-1;           % 时间序列

A = 2;               % 信号幅度
f0 = 100;            % 信号频率
alpha = 0.9;         % 滤波器系数
M = 20;              % 自适应滤波器阶数
muList = [0.01, 0.005, 0.001, 0.0005];  % 不同步长

% 1) 生成目标信号 s(n)
s = A * sin(2*pi*f0*n/fs);

% 2) 生成参考噪声 b2(n) ~ N(0,1)
b2 = randn(1, N);

% 3) 通过一阶滤波器得到 b1(n)
%    传递函数 = 1 / (1 - alpha z^-1)
%    可用滤波器系数实现：分子 = [1], 分母 = [1, -alpha]
b1 = filter([1],[1, -alpha], b2);

% 4) 构造观测信号
y1 = s + b1;  % 通道1含有目标信号和噪声
y2 = b2;      % 通道2即参考噪声

% 开始测试不同 mu
for mu = muList
% 初始化滤波器系数
h = zeros(M,1);  
% 为了记录输出和误差, 建立缓存
s_hat = zeros(1,N);
e = zeros(1,N);  % 这里 e(n) = s(n) - s_hat(n) (或 b1(n) - h^T b2(n))

% 自适应滤波迭代
for k = M:N-1
    % 取参考噪声段 b2(n) ... b2(n-M+1)
    b2_segment = flipud( b2(k-M+1:k)' );
    
    % 当前输出估计 y_hat 对噪声 b1(n) 的估计
    noise_est = h' * b2_segment;
    
    % 估计的 s_hat(n)
    s_hat(k) = y1(k) - noise_est;  % = s(n) + b1(n) - h^T b2(n)
    
    % 计算误差 e(k) = s(k) - s_hat(k)
    e(k) = s(k) - s_hat(k);
    
    % 梯度更新 h(n+1)
    h = h + mu * s_hat(k) * b2_segment;
end

% 画图或存储数据
% 比如, 题目要求 1:100 和 5000:5120 范围内的输出波形
figure;
subplot(211); plot(1:100, s_hat(1:100)); 
title(['Estimated s(n), \mu = ', num2str(mu), ', 1-100']);
subplot(212); plot(5000:5120, s_hat(5000:5120));
title(['Estimated s(n), \mu = ', num2str(mu), ', 5000-5120']);

% 误差 (s - s_hat)^2 全程波形
figure;
plot(e.^2); 
title(['(s - s\_hat)^2 over time, \mu = ', num2str(mu)]);

% 计算后 500 点的平均误差
steady_err = mean(e(N-499:N).^2);
disp(['\mu = ', num2str(mu), ', Steady-state error = ', num2str(steady_err)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

在上面的示例中，我们主要完成了题目 1—5 的要求：

• (1) 编写程序并测试：见上方示例代码，给出了完整的程序框架；可根据题意调参。
• (2) 绘制输出信号在时间点1到100以及5000到5120之间的波形：在示例中通过 plot(s_hat(1:100)) 与 plot(s_hat(5000:5120)) 完成。
• (3) 针对所有时间点绘制 (s - \hat{s})^2：在示例中用 plot(e.^2) 完成。
• (4) 讨论步长\mu的影响：
  • 若\mu过大，则收敛速度快，但稳态时误差可能更大，甚至不收敛；
  • 若\mu过小，则收敛速度慢，但稳态误差可能更低。
    在实际仿真中，可以看到随着\mu减小，误差曲线收敛得更平滑，但需要更多迭代才能达到稳态。
• (5) 通过对最后500个采样点的平均值计算稳态时的误差：在示例中 steady_err = mean(e(N-499:N).^2) 实现了这一功能。

\textbf{第二部分：在真实或更复杂信号上进行去眼部伪差}

接下来，对题目第6—7小问进行回答。

(6) 这里 s(n) 可是一串真实的癫痫发作间期尖波信号。根据题目描述，我们将通过 signals.mat 文件来加载下列三个信号：

• "s1"：间歇性尖峰 (interictal spike)
• "s2"：眼部运动 (ocular movement)
• "s3"：肌肉活动 (muscle activity)

我们将它们组合成一个矩阵：

(a) 加载这三个信号：

load('signals.mat');  % 得到 s1, s2, s3 三个向量

(b) 构建时空混合矩阵 X(2\times T)，其公式如下：

其中：

• \mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3 为列向量，且

  A = \begin{bmatrix}1 & 0.3 & 0.5 \\ -0.7 & 0.6 & 0.5\end{bmatrix}

  表示 A = [a_1,\,a_2,\,a_3]。

• \sigma_1=\sigma_2=\sigma，并且根据

  SNR = 20 \log_{10}\left(\frac{\sigma_s}{\sigma}\right), \quad \sigma_s = 1, \quad SNR = 5\, \text{dB}.

  可解得

  \sigma = \sigma_s \cdot 10^{-\frac{SNR}{20}} = 1 \cdot 10^{-\frac{5}{20}} = 10^{-0.25} \approx 0.5623.

因而：

注意，\|\mathbf{a}_i S(i,:)\| 表示向量 \mathbf{a}_i S(i,:) 的范数，作归一化处理以控制每个信号的贡献。

在 MATLAB 中实现时，可按以下示例思路（文字说明）：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 构建 A 矩阵
A = [1 0.3 0.5; -0.7 0.6 0.5];

% 将 s1, s2, s3 合并
S = [s1; s2; s3];  % 3 x T

% 计算 sigma
sigma_s = 1;
SNR = 5;  % dB
sigma = sigma_s * 10^(-SNR/20);

% 分别计算 a1*S(1,:), a2*S(2,:), a3*S(3,:)
a1 = A(:,1); 
a2 = A(:,2); 
a3 = A(:,3);

temp1 = a1 * S(1,:);  % 2 x T
temp2 = a2 * S(2,:);  % 2 x T
temp3 = a3 * S(3,:);  % 2 x T

% 归一化
temp1 = temp1 / norm(temp1(:));  
temp2 = temp2 / norm(temp2(:));
temp3 = temp3 / norm(temp3(:));

% 按照公式叠加
X = temp1 + sigma*temp2 + sigma*temp3;

% X(1,:) 和 X(2,:) 即 x1(n) 和 x2(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

(c) 绘制 s_1(n)、s_2(n)、s_3(n)，以及混合后的 x_1(n) 和 x_2(n)：

• 直接用 plot(s1), plot(s2), plot(s3), plot(X(1,:)), plot(X(2,:)) 即可。

\textbf{(7) 去眼部伪差：}

题目指出：每个混合信号 x_i(n) (i\in\{1,2\}) 都包含了 “间歇性尖峰 + 眼部活动 + 肌肉活动” 的混合。

我们现在的目标，是用自适应 LMS 算法去除眼部活动。参考噪声 b_2(n)（即题目中的眼部活动参考信号）已在 Ref_Oc1.mat 文件中给出。它与目标中要消除的眼部活动 b_1(n) 有较强相关性。这样，我们就可以把：

• y_1(n) = x_1(n) = \text{(尖波 + 眼部 + 肌肉混合)}
• y_2(n) = b_2(n) = \text{(参考眼部活动)}

在 MATLAB 中，可做以下操作：

1. 载入参考眼部信号 b_2(n)：

   • \text{load('Ref\_Oc1.mat');}

     这里的向量（假设名字叫 refOc）即为 b_2(n)。

2. 执行与第一部分类似的 LMS 步骤：

   • 初始化自适应滤波器的阶数M及步长\mu。
   • 令 \hat{x}_1(n) = x_1(n) - h^T b_2(n) 。
   • 更新公式：
     h(n+1) = h(n) + \mu \,\hat{x}_1(n)\, b_2(n).

3. 最终输出的 \hat{x}_1(n) 应当是去除了大部分眼部伪差后的信号。你可以对 x_2(n) 同理操作。

与之前的题目步骤第3和第4的要求一致，你需要：

• 绘制在一定时间范围内滤波后输出以及误差（若有真值可参考，或仅观察滤波输出），
• 对比不同步长\mu的收敛速度和稳态性能，
• 最后计算稳态时误差（若有真实参考即可做 s(n)-\hat{s}(n)，否则也可仅比较滤波器系数的收敛稳定性）。

\textbf{小结：}

1. 自适应滤波原理

   • 关键在于利用一个可训练的滤波器 h(n) 去估计噪声分量，从而将其从带噪信号中剔除，得到对目标信号的估计。
   • LMS 算法的收敛和稳态性能受步长\mu影响显著。

2. 不同步长\mu对稳态和收敛速度的影响

   • \mu越大，收敛越快，但稳态误差相对较大，若过大则可能发散；
   • \mu越小，稳态误差越小，但收敛慢，需更多迭代才能达到稳定。

3. 稳态时的误差度量

   • 常用做法是在足够长的采样时间后，取最后一段数据（如最后 500 点）的均值，估计稳定下来的均方误差(MSE)。

4. 在真实场景(EEG/癫痫发作间期尖波)中的应用

   • 首先需从混合模型中分离出眼部活动（伪差）；再考虑利用其他方法（例如 ICA）去除剩余的肌肉活动；最后即可得到相对纯净的癫痫尖波信号。
   • 此外，还可以把去噪后的信号用于后续的神经连接性分析。

通过以上步骤，你就可以完成题目中要求的所有子问题，包括在模拟与真实 EEG/癫痫尖波场景下对自适应滤波方法的完整测试、对滤波结果的绘图分析、对不同步长收敛性和稳态误差的讨论，以及最终去除眼部伪差的效果评估。

下面给出第三步：脑连接性部分的详细答案，包括相干函数性质的证明思路，以及在给定向量自回归模型或生理模型下如何估计两个神经群之间的功能连接性（即相干度），并附上在 MATLAB 中可实现的示例代码与结果讨论。

一、相干函数基本性质的证明

题目首先指出，对于两个中心化的随机过程 u(n) 和 v(n)，如果它们在二阶上是平稳的，那么有

其中

\gamma_{uu}(f), \gamma_{vv}(f)分别为 u(n) 和 v(n) 的功率谱密度 (power spectral density, PSD)，而 \gamma_{uv}(f) 为它们的互功率谱密度 (cross power spectral density, CPSD)。

为证明 |\rho_{uv}(f)|^2 \le 1，令

其中 \lambda 是一个复数，并在题目中取

1. 频域表示

令 U(f), V(f), Z(f) 分别为 u(n), v(n), z(n) 的傅里叶变换（或离散傅里叶变换），则有

2. z(n) 的功率谱密度非负

z(n) 的功率谱密度可写作

由于 \gamma_{uv}(f) = \overline{\gamma_{vu}(f)}，且取

可得

1. |\lambda| = 1（因为 \lambda 的模是 1），
2. \lambda\,\gamma_{uv}(f) = |\gamma_{uv}(f)|,
3. \overline{\lambda}\,\gamma_{vu}(f) = |\gamma_{uv}(f)|（同理）。

所以

这里因为 |\lambda|^2 = 1。显然 \gamma_{uu}(f) \ge 0, \gamma_{vv}(f) \ge 0，且 |\gamma_{uv}(f)| 是非负数，因此

3. 推得 |\rho_{uv}(f)|^2 \le 1

回到相干函数定义：

若 |\rho_{uv}(f)|^2 > 1，则构造的 \gamma_{zz}(f) 会变成负值（与非负性矛盾），因此必然 |\rho_{uv}(f)|^2 \le 1。又因为 PSD 与 CPSD 的定义可确保 |\rho_{uv}(f)|^2 \ge 0，故

二、基于相干函数的脑连接性分析：实验部分

题目给出了一个向量自回归 (VAR) 模型以及生理模型生成的 iEEG 信号，要求我们利用 Welch 方法 (pwelch) 和相干函数 (mscohere) 来估计这两个信号之间的功能连接性。

1. 向量自回归模型 (VAR) 及实验设置

题目中的 VAR 模型为

其中 w_1(n), w_2(n) 分别为零均值、方差为 1 的独立白噪声（可在数值仿真中用 randn 生成）。

该模型可理解为：y_1(n) 主要自回归于自身过去两时刻，而 y_2(n) 则部分地依赖于 y_1(n-1)，自身过去值 y_2(n-1) 以及本身的白噪声驱动。

2. 实验 a）：将 y_2(n) 视为 y_1(n) 的延迟版本

题目 a）先简化：只保留

(即只用方程 (1) 的第一行)，而 y_2(n) 不再来自 (1) 的第二行，而是简单地令

1. 生成并绘制各自的频谱
   • 可以采用 pwelch(y, win, overlap, nfft, fs) 命令得到每个信号的功率谱估计；
   • 题目要求用 4 个长度为 1024 点的块（segments），并且 50% 重叠。
2. 计算并绘制这两个信号的幅度平方相干性 (MSC)
   • 在 MATLAB 中用 [Cxy,F] = mscohere(y1, y2, win, overlap, nfft, fs)。
3. 结论：
   • 因为这里 y_2(n) 与 y_1(n) 只有一个简单的线性延迟关系，因此在某些频率上应当出现较高的相干性。
   • 你可以观察到：越是与 y_1(n) 的主能量集中频段吻合的地方（例如它的谐振峰或 AR 共振峰所在频带），两者的相干度越高。

下面给出一个示例 MATLAB 伪代码：

% 设定基本参数
N = 8192;           % 信号长度(>= 4096保证分段)
fs = 1;             % 采样率(若无明确说明,可取1)
nfft = 1024;
window = hanning(nfft);
noverlap = nfft/2;  % 50% 重叠

% 生成 y1(n)
y1 = zeros(1, N);
% 前两点随意初始化
y1(1) = randn; 
y1(2) = randn;
for i = 3:N
    y1(i) = 0.95*sqrt(2)*y1(i-1) - 0.9025*y1(i-2) + randn; 
end

% 生成 y2(n) = 0.5*y1(n-1)
y2 = zeros(1, N);
for i = 2:N
    y2(i) = 0.5*y1(i-1);
end

% 1) 计算并绘制功率谱
[Pyy1, f] = pwelch(y1, window, noverlap, nfft, fs);
[Pyy2, f] = pwelch(y2, window, noverlap, nfft, fs);
figure; 
subplot(2,1,1); plot(f, 10*log10(Pyy1)); title('Power spectrum of y1');
subplot(2,1,2); plot(f, 10*log10(Pyy2)); title('Power spectrum of y2');

% 2) 计算并绘制相干度
[Cxy, f] = mscohere(y1, y2, window, noverlap, nfft, fs);
figure;
plot(f, Cxy, 'LineWidth',1.2);
title('MSC between y1 and y2');
xlabel('Frequency'); ylabel('Coherence');
grid on;

% 根据结果进行讨论: 
% - 哪些频段 MSC 较高？
% - 与 y1(n) 的 AR 特征是否一致？

3. 实验 b）：采用完整的 VAR 方程 (1)

在第二步中，y_1(n) 和 y_2(n) 同时由方程 (1) 的两行共同生成：

1. 生成信号 y_1(n), y_2(n)：
   • 需要两重循环，每一步更新 y_1(n) 和 y_2(n)；
   • w_1(n), w_2(n) 可用 randn 生成。
2. 用 100 个长度为 1024 点的块、50% 重叠来估计频谱和 MSC：
   • 注意此时我们信号长度要足够大（如 N \ge 1024\times(1+\text{重叠系数})\times\text{块数}）。
   • 通过循环或直接使用 pwelch, mscohere 的参数控制。
3. 绘制平均 MSC 并分析结果：
   • 不同块间统计出平均值，然后可视化。
   • 可以观察到，相比于简化版 y_2(n)=0.5\,y_1(n-1)，现在 y_2(n) 同时依赖于 y_1(n-1)、自身过去值 y_2(n-1) 和白噪声 w_2(n)，因此其与 y_1(n) 的相干度分布将更复杂。

示例简要 MATLAB 伪代码（仅展示核心部分）：

N = 200000;  % 信号足够长
y1 = zeros(1,N); y2 = zeros(1,N);
w1 = randn(1,N); w2 = randn(1,N);

% 初始值
y1(1) = randn; 
y1(2) = randn;
y2(1) = randn; 

for i = 2:N
    y1(i) = 0.95*sqrt(2)*y1(i-1) - 0.9025*y1(i-2) + w1(i);
    y2(i) = -0.5*y1(i-1) + 0.25*sqrt(2)*y2(i-1) + w2(i);
end

% 计算 pwelch, mscohere, 并对 100 个块取平均(题目可自行实现)
[Px1, f] = pwelch(y1, 1024, 512, 1024, 1);
[Px2, f] = pwelch(y2, 1024, 512, 1024, 1);
[Cxy, f] = mscohere(y1, y2, 1024, 512, 1024, 1);

% 将 100 个分块计算结果取平均(也可用 pwelch/mscohere 自带的多段功能)
...
% 画图并讨论
figure;
plot(f, Cxy); 
title('Average MSC between y1 and y2 (VAR model)');

分析：

• 在该模型中，y_2(n) 与 y_1(n) 并非简单延迟关系，因而相干峰会更分散，也会出现与 y_2(n) 自身 AR 项关联的频带特征。
• 通过查看 MSC，可了解 y_1(n)、y_2(n) 在不同频率上的耦合强弱，反映了二者的 “功能连接” 程度。

4. 实验 c）：生理模型生成的两个信号

题目最后给出了一个由生理模型生成的两通道信号（Matlab 文件 “ScenarioK1500second1000a.mat” 中的 LFP_s 变量）。令这两个信号分别为 y_1(n) 和 y_2(n)。要求：

1. 估计两信号随时间变化的 MSC：
   • 可以在一定长度的滑动窗上（如 1024 点），并设置一定的步长（如每次移动 256 点），计算出相干度，然后随时间画出。
2. 绘制信号的频谱以及 MSC：
   • 同上，使用 pwelch / mscohere 即可。

示例思路（仅文字说明）：

load('ScenarioK1500second1000a.mat');  % 这里得到 LFP_s, 尺寸 (2 x N)
y1 = LFP_s(1,:);
y2 = LFP_s(2,:);

% 功率谱估计
[P1, f] = pwelch(y1, 1024, 512, 1024, fs);
[P2, f] = pwelch(y2, 1024, 512, 1024, fs);

% 平均 MSC
[Cxy, f] = mscohere(y1, y2, 1024, 512, 1024, fs);

% 若需要“随时间变化”曲线：用段式或滑窗的方式
% 可自行写一个循环, 每段 1024 点, overlap 512点, 对每段独立计算 mscohere
% 再将结果存储到一个矩阵中, 最后用 imagesc 或 surf 之类做时间-频率图

% 例如：
step = 512;  % 窗口步长
winlen = 1024;
idx_start = 1;
time_idx = 1;

while idx_start + winlen - 1 <= length(y1)
    seg1 = y1(idx_start : idx_start+winlen-1);
    seg2 = y2(idx_start : idx_start+winlen-1);

    [Cseg, f] = mscohere(seg1, seg2, [], [], 1024, fs);

    % 存到一个矩阵, 例如 Cmat(:, time_idx) = Cseg;
    Cmat(:, time_idx) = Cseg;

    idx_start = idx_start + step;
    time_idx = time_idx + 1;
end

% 然后可以用 imagesc 或 pcolor 画出随时间-频率的 MSC 分布：
figure; 
imagesc([1:time_idx-1], f, Cmat);
set(gca,'YDir','normal'); colorbar;
xlabel('Time index (sliding windows)'); ylabel('Frequency');
title('Time-varying MSC between y1 and y2');

结果分析：

• 查看在不同时段、不同频率上，两个通道的相干度如何变化；
• 如果信号是生理模型的模拟产物，通常会在特定频带上体现显著连接或耦合。

三、总结

1. 相干函数的定义和性质

   • 相干函数 |\rho_{uv}(f)|^2 衡量了两个平稳过程在频域的线性关系强度；
   • 它必然介于 0 和 1 之间。

2. VAR 模型下的连接性

   • y_1(n), y_2(n) 的生成方式直接决定了它们在频域上的耦合关系；
   • 简单的延迟关系 (如 y_2(n) = a\,y_1(n-1)) 会在某些频率上表现出高相干；
   • 真实场景或更复杂的生理/神经模型会使相干度分布更加丰富。

3. Welch 方法与 mscohere

   • 在实际应用中，常用段式平均 (segmented averaging) 的方式来估计功率谱和互功率谱，得到更平滑的结果；
   • 通过 mscohere 函数，可直接得到幅度平方相干性 (MSC)。
   • 若想分析随时间的连接性变化，可以采用滑动窗 (moving window) 的方法进行时频分析。

通过以上步骤，就可以完成第三步：脑连接性的所有实验及理论推导要求，包括相干函数的基本性质验证，以及在给定模型或真实/生理模型信号下对两通道 iEEG 的频谱和相干度进行估计与讨论。