【信号分析之心电信号处理】

信号分析之心电信号处理

任选下面的一组数据，利用自适应滤波中LMS和RLS方法进行处理，分析不同阶数、步长及指数加权因子对滤波结果的影响。同时给出迭代次数与滤波器系数，迭代次数与均方误差之间的关系曲线。（20分）
（1）消除心电图的电源干扰（数据源∶x1子目录）
（序号为偶数）
（2）胎儿心电监护（数据源∶ x2子目录）
（序号为奇数）
解答：胎儿心电监护

图 4.1 胎儿心电监护原理图
胎儿心电监护原理：监测胎儿心电，需要将母亲的心音及背景干扰去除。母亲的心音强度通常是胎儿心音的 2-10 倍，肌肉活动及胎儿运动产生的背景干扰 也常常大于胎儿心音的强度。所以采用自适应对消的方法来监护胎儿心电，将母亲胸部的信号作为参考信号，主要包括母亲的胸部心音和背景噪声，母亲腹部取出的信号作为输入信号，包括母亲、胎儿心音和背景干扰。
系统是一个有信号分量泄露到参考输入端的自适应对消系统。希望的结果是对消后母亲的心电信号的强度小于胎儿的心电信号强度，且胎儿的心电信号要比较清楚。利用 Matlab 进行数据分析，首先将原始心电数据导入Matlab并绘制如下图所示的心电信号波形图。

绘制母亲心电信号Matlab 代码：
 close all; % 清除工作空间的所有变量，函数，和MEX文件
clear all;% 加载数据文件，并命名为A

chest=load('chestt.txt');%读取胸部数据
abdomen=load('abdomenn.txt');%读取腹部数据
test_chest = chest.'; 
test_abdomen = abdomen.'; 

% 绘制txt文件第一列的数据
subplot(211);
plot(test_chest );
% 横坐标
xlabel('采样点数(个)');
% 纵坐标
ylabel('幅值(A)');
% 标题
title('胸部心电波形图');subplot(212);plot(test_abdomen );
% 横坐标
xlabel('采样点数(个)');
% 纵坐标
ylabel('幅值(A)');
% 标题
title('腹部心电波形图');

图 4.2 心电图

1）LMS滤波器

最小均方误差(LMS)算法是由Widrow和Hoff提出的，具有计算量小、易于实现等优点在实践中广泛应用。LMS算法的基本思想：调整滤波器自身的参数，使滤波器的输出信号与期望输出信号之间的均方误差最小，系统输出为有用信号的最佳估计。实质上LMS可以看成是一种随机梯度或者随机逼近算法，可以写成如下的基本迭代方程：

其中μ为步长因子，是控制稳定性和收敛速度的参量，X(k)表示第k时刻参考信号矢量，，k为迭代次数，M为滤波器的阶数。d(k)表示第k时刻的输入信号矢量，y(k)、e(k)分别表示第k时刻的输出信号与输出误差，W(k)表示k时刻权系数矢量，W(k)=[W(k，0)，W(k，1)…W(k，M-1)]。从上式可以看出其结构简单、计算量小且稳定性好等优点，但固定步长的LMS算法在收敛速度、跟踪速率及权失调噪声之间的要求是相互矛盾的，为了克服其缺点，人们提出了各种变步长的LMS改进算法，主要是采用减小均方误差或者以某种规则基于时变步长因子跟踪信号的时变，其中有正规LMS算法(NLMS)、梯度自适应步长算法、自动增益控制自适应算法、符号-误差LMS算法、符号-数据LMS算法、数据复用LMS算法等。
自适应滤波器收敛的条件是。其中是输入信号的自相关矩阵R的最大特征值。μ的选取必须在收敛速度和失调之间取得较好的折中，既要具有较快的收敛速度，又要使稳态误差最小。它控制了算法稳定性和自适应速度，如果μ很小，算法的自适应速度会很慢；如果μ很大，算法会变得不稳定。
为了分析不同阶数对滤波结果的影响，设置滤波器阶数分别为 10、20、
30，同时步长固定设置为 0.05，对数据进行处理分析。在画图分析中仅画出前三个系数与迭代次数的关系曲线图。
LMS算法Matlab代码：

