function Xf=FecgNotchFilt(X,fs,cName,graph,dbFlag)
% -----------------------------------------------------------------------------------
%   ECG power line canceling by notch filtering
%
% Xf=FecgNotchFilt(X,fs,cName,graph,dbFlag)
%
% X          : input signal matrix (one signal per column)
% fs         : sampling frequency
% cName      : record name
% graph      : flag enabling figure drawing
% Xf         : 
%
% Author: Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% -----------------------------------------------------------------------------------

if(nargin<3), cName=''; end
if(nargin<4), graph=0; end
if(nargin<5), dbFlag=0; end
graphD= graph &dbFlag ;
graphSpt= graph &dbFlag ;
%-------------------------------------------------------------
% recording duration
[ndt, ns]=size(X);
vtime= [1:ndt]/fs;
Xf=zeros(size(X));
for is=1:ns
    if(dbFlag && graphSpt), figure; pwelch(X(:,is),[],[],[],fs);
        title('Detrended ECG Welch spectrum');
    end
    [Px,Fv]=pwelch(X(:,is),[],[],[],fs);
    iF50w=find(49<Fv & Fv<51);  F50w=Fv(iF50w);
    [maxP50,imaxP50]=max(Px(iF50w));
    f50mP=F50w(imaxP50);
    iFes50w=find((46<Fv & Fv<49) | (51<Fv & Fv<54));
    pem50w=mean(Px(iFes50w)); pestd50w=std(Px(iFes50w));
    d50P=maxP50-pem50w;
    e50=d50P>5*pestd50w;
    
    iF60w=find(59<Fv & Fv<61);  F60w=Fv(iF60w);
    [maxP60,imaxP60]=max(Px(iF60w));
    f60mP=F60w(imaxP60);
    iFes60w=find((56<Fv & Fv<59) | (61<Fv & Fv<64));
    pem60w=mean(Px(iFes60w)); pestd60w=std(Px(iFes60w));
    d60P=maxP60-pem60w;
    e60=d60P>5*pestd60w;
    
    if(e50 || e60)
        if(d50P>d60P), fnotch=f50mP;
        else fnotch=f60mP;  end
        
        % Notch filter
        fnotchn=fnotch/fs; bwn=1/fs; grafSptn=0;
        
        Xf(:,is)= filtNotchFB(X(:,is),fnotchn,bwn, grafSptn);
        Xf(:,is)= filtNotchFB(Xf(:,is),2*fnotchn,bwn, grafSptn);
        Xf(:,is)= filtNotchFB(Xf(:,is),3*fnotchn,bwn, grafSptn);
        Xf(:,is)= filtNotchFB(Xf(:,is),4*fnotchn,bwn, grafSptn);
    else
        Xf(:,is)=X(:,is);
    end
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,Xf(:,is),'b');
        wgmimaV=mimaxscG(X(:,is),0,0,.1);
        ylim(wgmimaV);
        if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': original(r) & notch filtered(b) ECG']); end
    end
    shg
end

return
end %== function ================================================================
%