function Xd=FecgDetrFilt(X,fs,cName,graph,dbFlag)
% -------------------------------------------------------------------------------------------
% ECG detrending
%
% A baseline signal is estimated applying a low pass Butterworth filter
% in forward and backward direction. The filter cut frequency is 3.17Hz.
% Each detrended signal is obtained as difference between the original signal
% and the estimated baseline.
% In case of residual artifacts due to fast baseline movements, median filtering is used.
%
% function Xd=FecgDetrFilt(X,fs,cName,graph,dbFlag)
%
% X : signal matrix (one signal per column)
% fs : sampling frequency
% cName : record name
% graf : flag enabling figure drawing
% dbFlag : flag enabling figure drawing for tuning and debugging
%
% -------------------------------------------------------------------------------------------
% Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% -------------------------------------------------------------------------------------------
% This program is free software; you can redistribute it and/or modify it under the terms
% of the GNU General Public License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.
%
% This program is distributed "as is" and "as available" in the hope that it will be useful,
% but WITHOUT ANY WARRANTY of any kind; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
% -------------------------------------------------------------------------------------------
if(nargin<3), cName=''; end
if(nargin<4), graph=0; end
if(nargin<5), dbFlag=0; end
graphD= graph &dbFlag ;
graphSpt= graph &dbFlag ;
fprintf('\n --------------------------------------------------------- \n');
[progpath, progname] = fileparts(which(mfilename));
fprintf('Program: %s, record name: %s\n', progname, cName);
%-------------------------------------------------------------
% recording duration
[ndt, ns]=size(X);
vtime= [1:ndt]/fs;
%-------------------------------------------------------------------------------
fmaxd=5; % --> 3.17 Hz = filtfilt 3dB cut frequency
fmaxn = fmaxd/(fs/2);
[b,a]= butter(1,fmaxn,'low');
Xb=filtfilt(b,a,X); % estimated baseline
% remove baseline from original ecg
Xd=X-Xb;
%-------------------------------------------------------------------------------
RRpmean=1.;
ww=fix(RRpmean*1.5*fs); % wide windows (many contain at least one QRS)
for is=1:ns
x=Xd(:,is); n=length(x);
[xmeaMi,xmeaMa] = meanMiMaSc(x,ww,8,8);
thmima=1.3;
xthmi=xmeaMi*thmima; xthma=xmeaMa*thmima;
if(graphD)
figure, set(gcf,'Color','white'); hold on;
plot(x,'b.-');
plot(xmeaMi*ones(n,1),'m');
plot(xmeaMa*ones(n,1),'m');
plot(xthmi*ones(n,1),'r');
plot(xthma*ones(n,1),'r');
title(['EcgBaseline_atest:',' x(b), mi&ma (m), th (r)']);
end
iibadDetr=find(x<xthmi | xthma<x);
% iibadDetr(find(diff(iibadDetr)>1*fs));
if(length(iibadDetr)>1.3*fs)
fprintf('Switching to median filtering,rec=%s, is=%d, bads=%d\n',cName, is, length(iibadDetr));
wm=fix(0.26*fs); wm=2*floor(wm/2)+1; % window width
Xb(:,is) = medfilt1mit(X(:,is),wm,1);
x=X(:,is)-Xb(:,is); % remove baseline from original ECG
[xmeaMi,xmeaMa] = meanMiMaSc(x,ww,8,8);
thmima=1.3;
xthmi=xmeaMi*thmima; xthma=xmeaMa*thmima;
iinf=find(x<xthmi); x(iinf)=xthmi;
isup=find(x>xthma); x(isup)=xthma;
fmaxd=150; fmaxn=fmaxd/(fs/2); % Normalized cut-off frequency
[b,a]= butter(1,fmaxn,'low');
x=filtfilt(b,a,x); % low-pass filtered signal
Xd(:,is)=x;
if(graphSpt)
figure; freqz(conv(b,b(end:-1:1)),conv(a,a(end:-1:1)),[0:0.005:2],fs);
end
end
end
if(graphSpt)
for is=1:ns,
figure, set(gcf,'Color','white');
pwelch(X(:,is),[],[],[],fs);
title('Original ECG Welch spectrum');
shg
end
end
if(graphSpt)
for is=1:ns,
figure, set(gcf,'Color','white');
pwelch(Xd(:,is),[],[],[],fs);
title('Detrended ECG Welch spectrum');
shg
end
end
if(graphD)
figure, set(gcf,'Color','white');
for is=1:ns,
subplot(ns,1,is), hold on, plot(vtime,X(:,is),'r'), plot(vtime,Xb(:,is),'b');
wgmimaV=mimaxscG(X(:,is),0,0,.1);
ylim(wgmimaV);
if(is~=ns), set(gca,'XTickLabel',''); end
if(is==1), title('Original ECG (r) & estimated baseline(b)'); end
end
shg
end
if(graph)
figure, set(gcf,'Color','white');
for is=1:ns,
subplot(ns,1,is), hold on, plot(vtime,X(:,is),'r'), plot(vtime,Xd(:,is),'b');
wgmimaV=mimaxscG(Xd(:,is),0,0,.4);
ylim(wgmimaV);
if(is~=ns), set(gca,'XTickLabel',''); end
if(is==1), title([cName,': original(r) & detrended(b) ECG']); end
end
end
%if(dbFlag & learning) PlotSgnMrkNc(X', QRSa, fs, cName); end
end %== function ================================================================
%