%%%LMS滤波器%%%%%%%%%%%%%%%
%%2022年6月20日%%
close all;
chest=load('chestt.txt');%读取胸部数据数据
abdomen=load('abdomenn.txt');%读取腹部数据数据
test_chest = chest'; %转置
test_abdomen = abdomen'; %转置
M = 10; %阶数可改为10,20,40 
itr = 2000; %迭代次数2000 
en = zeros(itr,1); % 误差序列初值
W = zeros(M,itr); % 权系数 
mu = 0.05; % 步长因子可更改 
s = zeros(itr,1); % 均方误差初值 
% LMS 迭代 
for k = M:itr 
x = test_chest(k:-1:k-M+1); %胸部信号为参考信号
y = W(:,k-1).' * x;  %滤波处理后胸部信号
en(k) = test_abdomen(k) - y; %腹部信号减去滤波处理后胸部信号
W(:,k) = W(:,k-1) + 2 * mu * en(k) * x; %迭代
s(k) = s(k-1) + abs(en(k)).^2; %误差迭代
c(k) = (s(k)/k); 
end 
yn = inf * ones(size(test_chest)); 
for k = M:length(test_chest) 
x= test_chest(k:-1:k-M+1); 
yn(k) = W(:,end).' * x; 
end 
subplot(311); 
plot((1:2000),test_abdomen-yn,'r',(1:2000),test_chest,'g'); 
legend('胎儿心电图','母亲心电图');xlabel('采样点数');ylabel('幅度'); 
subplot(312); 
plot(c); 
xlabel('迭代次数');ylabel('均方误差'); 
subplot(313); 
plot(W(1,:),'r');hold on;
plot(W(2,:));
plot(W(3,:),'k'); 
xlabel('迭代次数');ylabel('滤波器系数'); 
legend('权限系 1','权限系 2','权限系 3');grid on;

当阶数M为10时：

图 4.3 a) 阶数为 10，步长为 0.05 时的胎儿心电图 b) 阶数为 10，步长为 0.05 时采样点数与均方误差关系曲线 c) 阶数为 10，步长为 0.05 时采样点数与滤波器系数关系曲线
当阶数M为20时：

图 4.4 a) 阶数为 20，步长为 0.05 时的心电图 b) 阶数为 20，步长为 0.05 时采样点数与均方误差关系曲线 c) 阶数为 20，步长为 0.05 时采样点数与滤波器系数关系曲线
当阶数M为40时：

图 4.5 a)阶数为 40，步长为 0.05 时的心电图 b) 阶数为 40，步长为 0.05 时采样点数与均方误差关系曲线 c) 阶数为40，步长为 0.05 时采样点数与滤波器系数关系曲线
由以上图可看出，当阶数增加时，胎儿心电信号幅值有一定增大，均方误差收敛速度变快，但是前期收敛稳定性较差，因此自适应滤波器的阶数取值并不是越大越好，不可任意选取,而是需要根据经验和实际仿真验证比较才能最终确定。当阶数大到一定程度时，收敛速度变化不明显了，而且还有可能引起系数迭代过程不收敛。
下面分析不同步长因子对滤波结果的影响，我们选取滤波器阶数固定为20 阶，同时步长分别设置为 0.05、0.01、0.1，对数据进行处理分析，步长为0.05时，图4.4已给出，其他两种情况如下图所示：
当步长为0.01时：

图 4.6 a) 阶数为 20，步长为 0.01 时的心电图 b) 阶数为 20，步长为 0.01 时采样点数与均方误差关系曲线 c) 阶数为20，步长为 0.01 时采样点数与滤波器系数关系曲线
当步长为0.1时：

图 4.7 a) 阶数为 20，步长为 0.1 时的心电图 b) 阶数为 20，步长为 0.1 时采样点数与均方误差关系曲线 c) 阶数为20，步长为 0.1 时采样点数与滤波器系数关系曲线
根据以上分析可以看出，当步长为0.1时，胎儿心电幅值更大，更为明显，步长影响系统的收敛速度和均方误差曲线的超调量。步长越大时从均方误差曲线图可以看出系统收敛速度变快，但是对应均方误差曲线超调量变大。另外步长越大，滤波性能越差。实际应用中应结合精度要求和收敛时间要求设置步长。

2）RLS滤波器

RLS（Recursive Least-Square）算法是通过矩阵求逆引理而导出的，它对输入信号的自相关矩阵的逆矩阵进行递推估计更新，其迭代公式如下：

