function [ecgf] = preprocessing(ecg,paramStruct,varargin)
% this function preprocess the ecg prior applying any source separation 
% technique.
%
% inputs  
%   ecg: 4x60000 (4 channels and 1min of signal at 1000Hz) matrix of
%   abdominal ecg channels.
%   paramStruct: the following parameters needs to be specified
%       paramStruct.fs
%       paramStruct.NotchCheck
%       paramStruct.fhigh
%       paramStruct.fbas
%       paramStruct.fhighNbCoeff
%       paramStruct.fbasNbCoeff
%   varargin
%       createFilters: recreate the filter coefficients?
%       path2filtersCoeff: path to where to save the filter 
%                          coefficients or load them from.
%       debug: enter debug mode
%
% outputs
%   ecgf: normalized and filtered ecg. Filtering
%                     includes Notch, high frequency, baseline wander removal
%                     and notmalisation in range [-1 1]. A non linear step is
%                     applied with the tanh so this might note be suited for
%                     morphological analysis but is usefull for FQRS detection.
%
% IMPORTANT NOTE: too many coefficients  for the IIR filters resulted in 
% ringing. So need to use as few as possible and avoid any useless filtering 
% step which is why I am looking whether there is a peak at 50/60Hz i.e. 
% whether or not to perform a Notch filtering step.
%
%
% 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.

createFilters = 0;
path2filtersCoeff = 'subfunctions/FiltersCoeff/';
debug = 0;
% == manage inputs
switch nargin
    case 2
    case 3
        createFilters=varargin{1};       
    case 4
        createFilters=varargin{1};
        path2filtersCoeff=varargin{2};
    case 5
        createFilters=varargin{1};
        path2filtersCoeff=varargin{2};
        debug=varargin{3};
    otherwise
        error('preprocessing: wrong number of input arguments \n');              
end
if size(ecg,2)>size(ecg,1); ecg = ecg'; end;
if ~exist(path2filtersCoeff,'dir'); mkdir(path2filtersCoeff); end;
currFold = pwd;

