基于盲源分离和半盲源分离的心电信号伪影消除方法（MATLAB 2018）

心电信号是通过测量放置在人体皮肤上的电极之间的电位差来获取的，其本身具有信号微弱、频段低、不稳定等特性。因此ECG信号在实际采集时极易受到不同噪声的影响，这会造成心电图本身的波形形态特征的失真，从而导致错误诊断和对患者的不当治疗。由于信号采集的环境、方法和设备等因素的干扰，影响心电信号质量的噪声主要包括：基线漂移、工频干扰、肌电噪声以及电极运动伪影。

基线漂移是一种低频噪声，主要由人体呼吸、身体运动、电极接触不良和皮肤与电极之间阻抗变化引起。漂移的幅度和持续时间取决于电极特性、电解质特性、皮肤阻抗和身体运动，通常基线漂移会造成心电信号ST段的失真，从而导致对心肌梗塞、Brugada综合征等与ST段异常相关的错误诊断。

工频干扰是一种窄带干扰噪声，主要由心电信号采集过程中50/60Hz电源线的电感和电容耦合引起。工频干扰与ECG信号混合会扭曲信号的形态，导致P波失真，从而导致对心房扩大和房颤的错误诊断。

肌电噪声是由肌肉在收缩期间的电活动或由于突然的身体运动引起的，其带宽通常在5-2000Hz之间，平均幅度为心电信号峰-峰值10%的水平，持续时间约为50ms，因此会导致ECG信号的局部波形失真。

电极运动伪影是由电极与皮肤之间阻抗随着电极运动变化引起的瞬时基线变化，其幅度能达到ECG信号峰-峰值的500%，持续时间约为300-500ms。电极运动伪影会严重干扰可穿戴式心电设备的心电信号采集，其较大的波动性和随机性会严重影响可穿戴设备采集到的心电信号质量。

鉴于此，采用盲源分离和半盲源分离方法对心电信号伪影消除，运行环境为MATLAB 2018。

clc;clear all;close all;
X_load=load("Ex1.mat");
X_load=X_load.X;
%% part a
X =X_load';
for dim =1:3
    X(:,dim) = X(:,dim) - sum(X(:,dim))/(size(X,1)); 
end
figure()
subplot (3,3,[1:6])
scatter3(X(:,1),X(:,2),X(:,3))
title("Visualized Data in 3D domain")
xlabel('X'); ylabel("Y"); zlabel("Z");
subplot(3,3,7)
scatter(X(:,1),X(:,2))
title("Visualized Data in 2D x-y plane")
xlabel('X'); ylabel("Y");
grid on
subplot(3,3,8)
scatter(X(:,1),X(:,3))
title("Visualized Data in 2D x-z plane")
xlabel('X'); ylabel("Z");
grid on
subplot(3,3,9)
scatter(X(:,2),X(:,3))
title("Visualized Data in 2D y-z plane")
xlabel('Y'); ylabel("Z");
grid on
%% part B.1 
Covariance_mat = cov(X)
[V,D] = eig(Covariance_mat)
D_wh = D^(-1/2);
C_final = D_wh*V'*Covariance_mat*V*D_wh;
kernel=-1:0.001:1;
figure()
subplot (3,3,[1:6])
scatter3(X(:,1),X(:,2),X(:,3));
hold on
plot3(kernel*V(1,2),kernel*V(2,2),kernel*V(3,2),'linewidth',2.5)
hold on
plot3(kernel*V(1,1),kernel*V(2,1),kernel*V(3,1),'linewidth',2.5)
hold on
plot3(40*kernel*V(1,3),40*kernel*V(2,3),40*kernel*V(3,3),'linewidth',2.5)
title("Visualized Data with eigen vectors")
xlabel("X"); ylabel("Y"); zlabel("Z");
subplot(3,3,7)
scatter3(X(:,1),X(:,2),X(:,3));
hold on
plot3(kernel*V(1,2),kernel*V(2,2),kernel*V(3,2),'linewidth',2.5)
hold on
plot3(kernel*V(1,1),kernel*V(2,1),kernel*V(3,1),'linewidth',2.5)
hold on
plot3(40*kernel*V(1,3),40*kernel*V(2,3),40*kernel*V(3,3),'linewidth',2.5)
title("Visualized Data with eigen vectors")
xlabel("X"); ylabel("Y"); zlabel("Z");
subplot(3,3,8)
scatter3(X(:,1),X(:,2),X(:,3));
hold on
plot3(kernel*V(1,2),kernel*V(2,2),kernel*V(3,2),'linewidth',2.5)
hold on
plot3(kernel*V(1,1),kernel*V(2,1),kernel*V(3,1),'linewidth',2.5)
hold on
plot3(40*kernel*V(1,3),40*kernel*V(2,3),40*kernel*V(3,3),'linewidth',2.5)
title("Visualized Data with eigen vectors")
xlabel("X"); ylabel("Y"); zlabel("Z");
subplot(3,3,9)
scatter3(X(:,1),X(:,2),X(:,3));
hold on
plot3(kernel*V(1,2),kernel*V(2,2),kernel*V(3,2),'linewidth',2.5)
hold on
plot3(kernel*V(1,1),kernel*V(2,1),kernel*V(3,1),'linewidth',2.5)
hold on
plot3(40*kernel*V(1,3),40*kernel*V(2,3),40*kernel*V(3,3),'linewidth',2.5)
title("Visualized Data with eigen vectors")
xlabel("X"); ylabel("Y"); zlabel("Z");
%% part B.2
Y = transpose(D_wh*V'*X'); 
Covariance_matY =cov(Y)
figure();
scatter3(Y(:,1),Y(:,2),Y(:,3));
title('Visualization of whitened Data');
xlabel('X');ylabel('Y');zlabel('Z'); 
%% Part C.1
COEF_pca = pca(X);
figure()
scatter3(X(:,1),X(:,2),X(:,3));
hold on
plot3(4*kernel*COEF_pca(1,2),kernel*COEF_pca(2,2),kernel*COEF_pca(3,2),'linewidth',2)
hold on
plot3(4*kernel*COEF_pca(1,1),kernel*COEF_pca(2,1),kernel*COEF_pca(3,1),'linewidth',2)
hold on
plot3(40*kernel*COEF_pca(1,3),40*kernel*COEF_pca(2,3),40*kernel*COEF_pca(3,3),'linewidth',2)
title('Visualized Data with eigen vectors (matlab)');
xlabel('X');ylabel('Y');zlabel('Z');
%% Part C.2
figure()
Y_matlab = COEF_pca'*X';
scatter3(Y_matlab(1,:),Y_matlab(2,:),Y_matlab(3,:));
title('Visuallization of whitened Data (using matlab pca)');
xlabel('X');ylabel('Y');zlabel('Z');
covarianceY_matlab =cov(Y_matlab')
%%
%完整代码可通过知乎学术咨询获得：https://www.zhihu.com/consult/people/792359672131756032?isMe=1
[U_svd,S_svd,V_svd] = svd(X);