I. 任务描述
目标
我们要解决医学信号处理中的两大问题：
- 从复杂的脑电信号中提取癫痫相关信号
癫痫患者的大脑活动记录通常会被许多其他信号干扰，比如眼睛的运动和肌肉活动。我们要从这些混合信号中分离出对研究癫痫发作最重要的信号。 - 研究癫痫相关脑信号之间的连接性
通过去噪后的数据，研究大脑中不同区域的 交互联系，以更好地了解癫痫发作时的脑网络行为。
数学模型
假设在观察时间 T 内，在一位癫痫患者头皮上使用的 N_c 个传感器（本实验只用两个传感器 N_c=2）来记录信号，最终观测到三种电生理信号的线性瞬时混合。所获取的数据在空间-时间矩阵中表示，记为 X(N_c \times T) ，它满足以下线性瞬时混合模型：
-
其中，A(N_c \times P) 是混合矩阵，表示源信号到传感器信号的传播关系；S(P \times T) 是源信号矩阵，其第 p 行表示第 p 个源信号的时间变化。
-
注意，在本实验中假设采用无仪器噪声的观测模型。
在第II和第III部分中，会涉及三种生理信号：
- s_1(n) ：癫痫发作间期尖波，是大脑中神经群活动的目标信号。
- s_2(n) ：眼部运动，作为伪差之一。
- s_3(n) ：肌肉活动，作为伪差之二。
其中，s_2(n) 和 s_3(n) 被视为伪差信号，它们干扰了对 s_1(n) 的提取，因此需要在后续处理中被分离并去除。
去伪差与分析步骤描述
- 眼部运动去除：
首先通过自适应滤波技术从时空混合信号 X 中分离并去除眼部运动伪差 s_2(n)。自适应滤波器能动态调整自身参数，以适应输入信号的变化，从而有效减小伪差的影响。 - 肌肉活动去除：
在去除眼部运动后，剩余信号中的肌肉活动伪差 s_3(n) 将通过独立成分分析（ICA）技术去除。ICA 是一种能够分解混合信号为独立源信号的统计方法，特别适合分离不同来源的生理信号。 - 大脑连接性研究：
在去除所有伪差后，提取的目标信号 s_1(n) 被用于研究模拟癫痫发作时的大脑连接性问题。这一部分探索了神经群之间的相互作用，为癫痫的病理研究提供了数据支持。
II. 使用带参考通道的自适应滤波算法去除观测信号中的眼部伪差
这一部分的目标是从两个传感器输出中去除眼部活动（伪差 s_2(n) ）。为此，需要考虑两套伪差去除系统(每个通道对应一套，即每个传感器输出对应一个伪差去除系统)，下面为了更清晰的描述，我们只详细介绍单个传感器的伪差去除系统，因此考虑另一个传感器的时候，采用相同的步骤即可。
降噪系统中可以观测到两个信号：
-
观测信号 y_1(n)，包含目标信号 s(n) 和噪声 b_1(n) （要从观测信号 y_1(n) 中去除），满足：
y_1(n) = s(n) + b_1(n) -
参考信号 y_2(n)，直接是噪声信号 b_2(n) 的观测：
y_2(n) = b_2(n)
这里的假设是：
- 参考噪声 b_2(n) 和目标噪声/伪差 b_1(n) 相关
- b_2(n) 是 b_1(n) 的某种滤波版本。
目标是设计一个 自适应滤波器 h，使得通过最小化噪声的估计误差 E\bigl[(b_1 - h(b_2))^2\bigr] ，即均方误差 E\bigl[\bigl(b_1(n) - h^T(n)\,b_2(n)\bigr)^2\bigr] ，最终得到去噪的目标信号估计 \hat{s}(n)。
自适应滤波器算法 数学描述及流程
在实际应用中，滤波器 h 是一个具有 M 个系数(本实验中 M=20 )的横向滤波器（FIR），其目的是利用已知的参考信号 b_2 来估计噪声/伪差 b_1 ，从而得到对目标信号 s 的最佳估计 \hat{s} 。下面给出了描述LMS算法的主要公式：
滤波后的输出信号：
自适应滤波器 h 的输出为：
- h^T(n)：滤波器当前的权值向量，表示为 [h_1(n), h_2(n), \dots, h_M(n)]
- b_2(n)：参考信号的延迟窗口，表示为 [b_2(n), b_2(n-1), \dots, b_2(n-M+1)]
- 利用参考信号 b_2(n) 的延迟窗口（即过去 M 个时刻的值）来估计目标噪声 b_1(n)
估计目标信号：
从观测信号中去除滤波器输出后，得到目标信号的估计：
- 其中，y_1(n) = s(n) + b_1(n)。
权值更新公式：
滤波器通过最小化估计误差 \hat{s}(n) 的均方误差（MSE），自适应更新权值：
- \mu：步长因子，控制更新速度和稳定性。
- \hat{s}(n) 和 b_2(n) 的乘积用于调整权值。
该降噪系统的示意图如下：
其中:
-
Reference noise 代表眼部活动这一参考噪声信号
-
Estimated signal 表示去除了眼部活动后的估计信号
问题
1. 在测试真实癫痫信号之前，我们先在完全模拟的数据上测试此自适应算法。给定的嘈杂模拟数据如下：
- 信号 s(n) 是振幅 A=2 、频率为100 Hz、采样率为1 kHz的正弦信号。
- b_2(n) 是均值为0、方差为1的白高斯噪声。
- b_1(n) 则通过一个传递函数为 1/(1-\alpha z^{-1}) 的滤波器对 b_2(n) 进行滤波得到。
请编写对应的程序，并在 \alpha=0.9 、 M=20 以及不同 \mu 取值（例如 \mu=0.01 、 \mu=0.0005 ）时进行测试。
我们先生成目标信号 s(n)
n = 0:N-1; % 生成样本索引序列，范围为 [0, N-1]
A = 2; % 振幅设为 2
s = A * sin(2*pi*f0*(n/fe)); % 目标信号为正弦信号
- f_0=100 Hz 为信号频率，f_e=1000 Hz 为采样率。
在生成噪声信号 b_2(n)，是一个标准白高斯噪声。
b2 = randn(1, N); % 生成长度为 N 的白高斯噪声，均值为 0，方差为 1
通过滤波生成 b_1(n)
b1 = filter([1], [1, -alpha], b2);
使用一个 IIR 滤波器对 b_2(n) 进行处理，其传递函数为 H(z) = 1 / (1 - \alpha z^{-1})。滤波器的系数为alpha=0.9 ，它决定了滤波器的“记忆”特性，也即 b_1(n) 的平滑程度。
现在我们要实现 LMS 自适应滤波算法
初始化参数
% 初始化自适应滤波器变量
h = zeros(M, 1); % 滤波器系数
s_hat = zeros(1, N); % 目标信号估计
e = zeros(1, N); % 误差
phi = zeros(1, N); % 均方误差（误差平方）
M=20是滤波器的阶数，决定了参考向量的长度- 误差
e的定义为 e(n) = s(n) - s_{hat}(n)
核心循环
% LMS 循环
for n = M:N
% 提取参考信号向量 b2_ref
b2_ref = flipud(b2(n-M+1:n).'); % 参考信号向量，长度为 M
% 噪声估计
b1_hat = h' * b2_ref;
% 目标信号估计
s_hat(n) = x(n) - b1_hat;
% 误差计算
e(n) = s(n) - s_hat(n);
% 更新滤波器系数
h = h + mu * s_hat(n) * b2_ref;
% 记录误差平方
phi(n) = e(n)^2;
end
首先我们从第 M 个样本开始，因为要生成参考向量 b_2^{ref}(n) 需要过去 M 个样本。
然后我们提取 从 b_2(n-M+1) 到 b_2(n) 的参考向量b2_ref，并将其翻转为列向量。
b2_ref = flipud(b2(n-M+1:n).'); % 参考信号向量，长度为 M
然后我们利用滤波器权重 h 和参考向量 b_2^{ref}(n) 的点积得到估计噪声 \hat{b}_1(n)
b1_hat = h' * b2_ref;
从 x(n) 中减去估计的 \hat{b}_1(n) ，得到估计的目标信号 \hat{s}(n)。
s_hat(n) = x(n) - b1_hat;
随后计算误差 e(n) = s(n) - \hat{s}(n)
e(n) = s(n) - s_hat(n);
使用 LMS 算法更新权重：
h = h + mu * s_hat(n) * b2_ref;
完整代码与结果
% 初始化参数
clear; close all; clc;
% 基本设置
fe = 1000; % 采样率
f0 = 100; % 信号频率
N = 20000; % 样本数
M = 20; % 自适应滤波器阶数
A = 2; % 振幅
alpha = 0.9; % 滤波器系数
mu = 0.001; % 步长
% 数据生成
s = A * sin(2 * pi * f0 * (0:N-1) / fe); % 目标信号 s(n)
b2 = randn(1, N); % 高斯白噪声 b2(n)
b1 = filter([1], [1, -alpha], b2); % 滤波后的噪声 b1(n)
x = s + b1; % 带噪信号 x(n)
% 初始化自适应滤波器变量
h = zeros(M, 1); % 滤波器系数
s_hat = zeros(1, N); % 目标信号估计
e = zeros(1, N); % 误差
phi = zeros(1, N); % 均方误差（误差平方）
% LMS 循环
for n = M:N
% 提取参考信号向量 b2_ref
b2_ref = flipud(b2(n-M+1:n).'); % 参考信号向量，长度为 M
% 噪声估计
b1_hat = h' * b2_ref;
% 目标信号估计
s_hat(n) = x(n) - b1_hat;
% 误差计算
e(n) = s(n) - s_hat(n);
% 更新滤波器系数
h = h + mu * s_hat(n) * b2_ref;
% 记录误差平方
phi(n) = e(n)^2;
end
% 绘图
figure('Name', 'LMS Filtering Results');
% (1) 原始目标信号 s(n)
subplot(6, 1, 1);
plot(s, 'LineWidth', 1.2);
title('s(n): 原始目标信号');
ylabel('Amplitude');
grid on;
% (2) 滤波后的噪声 b1(n)
subplot(6, 1, 2);
plot(b1, 'LineWidth', 1.2);
title('b_1(n): 滤波后的噪声');
ylabel('Amplitude');
grid on;
% (3) 带噪信号 x(n)
subplot(6, 1, 3);
plot(x, 'LineWidth', 1.2);
title('x(n): 带噪信号');
ylabel('Amplitude');
grid on;
% (4) 前 200 样本的估计信号 s_hat
subplot(6, 1, 4);
plot(1:200, s_hat(1:200), 'LineWidth', 1.2);
title('\hat{s}(n): 前 200 样本的估计信号');
ylabel('Amplitude');
grid on;
% (5) 最后 200 样本的估计信号 s_hat
subplot(6, 1, 5);
plot(1:200, s_hat(N-199:N), 'LineWidth', 1.2);
title('\hat{s}(n): 最后 200 样本的估计信号');
ylabel('Amplitude');
grid on;
% (6) 误差平方 \phi(n)
subplot(6, 1, 6);
plot(phi, 'LineWidth', 1.2);
title('\phi(n): 误差平方');
xlabel('Sample Index');
ylabel('Error^2');
grid on;
最后结果显示:
- 第一行：绘制原始目标信号 s(n)，用于观察正弦信号的特性。
- 第二行：绘制滤波后的噪声 b_1(n)，模拟眼动伪差。
- 第三行：绘制带噪信号 x(n)，是 s(n) 和 b_1(n) 的叠加。
- 第四行：绘制估计信号 \hat{s}(n) 的前 200 个样本，用于观察滤波器初期性能。
- 第五行：绘制估计信号 \hat{s}(n) 的最后 200 个样本，用于观察稳态性能。
- 第六行：绘制误差平方 \phi(n)，展示误差随时间的变化和收敛趋势。
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 的情况下，针对所有时间点，绘制 \bigl(s-\hat{s}\bigr)^2 这一误差项。
%% 带参考通道的自适应滤波 (LMS) 去噪测试
% 问题 2 & 问题 3 综合示例
clear; clc; close all;
%% 1. 参数设置
fe = 1000; % 采样率
f0 = 100; % 信号频率
N = 20000; % 样本数
M = 20; % 自适应滤波器阶数
A = 2; % 目标信号 s(n) 振幅
alpha = 0.9; % IIR 滤波器系数, 用于生成 b1(n)
% 要测试的不同步长 mu
muSet = [0.01, 0.005, 0.001, 0.0005];
%% 2. 数据生成
n = 0:N-1;
% 目标信号 s(n): 正弦波
s = A * sin(2 * pi * f0 * (n / fe));
% 白高斯噪声 b2(n)
b2 = randn(1, N);
% IIR 滤波器 H(z) = 1 / (1 - alpha z^-1), 生成 b1(n)
b1 = filter([1], [1, -alpha], b2);
% 带噪信号 x(n) = s(n) + b1(n)
x = s + b1;
%% 3. 针对每个 mu 进行 LMS 滤波，并绘图
for i = 1:length(muSet)
mu = muSet(i);
% (a) 初始化自适应滤波器变量
h = zeros(M, 1); % 滤波器系数，长度 M
s_hat = zeros(1, N); % 目标信号估计
e = zeros(1, N); % 误差
% (b) LMS 迭代
for k = M:N
% 取参考信号向量 (翻转成列向量)
b2_ref = flipud(b2(k - M + 1 : k).');
% 估计噪声 b1_hat
b1_hat = h' * b2_ref;
% 估计目标信号 s_hat(k) = x(k) - b1_hat
s_hat(k) = x(k) - b1_hat;
% 计算与真实目标信号 s(k) 的误差
e(k) = s(k) - s_hat(k);
% 更新滤波器系数
h = h + mu * s_hat(k) * b2_ref;
end
%% 问题 2: 绘制输出信号在 1~100 和 5000~5120 区间的波形
figure('Name', ['LMS Output - \mu = ', num2str(mu)]);
% 1) 估计信号前 100 个样本
subplot(2,1,1);
plot(1:100, s_hat(1:100), 'LineWidth', 1.2);
title(['\hat{s}(n) (前 100 个样本), \mu = ', num2str(mu)]);
xlabel('样本点'); ylabel('幅值'); grid on;
% 2) 估计信号在 5000~5120 区间
subplot(2,1,2);
plot(5000:5120, s_hat(5000:5120), 'LineWidth', 1.2);
title(['\hat{s}(n) (5000~5120 区间), \mu = ', num2str(mu)]);
xlabel('样本点'); ylabel('幅值'); grid on;
%% 问题 3: 针对所有时间点，绘制 (s - s_hat)^2
figure('Name', ['Error squared - \mu = ', num2str(mu)]);
plot((s - s_hat).^2, 'LineWidth', 1.2);
title(['(s - \hat{s})^2, 全部时间点, \mu = ', num2str(mu)]);
xlabel('样本点'); ylabel('(s - \hat{s})^2'); grid on;
drawnow; % 更新图像
end
4. 讨论步长 \mu 在稳态时和在收敛速度方面所产生的影响。
在自适应滤波中，步长 \mu 决定了滤波器系数更新幅度。
- 当 \mu 过大时，滤波器系数更新幅度较大，初始阶段的收敛速度很快，能迅速减小误差，但是由于更新过于剧烈，滤波器的输出有可能在误差附近振荡、不稳定，即稳态时往往残留的误差较大。
- 当 \mu 较小时，收敛速度变慢，需要更长时间（更多迭代次数）才能达到稳态，但一般会带来更好的稳态精度，获取更小的均方误差 (MSE)
因此步长 \mu 的选择需要平衡 收敛速度 和 稳态误差
5. 通过对最后500个采样点的平均，计算稳态时的误差。
为了观察自适应滤波在接近稳态时的表现，我们取最后500个采样点的误差取平均，从而估计稳态误差大小。真实信号是 s ，估计信号是 \hat{s} ，那么可用以下方式计算稳态误差（以均方误差 MSE 为例）：
在 MATLAB 中可直接写：
steady_error = mean((s(N-499:N) - s_hat(N-499:N)).^2);
根据收敛后期的平均误差情况，可以对不同步长 \mu 的的稳态性能进行比较。
6. 在本问题中，所关注的信号 s(n) 可以是一串癫痫发作间期尖波。
a) 使用 signals.mat 文件，通过 MATLAB 的 load 命令加载三个信号：
- 信号
s1表示 间歇性尖峰 (interictal spike)， - 信号
s2表示 眼部运动 (ocular movement)， - 信号
s3表示 肌肉活动 (muscle activity)。
dataStruct = load('signals.mat');
s1 = dataStruct.s1;
s2 = dataStruct.s2;
s3 = dataStruct.s3;
S = [s1(:)'; s2(:)'; s3(:)'];
信号矩阵 S，大小为 3 \times T，每一行对应一个信号。
b) 构建时空混合矩阵 X(2 \times T) ，其公式如下：
其中:
-
s_i(:,:) 是矩阵 S = [s_1^T, s_2^T, s_3^T]^T 的第 i 行
-
向量 A = \left[a_{1}, a_{2}, a_{3}\right] = \begin{bmatrix}1 & 0.3 & 0.5 \\ -0.7 & 0.6 & 0.5\end{bmatrix}
% 混合系数矩阵 A A = [1 0.3 0.5; -0.7 0.6 0.5]; % 提取矩阵 A 的列向量 a1 = A(:,1); % 第一列向量 a1 a2 = A(:,2); % 第二列向量 a2 a3 = A(:,3); % 第三列向量 a3 -
\sigma_1, \sigma_2 是用来调整信号噪声比 (SNR) 的标量
这里简化假设 \sigma_1 = \sigma_2 = \sigma ，其中：
- SNR = 20 \log_{10} (\frac{\sigma_s}{\sigma})
- \sigma_s 是信号的标准差
- 我们设 \sigma_s = 1 且 SNR = 5 dB
sigma_s = 1;
SNR_dB = 5;
这里缩放系数 \sigma 的计算公式为:
sigma = sigma_s / (10^(SNR_dB / 20));
构造混合信号矩阵 X
回顾上述混合信号矩阵 X 的构建公式，（别忘了对信号 S(1,:)、S(2,:) 和 S(3,:) 进行归一化）
tmp1 = a1 * S(1, :) / norm(a1 * S(1, :));
tmp2 = a2 * S(2, :) / norm(a2 * S(2, :));
tmp3 = a3 * S(3, :) / norm(a3 * S(3, :));
X = tmp1 + sigma * tmp2 + sigma * tmp3;
c) 绘制信号 s_1(n) 、 s_2(n) 、 s_3(n) ，以及混合信号 x_1(n) 和 x_2(n) （即 X(1,:) 和 X(2,:) ）
% 绘制源信号、带噪信号和去噪信号
figure('Name','5 行结果展示');
% 第一行：原始信号 s1 (b1)
subplot(5, 1, 1);
plot(s1, 'LineWidth', 1.2);
title('Signal s_1 (b1)');
ylabel('Amplitude');
grid on;
% 第二行：原始信号 s2 (spk)
subplot(5, 1, 2);
plot(s2, 'LineWidth', 1.2);
title('Signal s_2 (spk)');
ylabel('Amplitude');
grid on;
% 第三行：原始信号 s3 (ocl)
subplot(5, 1, 3);
plot(s3, 'LineWidth', 1.2);
title('Signal s_3 (ocl)');
ylabel('Amplitude');
grid on;
% 第四行：带噪声混合信号 x1
subplot(5, 1, 4);
plot(X(1,:), 'LineWidth', 1.2);
title('Mixed Signal x_1 (With Noise)');
ylabel('Amplitude');
grid on;
% 第五行：带噪声混合信号 x2
subplot(5, 1, 5);
plot(X(2,:), 'LineWidth', 1.2);
title('Mixed Signal x_2 (With Noise)');
ylabel('Amplitude');
grid on;
7. 根据上述混合模型，每个混合信号 x_i(n) ( i \in \{1, 2\} ) 表示一个 混合信号，它由 间歇性尖峰 和 记录于头皮表面上的眼肌活动 组成。
在本部分中，参考噪声信号（眼部活动） b_2(n) 是与目标眼部活动 b_1(n) （例如记录的实际眼部活动）相关的信号，并已提供于 MATLAB 文件 Ref_Oc1.mat 中
编写程序，通过自适应 LMS 算法来去除混合信号中眼部活动对间歇性尖峰和肌肉伪影的贡献。
从矩阵 X 中提取两个混合信号 x_1 和 x_2，分别对应第一通道和第二通道。
x1 = X(1,:);
x2 = X(2,:);
加载参考信号 ref_ocl.mat，变量名为 ocl_n，作为参考眼动信号 b_2(n)。
refData = load('ref_ocl.mat');
refOc = refData.ocl_n;
确保参考信号与混合信号长度一致，若长度不足则补齐或插值，若信号过长，仅截取前 T 个样本。非必要。
if length(refOc) < T
error('参考眼动信号长度不足，请检查并补齐或使用插值。');
elseif length(refOc) > T
refOc = refOc(1:T);
end
将参考信号归一化到 [-1, 1] 区间，避免过大的数值导致算法不稳定。
refOc = refOc / max(abs(refOc));
LMS 算法部分
**通道 1 的去噪过程 **
设置滤波器阶数 M = 20 和步长因子 \mu = 0.05。步长因子影响收敛速度和稳定性。
M = 20;
mu = 0.05;
初始化滤波器系数 h_1、去噪输出 x_{1D} 和误差信号 e_1。
h1 = zeros(M, 1);
x1_hat = zeros(1, T);
e1 = zeros(1, T);
下面进入循环部分
for n = M:T
取参考信号中最近 M 个样本，构成长度为 M 的窗口，按时间顺序反转为列向量。
refWindow = flipud(refOc(n-M+1:n).');
通过当前滤波器系数 h_1 与参考信号窗口相乘，得到估计的眼动噪声 noiseEst
noiseEst = h1' * refWindow;
从混合信号 x_1 中减去估计噪声，得到去噪信号 x1_hat(n)。
x1_hat(n) = x1(n) - noiseEst;
计算当前误差 e1(n)，定义为混合信号与去噪信号的差值。
e1(n) = x1(n) - x1_hat(n);
根据 LMS 算法公式，调整滤波器系数 h_1，以逐步减少误差。
h1 = h1 + mu * x1_hat(n) * refWindow;
通道 2 的去噪过程
初始化第二通道的滤波器参数、去噪输出和误差信号。
h2 = zeros(M, 1);
x2_hat = zeros(1, T);
e2 = zeros(1, T);
通道 2 的去噪过程与通道 1 类似，但针对混合信号 x_2 和滤波器 h_2 分别进行。
for n = M:T
refWindow = flipud(refOc(n-M+1:n).');
noiseEst = h2' * refWindow;
x2_hat(n) = x2(n) - noiseEst;
e2(n) = x2(n) - x2_hat(n);
h2 = h2 + mu * x2_hat(n) * refWindow;
end
完整代码及结果展示
x1 = X(1,:);
x2 = X(2,:);
% 加载参考眼动信号 (例如来自实际采集的眼动传感器)，此处仅做示例
% 假设文件 Ref_Oc1.mat 中包含一个变量 refOc，对应眼动参考信号
refData = load('ref_ocl.mat');
% 假设文件里面的变量名叫 refOc，也可能是别的名字，请根据实际情况替换
refOc = refData.ocl_n;
% 为了对齐，这里假设 refOc 的长度与 x1,x2 一样，若长度不同，可酌情裁剪/插值
if length(refOc) < T
error('参考眼动信号长度不足，请检查并补齐或使用插值。');
elseif length(refOc) > T
refOc = refOc(1:T); % 仅截取与 X 长度一致的部分
end
% 如果 refOc 过大或过小，也可以归一化或其他预处理
% 比如简单归一化:
refOc = refOc / max(abs(refOc));
% 使用自适应LMS算法去除眼动伪差 (逐通道处理)
% 对 x1 进行滤波，再对 x2 进行滤波
% 每个通道的流程与前述“带参考通道的LMS去噪”步骤相同
% ========== 设置自适应滤波器参数 ==========
M = 20; % 滤波器阶数
mu = 0.05; % 步长因子 (可多次测试不同 mu)
% ========== 通道1: x1 去眼动伪差 ==========
% 初始化
h1 = zeros(M, 1); % 滤波器系数
x1_hat = zeros(1, T); % 去噪后输出
e1 = zeros(1, T); % 误差信号
% LMS迭代
for n = M:T
% 取参考信号 refOc 的 M 长度向量，并翻转成列
refWindow = flipud(refOc(n-M+1:n).');
% 估计出的眼动噪声
noiseEst = h1' * refWindow;
% 去噪输出
x1_hat(n) = x1(n) - noiseEst;
% 误差: (真值 - 估值)，此处若只想最小化 x1 中眼动成分，可将误差定义成 x1_hat(n)本身
% 但如需对比真值 s1(n) 时，要看已知不已知 s1。这里示例按常见形式:
e1(n) = x1(n) - x1_hat(n);
% 更新系数
h1 = h1 + mu * x1_hat(n) * refWindow;
end
% ========== 通道2: x2 去眼动伪差 ==========
h2 = zeros(M, 1);
x2_hat = zeros(1, T);
e2 = zeros(1, T);
for n = M:T
refWindow = flipud(refOc(n-M+1:n).');
noiseEst = h2' * refWindow;
x2_hat(n) = x2(n) - noiseEst;
e2(n) = x2(n) - x2_hat(n);
h2 = h2 + mu * x2_hat(n) * refWindow;
end
X1D = x1_hat;
X2D = x2_hat;
figure('Name','7 行结果展示');
% 第一行：原始信号 s1 (b1)
subplot(7, 1, 1);
plot(s1, 'LineWidth', 1.2);
title('信号 s_1 (b1)');
ylabel('幅值');
grid on;
% 第二行：原始信号 s2 (spk)
subplot(7, 1, 2);
plot(s2, 'LineWidth', 1.2);
title('信号 s_2 (spk)');
ylabel('幅值');
grid on;
% 第三行：原始信号 s3 (ocl)
subplot(7, 1, 3);
plot(s3, 'LineWidth', 1.2);
title('信号 s_3 (ocl)');
ylabel('幅值');
grid on;
% 第四行：带噪声混合信号 x1
subplot(7, 1, 4);
plot(X(1,:), 'LineWidth', 1.2);
title('混合信号 x_1 (含噪声)');
ylabel('幅值');
grid on;
% 第五行：带噪声混合信号 x2
subplot(7, 1, 5);
plot(X(2,:), 'LineWidth', 1.2);
title('混合信号 x_2 (含噪声)');
ylabel('幅值');
grid on;
% 第六行：去噪后的 x1 (记为 x_{1D})
subplot(7, 1, 6);
plot(X1D, 'LineWidth', 1.2);
title('去噪后的信号 x_{1D}');
ylabel('幅值');
grid on;
% 第七行：去噪后的 x2 (记为 x_{2D})
subplot(7, 1, 7);
plot(X2D, 'LineWidth', 1.2);
title('去噪后的信号 x_{2D}');
xlabel('采样点');
ylabel('幅值');
grid on;
第二步：使用 ICA 去除肌肉运动
我们现在的目标是从结果混合信号 \hat{X}_{(2 \times N)} = [\hat{x}_1^T, \hat{x}_2^T]^T 中去除剩余的伪影（即肌肉运动）。 \hat{X} 是由 iEEG（脑电活动） 和 肌肉运动 两种不同的生理信号源混合而成，即 。由于这两种信号在统计上是相互独立的，可以使用 独立成分分析（ICA） 技术来分离它们。这种分离过程是解决盲源分离（Blind Source Separation, BSS）问题的核心任务。
具体实现中，我们使用 SOBI 方法（Second-Order Blind Identification，第二阶盲识别方法）。SOBI 利用信号在不同时间延迟下的统计特性，结合协方差矩阵的联合对角化来实现信号分离。
准备工作 -- 预白化
使用 PCA 对可用观测数据 \hat{X}_{(2 \times N)} 进行预白化：
为了让分离器更高效地工作，我们需要先对输入数据信号 \hat{X}_{(2 \times N)} 进行标准化处理，即预白化。预白化让信号变得更规整、简单，可以让数据变得统计上无关，且每个信号的方差归一化为 1，方便后续处理。具体过程如下：
a. 计算零均值观测 \hat{X}_c 的协方差矩阵：
我们从混合信号 \hat{X} 中减去其均值，得到零均值的信号矩阵 \hat{X}_c ，然后计算其协方差矩阵：
- E[\cdot] 表示期望值运算， N 是信号样本的数量，协方差矩阵 R_{\hat{X}} 描述了信号在不同通道之间的相关性。
零均值化
首先，将观测信号矩阵 \hat{X} 的每一列减去其均值向量 \mathbf{m}_{\hat{X}}：
m_X = mean(X, 2); % 每行（通道）的均值
X_c = X - m_X; % 零均值信号
协方差矩阵计算
协方差矩阵 R_{\hat{X}} 定义为：
R_X = (X_c * X_c') / N; % 协方差矩阵
这是一个 2 \times 2 的对称矩阵，表示通道之间的相关性和自相关性：
- \mathrm{Var}(\hat{x}_1) 是第 1 通道信号的方差。
- \mathrm{Cov}(\hat{x}_1, \hat{x}_2) 是两通道信号的协方差。
b. 计算协方差矩阵 \hat{R}_{\hat{X}} 的特征值分解（EVD）：
其中:
-
R_{\tilde{s}} 是新信号源向量 \tilde{s}(n) = [s_1(n), s_2(n)]^T 的协方差矩阵
-
U 是特征向量矩阵， \Lambda 是对角矩阵，其对角线上的元素为对应的特征值
-
同时，我们将矩阵 U 的列重新排序（按特征值从大到小的顺序），得到矩阵 \tilde{U}
在代码中，计算协方差矩阵的特征值和特征向量，用于对信号进行线性变换
[U, Lambda] = eig(R_X); % 计算特征值和特征向量
对特征值按降序排列，并重新排序对应的特征向量
[lambda_sorted, sort_idx] = sort(diag(Lambda), 'descend'); % 特征值降序排列
U_sorted = U(:, sort_idx); % 排序后的特征向量矩阵
Lambda_sorted = diag(lambda_sorted); % 排序后的特征值对角矩阵
显示排序后的特征值和特征向量
disp('排序后的特征值矩阵 Lambda_sorted:');
disp(Lambda_sorted);
disp('排序后的特征向量矩阵 U_sorted:');
disp(U_sorted);
c. 设置预白化算子：
- 其中 \dagger 是伪逆算子（在 Matlab 中可以使用
pinv命令完成此操作）， Q 是一个酉矩阵。
计算特征值矩阵的平方根倒数，用于标准化每个主成分的方差
Lambda_inv_sqrt = diag(1 ./ sqrt(lambda_sorted)); % 特征值的平方根倒数
计算排序后的特征向量矩阵的伪逆，用于构造预白化算子
U_dagger = pinv(U_sorted); % 特征向量矩阵的伪逆
结合特征值矩阵和特征向量矩阵构造预白化算子 W
W = Lambda_inv_sqrt * U_dagger; % 预白化算子
显示预白化算子
disp('预白化算子 W:');
disp(W);
d. 将原始数据矩阵 \hat{X}_c 投影到新子空间：
通过预白化算子 W ，将原始数据投影到新空间，得到预白化后的数据 Z ：
- 这里 Z 表示新的数据矩阵，已去除通道间的相关性，便于后续分离独立信号。
使用预白化算子将零均值化信号 \hat{X}_c 投影到新子空间，得到预白化信号 Z
Z = W * X_c; % 投影到新空间
绘制预白化信号的每个通道，观察其变化
figure('Name', '预白化后的信号');
subplot(2, 1, 1);
plot(Z(1, :));
title('预白化信号通道 1');
xlabel('样本点'); ylabel('幅值'); grid on;
subplot(2, 1, 2);
plot(Z(2, :));
title('预白化信号通道 2');
xlabel('样本点'); ylabel('幅值'); grid on;
验证预白化后的信号是否达到单位协方差矩阵，即统计无关性和方差归一化
R_Z = (Z * Z') / N; % 计算预白化信号的协方差矩阵
disp('预白化信号的协方差矩阵 R_Z (应接近单位矩阵):');
disp(R_Z);
完整代码及结果
%% 预白化全过程：从协方差矩阵到预白化信号投影
clear; clc; close all;
% 输入信号和初始化
% 模拟一个 2 通道的混合信号矩阵 (零均值化前的原始信号)
N = 1000; % 样本数量
X = [randn(1, N) + 2; 0.5 * randn(1, N) - 1]; % 模拟两个通道的混合信号
% 显示原始信号
figure('Name', '原始信号');
subplot(2, 1, 1);
plot(X(1, :));
title('原始信号通道 1');
xlabel('样本点'); ylabel('幅值'); grid on;
subplot(2, 1, 2);
plot(X(2, :));
title('原始信号通道 2');
xlabel('样本点'); ylabel('幅值'); grid on;
% 计算协方差矩阵
% 计算零均值信号
m_X = mean(X, 2); % 每行（通道）的均值
X_c = X - m_X; % 零均值信号
% 计算协方差矩阵
R_X = (X_c * X_c') / N; % 协方差矩阵
% 显示协方差矩阵
disp('协方差矩阵 R_X:');
disp(R_X);
% 特征值分解 (EVD)
% 计算协方差矩阵的特征值和特征向量
[U, Lambda] = eig(R_X);
% 排序特征值和特征向量（按特征值降序）
[lambda_sorted, sort_idx] = sort(diag(Lambda), 'descend'); % 特征值降序排列
U_sorted = U(:, sort_idx); % 排序后的特征向量矩阵
Lambda_sorted = diag(lambda_sorted); % 排序后的特征值对角矩阵
% 显示排序后的特征值和特征向量
disp('排序后的特征值矩阵 Lambda_sorted:');
disp(Lambda_sorted);
disp('排序后的特征向量矩阵 U_sorted:');
disp(U_sorted);
% 设置预白化算子
% 计算特征值矩阵的平方根倒数
Lambda_inv_sqrt = diag(1 ./ sqrt(lambda_sorted));
% 计算特征向量矩阵的伪逆
U_dagger = pinv(U_sorted);
% 构造预白化算子 W
W = Lambda_inv_sqrt * U_dagger;
% 显示预白化算子
disp('预白化算子 W:');
disp(W);
% 投影到新子空间
% 使用预白化算子将信号投影到新空间
Z = W * X_c;
% 显示预白化后的信号
figure('Name', '预白化后的信号');
subplot(2, 1, 1);
plot(Z(1, :));
title('预白化信号通道 1');
xlabel('样本点'); ylabel('幅值'); grid on;
subplot(2, 1, 2);
plot(Z(2, :));
title('预白化信号通道 2');
xlabel('样本点'); ylabel('幅值'); grid on;
% 验证预白化信号的协方差矩阵
R_Z = (Z * Z') / N;
% 显示预白化信号的协方差矩阵
disp('预白化信号的协方差矩阵 R_Z (应接近单位矩阵):');
disp(R_Z);
SOBI 方法去除肌肉活动伪影
在预白化之后，我们使用 SOBI 方法 从信号 Z 中分离出 iEEG 信号和肌肉运动信号。SOBI 的核心思想是利用信号在不同时间延迟下的协方差矩阵，找到一组能够分离信号的变换矩阵 Q 。以下是具体的实现过程：
a. 对时间延迟 \tau_k （ k=1, 2, \dots, 10 ）进行二阶矩估计，计算协方差矩阵：
- R_Z(\tau_k) 是时间延迟为 \tau_k 时的协方差矩阵
- z[n] 是第 n 个采样点的信号
为时间延迟范围 \tau_k = 1, 2, \dots, 10 初始化，准备存储延迟协方差矩阵 R_Z(\tau_k)
tau_max = 10; % 最大时间延迟
N_channels = size(Z, 1); % 信号通道数量
R_tau = zeros(N_channels, N_channels, tau_max); % 存储不同时间延迟下的协方差矩阵
循环计算每个时间延迟 \tau_k 下的协方差矩阵 R_Z(\tau_k)，通过信号时间偏移实现延迟计算
% 计算每个时间延迟的协方差矩阵
for tau = 1:tau_max
R_tau(:, :, tau) = (Z(:, 1:(n - tau)) * Z(:, (tau + 1):n)') / (n - tau);
end
b. 使用提供的联合对角化算法 jad.m 将上述所有协方差矩阵 R_Z(\tau_k) 变换为对角矩阵，从而估计混合矩阵 Q
Q_est = jad(R_tau); % 调用联合对角化算法，返回分离矩阵 Q_est
c. 通过估计得到的混合矩阵 \hat{Q} ，对预白化数据 Z 进行变换，得到包含 iEEG 活动和肌肉运动的源信号向量：
- \hat{Q} 是 Q 的估计。
- \hat{S} 是提取后的信号矩阵，其中每一行对应一个独立的信号源
- 一行是 iEEG 信号，另一行是肌肉运动信号。
% 使用分离矩阵提取源信号
S_hat = Q_est' * Z; % 提取源信号
disp('分离矩阵 Q_est:'); disp(Q_est);
d. 将提取的信号（即 \hat{S} 的各行）（即 iEEG 和肌肉运动）绘制出来，验证分离结果。通过观察信号的特性，可以大致判断 iEEG 和肌肉运动是否被成功分离。
绘制分离出的信号 \hat{S} 的每一行，观察其特性以判断分离效果
% 绘制分离后的信号
figure('Name', '分离后的信号');
subplot(2, 1, 1);
plot(S_hat(1, :)); title('提取信号 1 (iEEG)'); xlabel('样本点'); ylabel('幅值'); grid on;
subplot(2, 1, 2);
plot(S_hat(2, :)); title('提取信号 2 (肌肉运动)'); xlabel('样本点'); ylabel('幅值'); grid on;
e. 使用 Matlab 中的 corr 命令计算真实的 iEEG 信号 s_1(n) 和提取信号 \hat{s}_1(n) 之间的相关系数，并对结果进行解释说明。
使用 MATLAB 的 corr 函数计算真实信号 s_1(n) 和提取信号 \hat{s}_1(n) 的相关系数
% 计算真实信号与分离信号的相关系数
rho = corr(s1', S_hat(1, :)'); % 计算相关系数
disp(['真实信号与提取信号的相关系数 ρ: ', num2str(rho)]);
-
如果 \rho 接近 1，则说明分离效果良好，提取的信号 \hat{s}_1(n) 很好地反映了真实的 iEEG 信号。
-
如果 \rho 远小于 1，则可能是分离失败，需调整算法参数或检查预白化和协方差矩阵的计算。
第三步：脑连接性
这部分的目标是通过相干函数（coherence function）来分析和量化两个神经群之间的功能连接性。相干函数 \rho_{uv}(f) 衡量两个随机过程 u(n) 和 v(n) 在频率 f 上的线性相关性，是功能连接性的一个关键指标。
理论部分：相干函数的性质和证明
两个过程之间相干函数的定义：
- \gamma_{uu}(f) 和 \gamma_{vv}(f) 是过程 u(n) 和 v(n) 的功率谱密度（PSD），表示信号能量在不同频率的分布；
- \gamma_{uv}(f) 是互功率谱密度，描述了两个信号在频域上的相关性。
现在我们要证明，对于两个中心化的随机过程 u(n) 和 v(n) ，如果它们在二阶上是平稳的，那么有：
从"注意到"角度来讲: 因为分母的平方根是各信号能量的平方根之积，分子是它们之间的交互能量，分子不可能超过分母。
数学证明1 :
我们设计设计过程 z(n) 为：
- 其中 \lambda 是一个复数
根据线性性质和互功率谱的定义，有
该过程的频域表示为 Z(f)，则它的功率谱密度为
因此
展开后可得
由于二阶平稳的复共轭对称性，常见约定是
从而上式可写作
功率谱 \gamma_{zz}(f) 不可能为负，即对所有频率 f 有
下面的关键是选择合适的复数 \lambda 使得上式尽可能“紧”，从而给出对 \gamma_{uv}(f) 的强约束。通常做法是让 \lambda 的相位与 \gamma_{uv}(f) 的相位相抵消，使交叉项成为负的实数，从而逼近最小值。
令
其中 r 为实数且 r>0（当 \gamma_{uv}(f)=0 时，|\rho_{uv}(f)|本就等于0，直接成立）。这样一来，
将此代入 \gamma_{zz}(f)，得到
因为
所以交叉项之和为
于是
上式是关于 r 的一个“抛物线”形式。令
若要保证 \Phi(r)\ge 0 对所有实 r 都成立，则其判别式 \Delta 必须非正，即
因此
除以4后得到
由相干函数定义
因此
结合上面得到的
可以立刻推出
另外，由模的定义可知 |\rho_{uv}(f)|^2 \ge 0 始终成立，故最终得证
这就完成了对相干函数的下上界的证明。
数学证明2 : 没有本质差别，就是中间步骤更详细一些，我不舍得删除
我们先计算 z(n) 的功率谱密度并选择：
功率谱密度 \gamma_{zz}(f) 是信号 z(n) 的傅里叶变换后的功率密度分布。对于 z(n) ，它的功率谱密度定义为：
其中 R_{zz}(\tau) 是 z(n) 的自相关函数：
所以我们需要先计算 R_{zz}(\tau) 。
代入 z(n) = u(n) + \lambda v(n) ：
展开后：
下面我们对 z(n)z^*(n+\tau) 取期望，根据线性期望运算的性质：
利用功率谱和互功率谱的定义：
- \mathbb{E}[u(n)u^*(n+\tau)] 对应自相关函数 R_{uu}(\tau)
- \mathbb{E}[v(n)v^*(n+\tau)] 对应自相关函数 R_{vv}(\tau)
- \mathbb{E}[u(n)v^*(n+\tau)] 对应互相关函数 R_{uv}(\tau)
- \mathbb{E}[v(n)u^*(n+\tau)] 对应互相关函数的复共轭 R_{uv}^*(\tau)
因此：
下面我们对自相关函数取傅立叶变换的到功率谱密度:
对 R_{zz}(\tau) 的各项分别取傅里叶变换：
- \mathcal{F}\{R_{uu}(\tau)\} = \gamma_{uu}(f)，得到 u(n) 的功率谱密度
- \mathcal{F}\{R_{vv}(\tau)\} = \gamma_{vv}(f)，得到 v(n) 的功率谱密度
- \mathcal{F}\{R_{uv}(\tau)\} = \gamma_{uv}(f)，得到 u(n) 和 v(n) 的互功率谱密度
- \mathcal{F}\{R_{uv}^*(\tau)\} = \gamma_{uv}^*(f)，得到互功率谱密度的复共轭
因此：
后面做法同方法1。
实验
现在，考虑如下向量自回归模型（参见图 2）对应的两个 iEEG 信号：
其中， w_j, j = 1, 2 是具有零均值和单位方差的独立白噪声的实现。
在公式 (1) 中， y_1(n) 和 y_2(n) 可以视为 iEEG 信号。
实验步骤：
a. 在第一步中，仅将 y_1(n) 视为方程 (1) 中给出的 AR 模型，同时 y_2(n) 为 y_1(n) 的一个延迟版本，即：
计算并绘制两个信号在 4 个长度为 1024 点的块上、50% 重叠情况下的频谱（使用 pwelch 命令）。计算并绘制这两个信号在不同块上的幅度平方相干性（MSC），并得出结论（使用 mscohere 命令）。
对应代码
我们定义AR(2)模型的两个系数a1与a2，其中a1=0.95\sqrt{2}，a2=-0.9025。同时设置样本长度 N。
N = 8192;
a1 = 0.95*sqrt(2);
a2 = -0.9025;
产生均值0、方差1的白噪声w1(n)
w1 = randn(1, N);
接着初始化y1向量，并给定前两个点(起始条件)均为0。
y1 = zeros(1, N);
y1(1) = 0;
y1(2) = 0;
这里使用AR(2)模型的迭代公式生成y1(n)。
for n = 3:N
y1(n) = a1*y1(n-1) + a2*y1(n-2) + w1(n);
end
接下来，我们构造y2(n)=0.5y1(n-1)。
y2 = zeros(1, N);
for n = 2:N
y2(n) = 0.5 * y1(n-1);
end
为了使用 Welch 方法估计功率谱，我们设定分段长度 seg_len=1024以及50%重叠，也就是512点重叠。
seg_len = 1024;
overlap = seg_len/2;
再设置 nfft=seg\_len，并假设采样频率fs=1。
nfft = seg_len;
fs = 1;
使用 pwelch对y1(n)计算 Welch PSD，并将结果存储在 PSD_y1以及频率向量f里。
[PSD_y1, f] = pwelch(y1, hanning(seg_len), overlap, nfft, fs);
同样地，对y2(n)计算 Welch PSD，保存在 PSD_y2中。
PSD_y2 = pwelch(y2, hanning(seg_len), overlap, nfft, fs);
这里开始绘制两个功率谱图像，第一幅图绘制y1(n)的PSD。
figure('Name','Welch PSD of y1 & y2');
subplot(2,1,1);
plot(f, 10*log10(PSD_y1), 'LineWidth', 1.2);
title('Welch PSD of y_1');
xlabel('频率'); ylabel('PSD (dB)');
grid on;
第二幅图绘制y2(n)的PSD，保持相同坐标轴格式。
subplot(2,1,2);
plot(f, 10*log10(PSD_y2), 'LineWidth', 1.2);
title('Welch PSD of y_2');
xlabel('频率'); ylabel('PSD (dB)');
grid on;
接下来，用 mscohere函数计算y1和y2的幅度平方相干性(MSC)，返回结果Cxy和对应频率轴 f_c。
[Cxy, f_c] = mscohere(y1, y2, hanning(seg_len), overlap, nfft, fs);
这里我们将MSC随频率的变化绘制出来。
figure('Name','Magnitude-Squared Coherence (MSC)');
plot(f_c, Cxy, 'LineWidth', 1.2);
title('MSC between y_1 and y_2');
xlabel('频率'); ylabel('Coherence');
grid on;
完整代码及结果
clear; clc; close all;
% ========== 1. 生成 AR(2) 信号 y1(n) ==========
N = 8192; % 样本长度
a1 = 0.95*sqrt(2); % 系数
a2 = -0.9025; % 系数
w1 = randn(1, N); % 白噪声 w1(n), 均值0, 方差1
% 初始值
y1 = zeros(1, N);
y1(1) = 0;
y1(2) = 0;
% 根据 AR 模型生成 y1(n)
for n = 3:N
y1(n) = a1*y1(n-1) + a2*y1(n-2) + w1(n);
end
% ========== 2. 构造 y2(n) = 0.5 * y1(n-1) ==========
y2 = zeros(1, N);
for n = 2:N
y2(n) = 0.5 * y1(n-1);
end
% ========== 3. Welch 功率谱估计 ==========
% 按题意: 分块长度为 1024, 4 个块, 50% (512点) 重叠
seg_len = 1024;
overlap = seg_len/2;
nfft = seg_len; % 同步长, 以便可视化
fs = 1; % 设 fs=1
% (a) y1 的 PSD
[PSD_y1, f] = pwelch(y1, hanning(seg_len), overlap, nfft, fs);
% (b) y2 的 PSD
PSD_y2 = pwelch(y2, hanning(seg_len), overlap, nfft, fs);
% 绘图: y1, y2 的功率谱 (dB)
figure('Name','Welch PSD of y1 & y2');
subplot(2,1,1);
plot(f, 10*log10(PSD_y1), 'LineWidth', 1.2);
title('Welch PSD of y_1');
xlabel('频率'); ylabel('PSD (dB)'); grid on;
subplot(2,1,2);
plot(f, 10*log10(PSD_y2), 'LineWidth', 1.2);
title('Welch PSD of y_2');
xlabel('频率'); ylabel('PSD (dB)'); grid on;
% ========== 4. 计算 y1, y2 的幅度平方相干性 (MSC) ==========
[Cxy, f_c] = mscohere(y1, y2, hanning(seg_len), overlap, nfft, fs);
figure('Name','Magnitude-Squared Coherence (MSC)');
plot(f_c, Cxy, 'LineWidth', 1.2);
title('MSC between y_1 and y_2');
xlabel('频率'); ylabel('Coherence'); grid on;
b. 在第二步中，考虑 y_1(n) 和 y_2(n) 为方程 (1) 中给出的两个信号。计算两个信号在 100 个长度为 1024 点的块上、50% 重叠情况下的频谱，以及 MSC。绘制这些 100 个块上计算得到的平均 MSC，并对结果进行分析。
设定信号总长度N=1024\times100=102400，这样可以在 100 个长度为 1024 的块里做频谱分析，并且采用 50% (512点) 重叠。
N = 102400;
为方程 (1) 中的两个 iEEG 信号设置 AR 系数。方程 (1) 给出：
这里先定义好这些系数。
a1_1 = 0.95*sqrt(2); % y1(n)的一阶系数
a1_2 = -0.9025; % y1(n)的二阶系数
a2_1 = -0.5; % y2(n)关于y1(n-1)的系数
a2_2 = 0.25*sqrt(2); % y2(n)关于y2(n-1)的系数
生成长度为N的独立白噪声w_1(n)、w_2(n)，均为零均值、方差 1，用于激励这两个 AR 模型。
w1 = randn(1, N);
w2 = randn(1, N);
给两个信号y_1(n)和y_2(n)分配空间，并初始化前若干采样值，这里简单设为 0。之后，我们将依次根据方程 (1) 来迭代生成它们。
y1 = zeros(1, N);
y2 = zeros(1, N);
% 初始化
y1(1) = 0; y1(2) = 0;
y2(1) = 0; y2(2) = 0;
使用方程 (1) 的迭代形式生成y_1(n)和y_2(n)。注意：当
for n = 3:N
y1(n) = a1_1 * y1(n-1) + a1_2 * y1(n-2) + w1(n);
y2(n) = a2_1 * y1(n-1) + a2_2 * y2(n-1) + w2(n);
end
设置 Welch 法相关参数：分段长度seg_len= 1024，重叠长度 overlap= 512 (即 50%)，以及 FFT 点数nfft=seg\_len、采样频率fs=1。这样可以确保我们在整个信号上分到 100 个块（前后有部分重叠），可用于做频谱和相干性分析。
seg_len = 1024;
overlap = seg_len/2;
nfft = seg_len;
fs = 1;
利用
PSD_y1, PSD_y2，以及对应的频率向量f。
[PSD_y1, f] = pwelch(y1, hanning(seg_len), overlap, nfft, fs);
PSD_y2 = pwelch(y2, hanning(seg_len), overlap, nfft, fs);
绘制y_1和y_2的 Welch 功率谱，观察它们的幅度随频率分布 (以 dB 标度绘图)。
figure('Name','Welch PSD of y1 & y2');
subplot(2,1,1);
plot(f, 10*log10(PSD_y1), 'LineWidth',1.2);
title('Welch PSD of y_1 (Equation (1))');
xlabel('Frequency'); ylabel('PSD (dB)');
grid on;
subplot(2,1,2);
plot(f, 10*log10(PSD_y2), 'LineWidth',1.2);
title('Welch PSD of y_2 (Equation (1))');
xlabel('Frequency'); ylabel('PSD (dB)');
grid on;
使用\text{mscohere}(y1,y2,\dots)计算幅度平方相干性 (MSC)。该函数会在 100 个块的基础上，把每个块的相干性结果平均起来，从而得到最终的相干函数。
[Cxy, f_c] = mscohere(y1, y2, hanning(seg_len), overlap, nfft, fs);
绘制上面得到的 MSC 曲线，横轴是频率f_c，纵轴是 MSC 值 (介于 0 和 1 之间)。
figure('Name','Magnitude-Squared Coherence (MSC)');
plot(f_c, Cxy, 'LineWidth',1.2);
title('MSC between y_1 and y_2 (Equation (1))');
xlabel('Frequency'); ylabel('Coherence');
grid on;
完整代码及结果
N = 102400;
a1_1 = 0.95*sqrt(2); % y1(n)的一阶系数
a1_2 = -0.9025; % y1(n)的二阶系数
a2_1 = -0.5; % y2(n)关于y1(n-1)的系数
a2_2 = 0.25*sqrt(2); % y2(n)关于y2(n-1)的系数
w1 = randn(1, N); w2 = randn(1, N);
y1 = zeros(1, N); y2 = zeros(1, N);
% 初始化
y1(1) = 0; y1(2) = 0;
y2(1) = 0; y2(2) = 0;
for n = 3:N
y1(n) = a1_1 * y1(n-1) + a1_2 * y1(n-2) + w1(n);
y2(n) = a2_1 * y1(n-1) + a2_2 * y2(n-1) + w2(n);
end
seg_len = 1024;
overlap = seg_len/2;
nfft = seg_len;
fs = 1; % 若无明确采样率信息，可先假设为 1
[PSD_y1, f] = pwelch(y1, hanning(seg_len), overlap, nfft, fs);
PSD_y2 = pwelch(y2, hanning(seg_len), overlap, nfft, fs);
figure('Name','Welch PSD of y1 & y2');
subplot(2,1,1);
plot(f, 10*log10(PSD_y1), 'LineWidth',1.2);
title('Welch PSD of y_1 (Equation (1))');
xlabel('Frequency'); ylabel('PSD (dB)');
grid on;
subplot(2,1,2);
plot(f, 10*log10(PSD_y2), 'LineWidth',1.2);
title('Welch PSD of y_2 (Equation (1))');
xlabel('Frequency'); ylabel('PSD (dB)');
grid on;
[Cxy, f_c] = mscohere(y1, y2, hanning(seg_len), overlap, nfft, fs);
figure('Name','Magnitude-Squared Coherence (MSC)');
plot(f_c, Cxy, 'LineWidth',1.2);
title('MSC between y_1 and y_2 (Equation (1))');
xlabel('Frequency'); ylabel('Coherence');
grid on;
c. 接下来，引入由生理模型生成的两个信号（提供于 Matlab 文件 ScenarioK1500second1000a.mat 中，并使用变量 LFP_s）。估计这两个信号随时间变化的 MSC。绘制信号的频谱以及 MSC。