function NeedFiltering = check_notch(ecg,fc,fs,debug)
% this function checks if there is an energy peak in the frequency band
% +/-1 Hz around fc. The reason for this function is to avoid doing too much
% filtering if we can, i.e. decide to Notch filter or not the signal, 
% because this might create ripples that affect the ecg.
% 
% inputs
%   ecg:    a single channel ecg
%   fc:     the central frequency we are studying
%   fs:     sampling frequency
% 
% ouputs
%   NeedFiltering: boolean. 1 if there is an energy peak and 0 otherwise
%
%
% FECG extraction toolbox, version 1.0, Sept 2013
% Released under the GNU General Public License
%
% Copyright (C) 2013  Joachim Behar
% Oxford university, Intelligent Patient Monitoring Group - Oxford 2013
% joachim.behar@eng.ox.ac.uk
%
% Last updated : 02-08-2013
%
% 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 in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

% == managing inputs
if nargin<2; error('check_notch: wrong number of input arguments \n'); end;
if nargin==2; fs=1000; debug=0; end;
if nargin==3; debug=0; end;
if size(ecg,2)>size(ecg,1); ecg=ecg'; end;

% == general parameters
NFFT = 10000;        % number of NFFT coefficients
pace = fs/NFFT;      % one data point corresponds to pace Hz
df = 4;              % delta frequency: we look at +/-df around fc
NeedFiltering = 0;

% == core function
try
    % == define spectral band we are looking for a peak into
    specBand = round(fc/pace)-round(df/pace):round(fc/pace)+round(df/pace);

    % == PSD
    [Pxx F1] = pwelch(ecg(:,1),5*fs,fs,NFFT,fs,'onesided');
    PxxdB = 10*log10(Pxx/(mean(Pxx)));

    % == look for local max in specBand
    [~,pos] = max((PxxdB(specBand)));

    % == and the frequency it corresponds to
    CorrespondingFreq = F1(round(fc/pace)-round(df/pace)+pos-1);

    % == is it in +/-1Hz around fc?
    if CorrespondingFreq<fc+1 && CorrespondingFreq>fc-1; NeedFiltering = 1; end;
catch ME
    for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
    NeedFiltering = 1;
end

% == debug
if debug
    plot(F1,PxxdB); hold on, plot(F1(specBand),1,'+r');
    hold on, plot(CorrespondingFreq,PxxdB(round(fc/pace)-round(df/pace)+pos-1),'+g');
    set(findall(gcf,'type','text'),'fontSize',14,'fontWeight','bold');
end

end