【心电信号】基于matlab心电信号去除基线漂移【含Matlab源码 955期】
一、简介
1 心音：
心脏收缩舒张时产生的声音，可用耳或听诊器在胸壁听到，亦可用电子仪器记录下来（心音图）。可分为第一心音（S1）第二心音（S2）。（正常情况下均可听到）。第三心音(S3通常仅在儿童及青少年可听到)，第四心音（S4正常情况很少听到）。从心脏产生的心音经过组织的介导传到胸壁表面，其中以骨传导最好。
心音是心脏及心血管系统机械运动状况的反映，其中包含着心脏各个部分本身及相互之间作用的生理和病理信息。心音信号的识别与分类对心血管系统疾病的诊断具有重要的意义，其准确性、可靠性的好坏决定着诊断与治疗心脏病患者的效果。早期的心音识别与分类是医生根据听诊结果来完成的，显然这一过程具有一定的主观性且可靠性不高。随着信号处理与分析技术的不断发展，对心音的研究也逐步由定性分析进入了定量分析的阶段。国内外许多生物医学工程研究人员将传统的模式识别方法，以及神经网络方法用于心音的识别与分类，期望实现心音的自动解释和自动诊断，以便向临床医生提供实用的辅助诊断信息。此外，心音的识别与分类还有助于对心音产生机制的理解。
心音信号研究主要是采用微电子技术，检测技术，现代数字信号处理技术和生物医学工程技术，研究和揭示心音与心脏病之间的关系。
无论是多大单位的记录，人的体重，心率，血压等生理参数，都是时变的，称为心率或者血压的变异性。
心音是在心动周期中，由于心肌收缩和舒张，瓣膜启闭，血流冲击心室壁和大 动脉等因素引起的机械振动，通过周围组织传到胸壁，将耳紧贴胸壁或将听诊器放在胸壁 定部位，听到的声音 通常很容易听到第一和第二心音，有时在某些情况下听到第三或第四心音。
第一心音：S1,发生在心脏收缩期开始，音调低沉，持续时间较长(约0.15秒)。产生的原因包括心室肌的收缩，房室瓣突然关闭以及随后射血入主动脉等引起的振动。第一心音的最佳听诊部位在锁骨中线第五肋间隙或在胸骨右缘。相对于心电图上QRS波开始后0.020.04s，占时0.080.15s
第二心音：S2,发生在心脏舒张期的开始，频率较高，持续时间较短(约0.08秒)。产生的原因是半月瓣关闭，瓣膜互相撞击以及大动脉中血液减速和室内压迅速下降引起的振动。第二心音的最佳听诊部位在第二肋间隙右侧的主动脉瓣区和左侧的肺动脉瓣区。相对于T波终末部。
第三心音和第四心音：
第三心音S3发生在第二心音后0.1～0.2秒，频率低，它的产生与血液快速流入心室使心室和瓣膜发生振动有关，通常仅在儿童能听到，因为 较易传导到体表。相当于T波后距第二心音0.12~0.20s。
第四心音S4由心房收缩引起，也称心房音,相当于心电图上P波后0.15~0.18s，振幅低。
心音杂音对正常心音形成了一定的干扰，但心音杂音的出现对心音信号分析具有实际应用价值和临床意义。
根据杂音出现的时间与S1,S2心音的关系，可分为早，中，晚期杂音。杂音的强度一般可视其振幅与S1比较分类。大于S1,强。小于S1,大于1/3S1，中。低于S1的1/3则为低，仅有轻微振动则为极低振幅杂音。
由于心音信号属于强噪声背景下的人体微弱生物信号，由于心音信号是由复杂的生命体发出的不稳定的自然信号。心音的改变和杂音的出现往往是心脏产生器质性病变的早期症状，心脏内部的物理结构发生变化将直接影响和改变心音信号。
目前广泛应用 的心电图检查是心脏变时性和变传导性的最佳检测方法，但不能用来检测心脏的变力性先天心脏瓣膜受损。心脏传导组织病变引起的心脏机械活动障碍不会首先反应在心电图上，却能首先反应在心音信号上。当冠心病的阻塞达到70%以上时才能引起心电图信号的改变，实际上，达到25%就可以改变心音信号。
2 异常心音
包括S2、S2的异常及收缩期，舒张期的附加音（或额外音）。
① 第一心音异常。指S1增强，减弱或分裂。估计S1的响度，最好是用听诊，心音图判断能力有限。S1增强、减弱或强弱不等的临床情况见表1。S1分裂指M1与T1相距>0.04sec，可见于正常儿童、青年及体瘦者，无重要意义。S1异常分裂时M1与T1相距可大于0.06秒，见于电激动延迟(如右束支传导阻滞等)、机械活动延迟（如房间隔缺损、严重二尖瓣狭窄等。听诊S1分裂在二尖瓣及三尖瓣区最清楚坐位及呼气时更清楚。
② 第二心音异常。包括S2增强，减弱或分裂。
S2增强又分P2增强及A2增强，P2增强见于肺血流量增多（如间隔缺损），肺血管阻力增加；肺静脉压力增高(如二尖瓣狭窄)。P2亢进一般在肺动脉瓣区听到。A2增强见于体循环阻力增高或血流量增多，向肺动脉瓣及心尖区传导，见于 高血压等。
S2减弱又分P2减弱及A2减弱。P2减弱见于肺动脉压低、肺血流量减少或 肺动脉瓣狭窄等。A2减弱见于体循环阻力低、血流减少、血压低、 主动脉瓣狭窄或严重关闭不全。
S2分裂可为生理性。右室射血结束稍晚于左室，吸气时P2延迟出现，此时可听到S2分裂。但呼气时A2与P2接近或重叠，分裂消失。这见于青少年，在肺动脉瓣区听诊明显，坐位时可消失。
S2异常分裂包括宽分裂、固定分裂、逆分裂、分裂减窄。宽分裂是呼气时S2分裂不消失，见于右室射血时间延长或左室射血时间缩短等情况。S2固定分裂指呼吸时A2 -P2间期无明显改变或变动<0.02sec，见于分流量较大的房间隔缺损等。S2逆分裂指A2在P2之后，吸气时A2-P2分裂不显，呼气时P2提早出现，分裂增宽。见于主动脉瓣关闭延迟。S2分裂减窄常由于严重肺动脉高压P2较早出现所致。
3 心音信号分析方法：
传统的谱分析方法通过快速傅立叶变换将时，频域关联起来。但FFT时频域分离，并以信号的频率特性时不变，或统计特性稳定为前提。
传统的稳态分析方法反映的是信号的静态频谱特征，对于包括人体心音信号在内的生物学生理信号，由于环境的影响而表现为非平稳时变特性。因此采用经典谱分析方法难以准确反映心音信号的动态变化。
在传统心音分析的基础上，提出了很多方法：
1.短时傅立叶变换（STFT）
2.小波变换
3.其他时频分析方法
二、源代码
clc;
clear;
close all;
%% 提取信号
M = importdata(‘3.txt‘);
fsample=1000;%采样率为1KHz
[mx,my]=size(M);
Signal=M(:,2);%M的第一列为时间，第二列为信号
length=floor(mx/2);%取原始信号的一半。
S=Signal(1:length);
%% 高通滤波，去除基线漂移的影响
disp(‘-------------------------------------------‘);
disp(‘1：工具箱巴特沃斯高通滤波器‘);
disp(‘2：IIR高通滤波‘);
disp(‘3：FIR高通滤波‘);
disp(‘4：中值滤波‘);
disp(‘5：稀疏小波滤波‘);
disp(‘6：中值+小波滤波‘);
disp(‘-------------------------------------------‘);
choose=input(‘选择滤波方式choose=‘);
function [s, err_mse, iter_time]=sparseOMP(x,A,m,varargin)
[n1 n2]=size(x);
if n2 == 1
n=n1;
elseif n1 == 1
x=x‘;
n=n2;
else
error(‘x must be a vector.‘);
end
sigsize = x‘*x/n;
initial_given=0;
err_mse = [];
iter_time = [];
STOPCRIT = ‘M‘;
STOPTOL = ceil(n/4);
MAXITER = n;
verbose = false;
s_initial = zeros(m,1);
GradSteps = ‘auto‘;
alpha = 1;
weakness = 1;
if verbose
display(‘Initialising...‘)
end
switch nargout
case 3
comp_err=true;
comp_time=true;
case 2
comp_err=true;
comp_time=false;
case 1
comp_err=false;
comp_time=false;
case 0
error(‘Please assign output variable.‘)
otherwise
error(‘Too many output arguments specified‘)
end
% Put option into nice format
Options={};
OS=nargin-3;
c=1;
for i=1:OS
if isa(varargin{i},‘cell‘)
CellSize=length(varargin{i});
ThisCell=varargin{i};
for j=1:CellSize
Options{c}=ThisCell{j};
c=c+1;
end
else
Options{c}=varargin{i};
c=c+1;
end
end
OS=length(Options);
if rem(OS,2)
error(‘Something is wrong with argument name and argument value pairs.‘)
end
for i=1:2:OS
switch Options{i}
case {‘stopCrit‘}
if (strmatch(Options{i+1},{‘M‘; ‘corr‘; ‘mse‘; ‘mse_change‘},‘exact‘));
STOPCRIT = Options{i+1};
else error(‘stopCrit must be char string [M, corr, mse, mse_change]. Exiting.‘); end
case {‘stopTol‘}
if isa(Options{i+1},‘numeric‘) ; STOPTOL = Options{i+1};
else error(‘stopTol must be number. Exiting.‘); end
case {‘P_trans‘}
if isa(Options{i+1},‘function_handle‘); Pt = Options{i+1};
else error(‘P_trans must be function _handle. Exiting.‘); end
case {‘maxIter‘}
if isa(Options{i+1},‘numeric‘); MAXITER = Options{i+1};
else error(‘maxIter must be a number. Exiting.‘); end
case {‘wf‘}
if isa(Options{i+1},‘numeric‘); alpha = Options{i+1};
if alpha <1 weakness =0; else alpha =1; weakness = 1; end
else error(‘wf must be a number. Exiting.‘); end
case {‘verbose‘}
if isa(Options{i+1},‘logical‘); verbose = Options{i+1};
else error(‘verbose must be a logical. Exiting.‘); end
case {‘start_val‘}
if isa(Options{i+1},‘numeric‘) && length(Options{i+1}) == m ;
s_initial = Options{i+1};
initial_given=1;
else error(‘start_val must be a vector of length m. Exiting.‘); end
case {‘GradSteps‘}
if isa(Options{i+1},‘numeric‘) || strcmp(Options{i+1},‘auto‘) ;
GradSteps = Options{i+1};
else error(‘start_val must be a vector of length m. Exiting.‘); end
otherwise
error(‘Unrecognised option. Exiting.‘)
end
end
if strcmp(STOPCRIT,‘M‘)
maxM=STOPTOL;
else
maxM=MAXITER;
end
if nargout >=2
err_mse = zeros(maxM,1);
end
if nargout ==3
iter_time = zeros(maxM,1);
end
if isa(A,‘float‘) P =@(z) A*z; Pt =@(z) A‘*z;
elseif isobject(A) P =@(z) A*z; Pt =@(z) A‘*z;
elseif isa(A,‘function_handle‘)
try
if isa(Pt,‘function_handle‘); P=A;
else error(‘If P is a function handle, Pt also needs to be a function handle. Exiting.‘); end
catch error(‘If P is a function handle, Pt needs to be specified. Exiting.‘); end
else error(‘P is of unsupported type. Use matrix, function_handle or object. Exiting.‘); end
if initial_given ==1;
IN = find(s_initial);
Residual = x-P(s_initial);
s = s_initial;
oldERR = Residual‘*Residual/n;
else
IN = [];
Residual = x;
s = s_initial;
sigsize = x‘*x/n;
oldERR = sigsize;
end
mask=zeros(m,1);
mask(ceil(rand*m))=1;
nP=norm(P(mask));
if abs(1-nP)>1e-3;
display(‘Dictionary appears not to have unit norm columns.‘)
end
if verbose
display(‘Main iterations...‘)
end
tic
t=0;
normA=(sum(A.^2,1)).^0.5;
NI = 1:size(A,2);
p=zeros(m,1);
DR=Pt(Residual);
[v I]=max(abs(DR(NI))./normA(NI)‘);
INI = NI(I);
IN=[IN INI];
NI(I) = [];
if weakness ~= 1
[vals inds] = sort(abs(DR),‘descend‘);
I=inds( find( vals >= alpha * v ) );
end
IN = union(IN,I);
if strcmp(STOPCRIT,‘M‘) & length(IN) >= STOPTOL
IN=IN(1:STOPTOL);
end
MASK=zeros(size(DR));
pDDp=1;
done = 0;
iter=1;
while ~done
% Select new element
if isa(GradSteps,‘char‘)
if strcmp(GradSteps,‘auto‘)
% finished=0;
% while ~finished
% Update direction
if iter==1
p(IN)=DR(IN);
Dp=P(p);
else
MASK(IN)=1;
PDR=P(DR.*MASK);
b=-Dp‘*PDR/pDDp;
p(IN)=DR(IN) +b*p(IN);
Dp=PDR +b* Dp;
end
% Step size
% Dp=P(p); % =P(DR(IN)) +b P(p(IN));
pDDp=Dp‘*Dp;
a=Residual‘*Dp/(pDDp);
% Update coefficients
s=s+a*p;
% New Residual and inner products
Residual=Residual-a*Dp;
DR=Pt(Residual);
% select new element
[v I]=max(abs(DR(NI))./normA(NI)‘);
INI = NI(I);
if weakness ~= 1
[vals inds] = sort(abs(DR),‘descend‘);
I=inds( find( vals >= alpha * v ) );
end
IN = union(IN,INI);
NI(I) = [];
if strcmp(STOPCRIT,‘M‘) & length(IN) >= STOPTOL
IN=IN(1:STOPTOL);
end
else
error(‘Undefined option for GradSteps, use ‘‘auto‘‘ or an integer.‘)
end
elseif isa(GradSteps,‘numeric‘)
% Do GradSteps gradient steps
count=1;
while count<=GradSteps
% Update direction
if iter==1
p(IN)=DR(IN);
Dp=P(p);
else
MASK(IN)=1;
PDR=P(DR.*MASK);
b=-Dp‘*PDR/pDDp;
p(IN)=DR(IN) +b*p(IN);
Dp=PDR +b* Dp;
end
% Step size
% Dp=P(p);
pDDp=Dp‘*Dp;
a=Residual‘*Dp/(pDDp);
% Update coefficients
s=s+a*p;
% New Residual and inner products
Residual=Residual-a*Dp;
DR=Pt(Residual);
count=count+1;
end
% select new element
[v I]=max(abs(DR(NI))./normA(NI)‘);
INI = NI(I);
if weakness ~= 1
[vals inds] = sort(abs(DR),‘descend‘);
I=inds( find( vals >= alpha * v ) );
end
IN = union(IN,INI);
NI(I) = [];
if strcmp(STOPCRIT,‘M‘) & length(IN) >= STOPTOL
IN=IN(1:STOPTOL);
end
else
error(‘Undefined option for GradSteps, use ‘‘auto‘‘ or an integer.‘)
end
ERR=Residual‘*Residual/n;
if comp_err
err_mse(iter)=ERR;
end
if comp_time
iter_time(iter)=toc;
end
if strcmp(STOPCRIT,‘M‘)
if length(IN) >= STOPTOL
done =1;
elseif verbose && toc-t>10
display(sprintf(‘Iteration %i. --- %i selected elements‘,iter ,length(IN)))
t=toc;
end
elseif strcmp(STOPCRIT,‘mse‘)
if comp_err
if err_mse(iter)<STOPTOL;
done = 1;
elseif verbose && toc-t>10
display(sprintf(‘Iteration %i. --- %i mse‘,iter ,err_mse(iter)))
t=toc;
end
else
if ERR<STOPTOL;
done = 1;
elseif verbose && toc-t>10
display(sprintf(‘Iteration %i. --- %i mse‘,iter ,ERR))
t=toc;
end
end
elseif strcmp(STOPCRIT,‘mse_change‘) && iter >=2
if comp_err && iter >=2
if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL);
done = 1;
elseif verbose && toc-t>10
display(sprintf(‘Iteration %i. --- %i mse change‘,iter ,(err_mse(iter-1)-err_mse(iter))/sigsize ))
t=toc;
end
else
if ((oldERR - ERR)/sigsize < STOPTOL);
done = 1;
elseif verbose && toc-t>10
display(sprintf(‘Iteration %i. --- %i mse change‘,iter ,(oldERR - ERR)/sigsize))
t=toc;
end
end
elseif strcmp(STOPCRIT,‘corr‘)
if max(abs(DR)) < STOPTOL;
done = 1;
elseif verbose && toc-t>10
display(sprintf(‘Iteration %i. --- %i corr‘,iter ,max(abs(DR))))
t=toc;
end
end
% Also stop if residual gets too small or maxIter reached
if comp_err
if err_mse(iter)<1e-16
display(‘Stopping. Exact signal representation found!‘)
done=1;
end
else
if iter>1
if ERR<1e-16
display(‘Stopping. Exact signal representation found!‘)
done=1;
end
end
end
if iter >= MAXITER
display(‘Stopping. Maximum number of iterations reached!‘)
done = 1;
end
if ~done
iter=iter+1;
oldERR=ERR;
end
end
if nargout >=2
err_mse = err_mse(1:iter);
end
if nargout ==3
iter_time = iter_time(1:iter);
end
if verbose
display(‘Done‘)
end
三、运行结果
四、备注
版本：2014a
完整代码或代写加1564658423