clear all ;
load foetal_ecg.dat
abdominal =foetal_ecg(:,2:6); %loading abdominal signals
thoraic =foetal_ecg(:,7:9) ; %loading thoraic signals
d= mean(abdominal,2); %Average of abdominal signals
x=mean(thoraic,2); %Average of thoraic signals
M=2500; %数据长度
N = 8; %order of the filter
mu = .0000006; %value of step size
iter=length(x); %迭代次数默认为输入信号长度
w=zeros(N,M); %每一行代表一个加权参量，每一列代表一次迭代，初始为 0
e=zeros(M,1); %误差序列， e(k) 表示第 k 次迭代时预期输出与实际输入的误差
%将 y 信号右移 N位，以便进行迭代时的加权运算
x1=zeros(N,1);
for i=N+1:(N+M-1)
x1(i)=x(i+1-N);
end
%调整滤波器系数的 LMS算法
for n=2:iter;
x2=x1(n:1:n+N-1);
y(n)=w(:,n-1).'*x2; %每次取 N个参考信号与权系数相乘得到噪声估计
e(n)=d(n)-y(n); %噪声抵消的输出
w(:,n)=w(:,n-1)+2*mu*x2*e(n); %权矢量迭代
end ;
%最后稳定后的胎儿心电信号
for n=2:iter;
x2=x1(n:1:n+N-1);
y(n)=w(:,end).'*x2; %每次取 N个参考信号与最后稳定权系数相乘得到噪声估计
e(n)=d(n)-y(n); %噪声抵消的输出
end
%输出三个心电波形
figure(1);
subplot(3,1,1);
plot(x);
title( ' 母亲心电信号 ' );
subplot(3,1,2);
plot(d);
title( ' 混合心电信号 ' );
subplot(3,1,3)
plot(e);
title( ' 滤波后的胎儿心电信号 ' );
