一、获取代码方式
获取代码方式1：
完整代码已上传我的资源：【心电信号】基于matlab心电信号去除基线漂移【含Matlab源码 955期】
获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。
备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；
二、心电信号中基线漂移的去除方法比较及算法实现简介
0 引言
心电信号是生物医学信号中一种低幅低频的微弱信号。基线漂移一般是由于人体呼吸、电极运动等引起, 其一般小于1 Hz。常用的去除方法有中值滤波法、小波变换法、形态学滤波法等。本文简单介绍了几种常用方法并用Matlab实现。以下用于处理的原始心电信号是截取自MIT-BIH数据库的序号为103和114两种数据中的一段数据, 我们通过对比基线漂移严重和微弱两种情况, 比较各种方法的优劣之处, 其中103为较严重漂移的心电数据, 114为漂移较为微弱的心电数据。
1 方法
1.1 中值滤波法
中值滤波的定义就是取信号中的某一点前的N (N为奇数) 个数据, 将这N个数据进行排序, 取数据的中值, 将数据的中值作为该点的数据。将整段信号都做以上处理, 提取出漂移信号, 将原始信号减去漂移信号, 便得到滤波后的效果。测试效果如图1、图2所示。
图1 中值滤波去除严重基线漂移
1.2 高通滤波法
基线漂移的频率成分较低, 所以用高通的方法对基线漂移的去除也有一定的效果。IIR滤波器能用较低的阶数获得较高的频率选择性, 但不具有线性相位。由于心电信号的处理对线性相位的要求不是很高, 所以可以选择IIR滤波器。这里我们选用butterworth滤波器, 截止频率为1 Hz。测试效果如图3、图4所示。
图2 中值滤波去除微弱基线漂移
图3 高通滤波去除严重基线漂移
图4 高通滤波去除微弱基线漂移
1.3 小波变换法
根据小波变换多分辨分析理论, 对含噪信号进行多尺度小波分解后, 可以利用信号和噪声在频谱和能量分布上的不同, 直接将噪声所对应的小波分解尺度上的细节分量去除, 然后利用小波逆变换进行信号重构, 就可以有效去除信号中的噪声成份。经过试验得到, 采用小波函数db8对原始信号进行8个尺度的分解, 第8层分解的信号频率与基线漂移的频率接近, 因此将前7层信号重构, 去除第8层信号, 即可得到去除基线漂移的信号。测试效果如图5、图6所示。
图5 小波变换去除严重基线漂移
图6 小波变换去除微弱基线漂移
三、部分源代码
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 == 1n=n1;
elseif n1 == 1x=x';n=n2;
elseerror('x must be a vector.');
endsigsize = 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 verbosedisplay('Initialising...')
endswitch nargout case 3comp_err=true;comp_time=true;case 2 comp_err=true;comp_time=false;case 1comp_err=false;comp_time=false;case 0error('Please assign output variable.') otherwiseerror('Too many output arguments specified')
end% Put option into nice format
Options={};
OS=nargin-3;
c=1;
for i=1:OSif isa(varargin{i},'cell')CellSize=length(varargin{i});ThisCell=varargin{i};for j=1:CellSizeOptions{c}=ThisCell{j};c=c+1;endelseOptions{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:OSswitch 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.'); endcase {'P_trans'} if isa(Options{i+1},'function_handle'); Pt = Options{i+1}; else error('P_trans must be function _handle. Exiting.'); endcase {'maxIter'}if isa(Options{i+1},'numeric'); MAXITER = Options{i+1}; else error('maxIter must be a number. Exiting.'); endcase {'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.'); endcase {'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.'); endcase {'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.'); endotherwiseerror('Unrecognised option. Exiting.') end
endif strcmp(STOPCRIT,'M') maxM=STOPTOL;
elsemaxM=MAXITER;
endif nargout >=2err_mse = zeros(maxM,1);
end
if nargout ==3iter_time = zeros(maxM,1);
endif 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') tryif isa(Pt,'function_handle'); P=A;else error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); endcatch 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.'); endif initial_given ==1;IN = find(s_initial);Residual = x-P(s_initial);s = s_initial;oldERR = Residual'*Residual/n;
elseIN = [];Residual = x;s = s_initial;sigsize = x'*x/n;oldERR = sigsize;
endmask=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.')endif verbosedisplay('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 ) );
endIN = union(IN,I);
if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOLIN=IN(1:STOPTOL);
end
MASK=zeros(size(DR));
pDDp=1;
done = 0;
iter=1;while ~done% Select new elementif isa(GradSteps,'char')if strcmp(GradSteps,'auto')% finished=0;
% while ~finished% Update direction if iter==1p(IN)=DR(IN);Dp=P(p);elseMASK(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 productsResidual=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 ) );endIN = union(IN,INI);NI(I) = [];if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOLIN=IN(1:STOPTOL);endelseerror('Undefined option for GradSteps, use ''auto'' or an integer.')endelseif isa(GradSteps,'numeric') % Do GradSteps gradient stepscount=1;while count<=GradSteps% Update direction if iter==1p(IN)=DR(IN);Dp=P(p);elseMASK(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 productsResidual=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 ) );endIN = union(IN,INI);NI(I) = [];if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOLIN=IN(1:STOPTOL);endelseerror('Undefined option for GradSteps, use ''auto'' or an integer.')endERR=Residual'*Residual/n;if comp_errerr_mse(iter)=ERR;endif comp_timeiter_time(iter)=toc;endif strcmp(STOPCRIT,'M')if length(IN) >= STOPTOLdone =1;elseif verbose && toc-t>10display(sprintf('Iteration %i. --- %i selected elements',iter ,length(IN))) t=toc;endelseif strcmp(STOPCRIT,'mse')if comp_errif err_mse(iter)<STOPTOL;done = 1; elseif verbose && toc-t>10display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter))) t=toc;endelseif ERR<STOPTOL;done = 1; elseif verbose && toc-t>10display(sprintf('Iteration %i. --- %i mse',iter ,ERR)) t=toc;endendelseif strcmp(STOPCRIT,'mse_change') && iter >=2if comp_err && iter >=2if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL);done = 1; elseif verbose && toc-t>10display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) t=toc;endelseif ((oldERR - ERR)/sigsize < STOPTOL);done = 1; elseif verbose && toc-t>10display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize)) t=toc;endendelseif strcmp(STOPCRIT,'corr') if max(abs(DR)) < STOPTOL;done = 1; elseif verbose && toc-t>10display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR)))) t=toc;endend% Also stop if residual gets too small or maxIter reachedif comp_errif err_mse(iter)<1e-16display('Stopping. Exact signal representation found!')done=1;endelseif iter>1if ERR<1e-16display('Stopping. Exact signal representation found!')done=1;endendendif iter >= MAXITERdisplay('Stopping. Maximum number of iterations reached!')done = 1; endif ~doneiter=iter+1;oldERR=ERR;end
endif nargout >=2err_mse = err_mse(1:iter);
end
if nargout ==3iter_time = iter_time(1:iter);
end
if verbosedisplay('Done')
end四、运行结果
五、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]陈苹,叶继伦,张旭,陈刚.心电信号中基线漂移的去除方法比较及算法实现[J].中国医疗器械杂志. 2018,42(05)
支付宝扫一扫
微信扫一扫
.png)
.png)
.png)
.png)