% == core function
try 
    % some AECG (the artificial ones) saturate sometimes. Because they were
    % rescaled (after being saved in WFDB format) then this leads to weird
    % disptoportionated saturation. The lines below are aimed at solving
    % this problem. (If I open the records in physical units and run the
    % filters, I am gettnig some NaN..)
    % - manage saturation for wfdb records
    % NOTE: this is wrong but works by luck somehow- will need to be fixed
    % at some point!
    [ind_t,ind_ch] = find(ecg>32750);
    ecg(ind_t,ind_ch) = median(ecg(ind_t-5:ind_t-5,ind_ch));
    [ind_t,ind_ch] = find(ecg<-32750);
    ecg(ind_t,ind_ch) = median(ecg(ind_t-5:ind_t-5,ind_ch)); % should be -32767 but they have rescaled the artificial signals

    % - manage NaN for .cvs files
    ind_nan = isnan(ecg);
    for ind_ch_nan=1:4
        if sum(ind_nan(:,ind_ch_nan))>0
            nan_ch_pos = find(ind_nan(:,ind_ch_nan)); % position of the NaN
            if ~isempty(nan_ch_pos)
                for jj=1:length(nan_ch_pos)
                    
                    start = nan_ch_pos(jj)-5;
                    stop = nan_ch_pos(jj)+5;
                    if start<0; start=0; end;
                    if stop>length(ecg); stop=length(ecg); end;
                    
                    vecTrans = ecg(start:stop,ind_ch_nan);
                    vecTrans(isnan(vecTrans)) = [];
                    if ~isempty(vecTrans)
                        ecg(nan_ch_pos(jj),ind_ch_nan) = median(vecTrans);
                    end
                end 
            end
        end
    end
    
    % == Filtering
    if createFilters
        % == 60Hz Notch
        wo = 60/(paramStruct.fs/2);  bw = wo/35;
        [bn_60Hz,an_60Hz] = iirnotch(wo,bw);
        dlmwrite(strcat(path2filtersCoeff,'bn_60Hz.txt'),bn_60Hz,'precision',20);
        dlmwrite(strcat(path2filtersCoeff,'an_60Hz.txt'),an_60Hz,'precision',20);    
        % == 50Hz Notch
        wo = 50/(paramStruct.fs/2);  bw = wo/35;
        [bn_50Hz,an_50Hz] = iirnotch(wo,bw);
        dlmwrite(strcat(path2filtersCoeff,'bn_50Hz.txt'),bn_50Hz,'precision',20);
        dlmwrite(strcat(path2filtersCoeff,'an_50Hz.txt'),an_50Hz,'precision',20);
        % == high freq
        [b_lp,a_lp] = butter(paramStruct.fhighNbCoeff,paramStruct.fhigh/(paramStruct.fs/2),'high');
        dlmwrite(strcat(path2filtersCoeff,'b_lp.txt'),b_lp,'precision',20);
        dlmwrite(strcat(path2filtersCoeff,'a_lp.txt'),a_lp,'precision',20);
        % == baseline
        [b_bas,a_bas] = butter(paramStruct.fbasNbCoeff,paramStruct.fbas/(paramStruct.fs/2),'high');
        dlmwrite(strcat(path2filtersCoeff,'b_bas.txt'),b_bas,'precision',20);
        dlmwrite(strcat(path2filtersCoeff,'a_bas.txt'),a_bas,'precision',20);
        %fvtool(b,a);     
    else
        cd(path2filtersCoeff);
            an_60Hz = importdata('an_60Hz.txt');
            bn_60Hz = importdata('bn_60Hz.txt');
            an_50Hz = importdata('an_50Hz.txt');
            bn_50Hz = importdata('bn_50Hz.txt');
            b_lp = importdata('b_lp.txt');
            a_lp = importdata('a_lp.txt');
            b_bas = importdata('b_bas.txt');
            a_bas = importdata('a_bas.txt');
        cd(currFold)
    end

    % == Notch fileting
    if paramStruct.NotchCheck
        % == remove 60Hz?
        NeedFiltering = check_notch(ecg(:,1),60);
        if NeedFiltering
            ecgn = filtfilt(bn_60Hz,an_60Hz,ecg);
        else
            ecgn = ecg;
        end
        % == remove 50Hz?
        NeedFiltering = check_notch(ecg(:,1),50);
        if NeedFiltering
            ecgn = filtfilt(bn_50Hz,an_50Hz,ecgn);
        end
    else
        % == Notch fileter at 50Hz AND 60Hz
        ecgn = filtfilt(bn_60Hz,an_60Hz,ecg);
        ecgn = filtfilt(bn_50Hz,an_50Hz,ecgn); 
    end

    % == remove high freq and baseline
    ecgnb  = ecgn-filtfilt(b_lp,a_lp,ecgn);
    ecgnbh = filtfilt(b_bas,a_bas,ecgnb);

    % == normalize the data
    [ecgf,~] = normalise_ecg(ecgnbh,1*paramStruct.fs,5*paramStruct.fs);
    % reason for starting at 1sec is because of filtering initialisation at the
    % border which takes some datapoints

    ecgf(isnan(ecgf))=0; %FIXME
    
catch ME
    for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
    ecgf = ecg;
end

% == plots
if debug
    % plot first channel before and after preprocessing
    tm = 1/paramStruct.fs:1/paramStruct.fs:length(ecg)/paramStruct.fs; figure(1);
    ax(1) = subplot(211); plot(tm,ecg(:,1)); legend('before filtering');
    ax(2) = subplot(212); plot(tm,ecgf(:,1),'r','LineWidth',2); 
    legend('after filtering'); linkaxes(ax,'x'); xlabel('Time [sec]'); ylabel('Amplitude (NU)');
    set(findall(gcf,'type','text'),'fontSize',14,'fontWeight','bold');
end
end


%%Plot PSD       
%     plot_psd(ecg(:,1),ecgnbh(:,1),paramStruct.fs,'welch');
%     global NameCurrent
%     saveas(gcf,NameCurrent(1:3));