心电信号基于matlab NLM时间序列心电信号去噪含Matlab源码 1547期

Posted 2021-12-24 紫极神光

一、心电信号简介

0 引言
心电信号是人类最早研究的生物信号之一, 相比其他生物信号更易于检测, 且具有直观的规律。心电图的准确分析对心脏病的及早治疗有重大的意义。人体是一个复杂精密的系统, 有许多不可抗的外界因素, 得到纯净的心电信号非常困难。可以采用神经网络算法去除心电信号的噪声, 但这种方法存在训练难度大、耗时长的缺点。小波变换在处理非线性、非平稳且奇异点较多的信号时具有一定的优越性, 近年来许多学者使用其对心电信号进行研究。

1 心电信号简介
心电信号由以下几个波段组成, 一个典型的心电图如图1所示。

图1 典型心电图
(1) P波:反映心房肌在除极过程中的电位变化过程;
(2) P-R间期:反映的是激动从窦房结通过房室交界区到心室肌开始除极的时限;
(3) QRS波群:反映心室肌除极过程的电位变化;
(4) T波:代表心室肌复极过程中所引起的电位变化;
(5) S-T段:从QRS波群终点到达T波起点间的一段水平线[2];
(6) Q-T间期:心室从除极到复极的时间[3];
(7) U波:代表动作电位的后电位。
由于心电信号十分微弱, 且低频, 极易受到干扰, 不同的干扰源的噪声虽是随机的, 但来自同一个干扰源的噪声往往具有同一类特征。分析干扰的来源, 针对不同的来源使用合适的处理方法, 是数据采集重点考虑的一个问题。常见干扰有3种: (1) 工频干扰; (2) 基线漂移; (3) 肌电干扰。其中已经证明小波变换在抑制心电信号的工频干扰方面具有较大优势。具体噪声频带如表1所示。
表1 心电信号以及主要噪声频带

二、部分源代码

% EXAMPLES of calling NLM code


% load the data
ecg100=load('ecg100.txt');
ix=1:length(ecg100);  % vector of time samples, for plotting

% set parameters - see Tracey & Miller for a discussion
PatchHW=10;  % patch half-width; ~ size of smallest feature, in samples
P = 1000;    % neighborhood search width; wider=more computation but more
             % channces to find a similar patch

% set bandwidth for NLM - here set by an "eyeball estimate" 
lambda = 0.6*.02;

% now denoise 
[dnEcg,debugEcg]= NLM_1dDarbon(ecg100,lambda,P,PatchHW);


% now, create a signal with 10 dB SNR
[noisySig,targetNoiseSigma] = createSignalPlusNoise(ecg100,10);

% since we know sigma here, use that to set NLM bandwidth parameter
lambda=0.6*targetNoiseSigma;  

% and denoise....
dnEcg2= NLM_1dDarbon(noisySig,lambda,P,PatchHW);

%% plot results

xlim_vals = [1000 2000];

figure,
subplot(221),
plot(ix,ecg100)
xlim(xlim_vals)
title('Original ECG100 data')

subplot(223),
plot(ix,dnEcg,'r')
xlim(xlim_vals)
xlabel('Time, samples')
title('Denoised ECG100')

subplot(222),
plot(ix,noisySig)
xlim(xlim_vals)
title('ECG100 +  noise')


subplot(224),
plot(ix,dnEcg2,'r')
xlim(xlim_vals)
xlabel('Time, samples')
title('Denoised ECG100 + noise')
function [denoisedSig,debug] = NLM_1dDarbon(signal,lambda,P,PatchHW)
% function [denoisedSig,debug] = NLM_1dDarbon(signal,lambda,P,PatchHW)
% Implements fast NLM method of Darbon et al, for a 1-D signal
% INPUTS:
% signal: input signal (vector)
% lambda: Gaussian scale factor
% P: max search distance
% PatchWH: patch half-width
% OUTPUTS:
% denoisedSig: the NLM-denoised signal
% debug: structure containing various quantitities that can help debug
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =========================================================================
% PAPER INFO:
%       Brian Tracey and Eric Miller, "Nonlocal means denoising of ECG signals",  
%       IEEE Transactions on Biomedical Engineering, Vol 59, No 9, Sept
%       2012, pages 2383-2386
% -------------------------------------------------------------------------
%       PLEASE CITE THIS PAPER, IF YOU USE THIS CODE FOR ACADEMIC PURPOSES
% -------------------------------------------------------------------------
%       For all inquiries, please contact author Brian Tracey(btracey[at]alum.mit.edu)
%
%       Last Update 05/09/2013
% =========================================================================
%

if length(P)==1,  % scalar has been entered; expand into patch sample index vector
    Pvec = -P:P;
else
   Pvec = P;  % use the vector that has been input  
end
debug=[];
N = length(signal);

denoisedSig = NaN*ones(size(signal));

% to simpify, don't bother denoising edges
iStart=1+PatchHW+1;
iEnd = N-PatchHW;
denoisedSig(iStart:iEnd) = 0;

debug.iStart = iStart;
debug.iEnd = iEnd;

% initialize weight normalization

% convert lambda value to 'h', denominator, as in original Buades papers
Npatch = 2*PatchHW+1;
h = 2*Npatch*lambda^2;

for idx = Pvec  % loop over all possible differences: s-t
 
    % do summation over p  - Eq. 3 in Darbon
    k=1:N;
    kplus = k+idx;
    igood = find(kplus>0 & kplus<=N);  % ignore OOB data; we could also handle it
    SSD=zeros(size(k));
    SSD(igood) = (signal(k(igood))-signal(kplus(igood))).^2;
    Sdx = cumsum(SSD);
   

    for ii=iStart:iEnd  % loop over all points 's'
        distance = Sdx(ii+PatchHW) - Sdx(ii-PatchHW-1); % Eq 4; this is in place of point-by-point MSE
        % but note the -1; we want to icnlude the point ii-iPatchHW

      
        end

    end
end % loop over shifts

% now apply normalization
denoisedSig = denoisedSig./(Z+eps);
debug.Z = Z;

return

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]焦运良，邢计元，靳尧凯.基于小波变换的心电信号阈值去噪算法研究[J].信息技术与网络安全. 2019,38(05)