function xc=ImpArtElimS(x,thE,wm,pvsc,plotFlag)
% -------------------------------------------------------------------------------------------
%   Impulsive Artifact Canceling from a signal.
%     xc=ImpArtElimS(x,thE,wm,pvsc,plotFlag)
%   Input parameters:
%     x        : array of data to be cleaned from wild points
%     thE      : threshold on absolute error
%     wm       : window length for "deviation" estimation
%     pvsc     : % of value to discard in deviation estimation
%     plotFlag :
%
%   Ouput parameters:
%     xc       : corrected array of data
%
% -------------------------------------------------------------------------------------------
% 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
% -------------------------------------------------------------------------------------------
% 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 < 1), [x, xv, thE, wm,pvsc,plotFlag]=ImpArtElim_atest; testFlag=1; 
else
    if (nargin < 2 || isempty(thE)), thE=4;  end
    if (nargin < 3 || isempty(wm)),  wm=min(33, length(x)-1);  end  % default values for RR interval series
    if (nargin < 4 ),  pvsc=10;  end
    if (nargin < 5 ),  plotFlag=0;  end
    testFlag=0; 
end
[nr,nc]=size(x);
if(nc>nr), n=nc; x=x'; else n=nr; end

debugFlag=0;

if(wm>n), fprintf('Error: parameter wm must be smaller than length(x)\n'); xc=[]; return; end
% --------
xc=x;

if(debugFlag)
    wLd=60*5000;
    xin=x(1:min(wLd,length(x)));
    xmed= median(xin);    % median value
    xad=abs(xin-xmed);    % vector of absolute error respect to median
    admed=median(xad);    % median of absolute error (deviation)
    admedsd=1.483 *admed;  % robust estimate of standard deviation
    % for a Gaussian signal
    figure; plot(xin);   axis tight;
    figure; plot(xad);   axis tight;
    hold on; plot([1; n],admedsd*[thE; thE],'r');
end

wm=2*floor(wm/2)+1;

xmed = medfilt1mit(xc,wm,1);
xad=abs(xc-xmed);
xadm=maxsc(xad(find(xad>0)),pvsc);

if(plotFlag)
    figure;  hold on;
    plot(x,'r.-'); plot(xmed,'b.-');
    ht=title(['ImpArtElim_atest:',' x (r) - xmed (b)']);
    set(ht,'Interpreter','none');
    shg;   % it need in Linux
    figResize(0, 1, 1, .35); axis tight;
    figure;  hold on;
    plot(xad,'b.-');
    plot(xadm*ones(n,1),'m');
    plot(thE*xadm*ones(n,1),'r');
    ht=title(['ImpArtElim_atest:',' xad (b), xadm (m), th (r)']);
    set(ht,'Interpreter','none');
    shg;   % it need in Linux
    figResize(0, 1, 1, .35); axis tight;
end

kv=find(xad - thE*xadm > 0);
if(~isempty(kv))
    i=1;
    while i<=length(kv)
        iivx=kv(i);
        while (i<length(kv) && kv(i+1)==kv(i)+1), i=i+1; end
        ifvx=kv(i);
        xck=(xc(max(iivx-1,1))+xc(min(ifvx+1,n)))/2;
        xc(iivx:ifvx)= xck;
        i=i+1;
    end
end

if(plotFlag)
    figure;  hold on;
    if(testFlag), plot(xv,'b.-');  end
    plot(x,'r.-'); plot(xc,'g.-');
    ht=title(['ImpArtElim_atest:',' xv (b) - x (r) - xc (g)']);
    set(ht,'Interpreter','none');
    shg;   % it need in Linux
    figResize(0, 1, 1, .35); axis tight;
end

if(nc>nr), xc=xc'; end
if(testFlag),  xc=[];  end
%
end %== function ================================================================

% -------------------------------------------------------------------------------------------------
function [x, xv, thE, wm, pvsc, plotFlag]=ImpArtElim_atest
close all
plotFlag=1;
thE=1.4;  % thE=2;
wm=33;  % wm=33;
pvsc=10;
% --   artificial RR series
% fc_bpm=75bpm => fc=fc_bpm/60= 1.25Hz  => Tc= 0.8s
% fHF=0.25Hz => fHF=fHF/fc
fc=1.25; Tc=1/fc;
fHF=0.25; aHF=0.03*Tc;
fLF=0.1;   aLF= 0.02*Tc;
fHFn=fHF/fc; fLFn=fLF/fc;  % normalized frequencies (cycles/point)
t = (0:Tc:3*60)';
xv = Tc+ aLF*sin(2*pi*fLF*t)+ aHF*sin(2*pi*fHF*t)+ 0.005*Tc*randn(size(t));
fprintf('--- Test for "ImpArtElim" routine ---\n');
fprintf('Artificial  RR interval series - ');
fprintf(' normalized frequencies: LF=%3.2g, HF=%3.2g\n', fLFn, fHFn);
fprintf('n=%d,  thE=%3.2g, wm=%d\n', length(xv), thE, wm);
% --------------------------
% add artifact to the x series
x=xv;
xmed=median(x); mdaerr=median(abs(x-xmed));
fprintf('median= %5.1f,  medAbsErr= %5.2f\n',xmed, mdaerr);
ia=floor(length(x)/4); x(ia)= x(ia)+14*mdaerr; x(ia+2)= x(ia+2)+14*mdaerr;
ia=floor(length(x)/3); x(ia:ia+5)= x(ia:ia+5)+14*mdaerr;
ia=floor(length(x)/2); x(ia:ia+3)= x(ia:ia+3)+8*mdaerr;
ia=floor(length(x)*2/3); x(ia:ia+2)= x(ia:ia+2)+4*mdaerr;
ia=floor(length(x)*3/4); x(ia)= x(ia)-8*mdaerr;
return
end %== function ================================================================
%