上式中，n是递推次数，输入为x(n)，P(n)为信号自相关矩阵的逆，g(n)为增益矢量，误差为e(n)，期望输出为d(n)，权重因子为λ。RLS算法能实现快速收敛，且其收敛速率不随输入向量x(n)平均相关矩阵R特征值的扩散度而改变。当工作于时变环境中时，这类算法具有极好的性能，但都以增加计算复杂度和稳定性为代价。
Matlab仿真分析：滤波器的阶数为20，初始化设置为P(0)=δ-1I，其中ߜδ=0.001，w(0)=wi=0，同时改变遗忘因子la，对la分别取 0.1，0.5，0.75、0.9，分别观察滤波结果。同 LMS 算法时一样，由于滤波器系数过多，在画图分析时只画出前三个滤波器的系数与迭代次数的曲线。分析结果如下图，可修改la值。

```c
RLS滤波Matlab代码：
%%%RLS滤波器%%%
%%2022.6.21%%
close all;
chest=load('chestt.txt');%读取胸部数据
abdomen=load('abdomenn.txt');%读取腹部数据
test_chest = chest.'; 
test_abdomen = abdomen.'; 

M = 20; %%%阶数
itr = 2000; 
en = zeros(itr,1); % 误差序列 
W = zeros(M,itr); % 权系数 
s = zeros(itr,1); % 均方误差 
la = 0.5; % %%可选取不同遗忘因子0.1,0.5,0.9 
Rxx = eye(M) * 0.001; 
for k = M:itr 
x = test_chest(k:-1:k-M+1);y(k) = W(:,k).' * x;  
K(:,k) = (Rxx * x)./(la + x.' * Rxx * x); 
e(k) = test_abdomen(k) - W(:,k-1).' * x; 
W(:,k) = W(:,k-1)+K(:,k) * e(k); %%权系数更新
Rxx = (Rxx - K(:,k) * x.' * Rxx)/la;
s(k) = sum(abs(e).^2); %均方误差
c(k) = (s(k)/k); 
end 

```c
```c

subplot(311); 
plot((1:2000),test_

abdomen - y.','r',(1:2000),test_chest,'g'); 
legend('胎儿心电图','母亲心电图');xlabel('采样点数');ylabel('幅度'); grid on;
subplot(312); 
plot(c); 
xlabel('迭代次数');ylabel('均方误差'); grid on;
subplot(313); 
plot(W(1,:),'r');hold on;plot(W(2,:));plot(W(3,:),'k');grid on;

采用RLS滤波器分析如下所示，遗忘因子为la=0.5。

图 4.4  a) 遗忘因子为la=0.5心电图  b)采样点数与均方误差关系   c）采样点数与滤波器权系数关系
下面分析不同遗忘因子对滤波的影响：
遗忘因子为la=0.1时：

图 4.9  a) 遗忘因子为0.1时的心电波形  b) 遗忘因子为0.1时均方误差与迭代次数关系
c) 遗忘因子为0.1时系数与迭代次数关系
当遗忘因子为la=0.75时：

图 4.9 遗忘因子为0.75时心电图  b) 遗忘因子为0.75时均方差与迭代次数关系
c) 遗忘因子为0.75时系数与迭代次数关系
当遗忘因子为la=0.9时：



图 4.10 遗忘因子为0.9时心电图  b) 遗忘因子为0.9时均方差与迭代次数关系
c) 遗忘因子为0.9时系数与迭代次数关系
根据以上仿真分析可看出，当遗忘因子为0.1、0.5、0.75和0.9时，胎儿心电信号幅值变化不大，都在0.2左右，但均方误差曲线不同。当遗忘因子在 0.9 时，均方误差逐渐收敛，为0.1时明显波动较大。总之，遗忘因子越小，系统的跟踪能力越强，同时对噪声越敏感；其值越大，系统跟踪能力减弱，但对噪声不敏感，收敛时估计误差也越小。根据经验，采用RLS算法时，一定要对la有个准确的取值，才能保证系统性能的最佳状态，所以一般文献要求la的取值范围为0.95-0.995。
总结：相比 LMS 算法和 RLS 算法，RLS 算法在迭代过程中产生的误差明显小于 LMS 算法，由此可以得到 RLS 在提取信号时，收敛速度快，估计精度高而且稳定性好，而 LMS 算法收敛速度慢，估计精度低而且权系数估计值因瞬时梯度估计围绕精确值波动较大，权噪声大，稳定性差。

```bash