function ecgout = bss_step(ecg,icaecgf,rMQRS,paramStruct,debug)
% Blind Source Separation (BSS) step. In the case of the DEFLATION method
% BSS is applied to the ABD signals and MECG cancellation takes place in
% the source domain. In the other cases then BSS is applied on the ecg
% matrix containing the residuals of the MECG cancellation perfomed in the
% time domain.
%
% inputs
%   ecg:         ecg matrix onto which to perform BSS
%   icaecgf:     ica applied on filtered ecg stack
%   rMQRS:       reference MQRS location in samples
%   paramStruct: important algo parameters
%
% output
%   ecgout:      output ecg that came through the BSS step
%                case DEFLATION: [ICA-TS], [ICA-TS,ICA-TS-ICA]
%                case TS: [TS-ICA]
%
%
% 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<4; error('bss_step: wrong number of input arguments \n'); end;
if size(ecg,2)<size(ecg,1); ecg = ecg'; end;
if size(icaecgf,2)>size(icaecgf,1); icaecgf = icaecgf'; end;
if nargin<5; debug=0; end;

% == general
NB_ABD_CHANNELS = size(icaecgf,2);
MQRS = cell(NB_ABD_CHANNELS,1);

% == core function
try
    if isempty(find(strcmp(paramStruct.method,{'DEFLATION','FUSE'}),1))
        % == ICA on residual in time domain
        [~,icaecgr] = jade(ecg);
        ecgout = icaecgr;
    end

    if find(strcmp(paramStruct.method,{'DEFLATION','FUSE'}))
        % == ICA on residual in source space
        recgs = zeros(size(icaecgf')); % residual ecg in source domain
        for i=1:NB_ABD_CHANNELS
            MQRS{i} = adjust_mqrs_location(icaecgf(:,i),rMQRS,paramStruct.fs,0);
            recgs(i,:) = mecg_cancellation(MQRS{i},icaecgf(:,i)','TS',20);
            %[recgs(i,:),~,~] = run_kalman(icaecgf(:,i),MQRS{i},paramStruct);
        end
        [~,icaecgr] = jade(recgs);
        if strcmp(paramStruct.methodComplexity,'COMPLEX')
            ecgout = [recgs;icaecgr]; 
        else
            ecgout = recgs;
        end
    end
    
    if size(ecg,1)>size(ecg,2)
        ecgout = ecgout';
    end
    
catch ME
    for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
    ecgout = ecg;
end

% == plots
if debug
   ax(1) = subplot(211); plot(ecg); 
   title('ecgin');
   ax(2) = subplot(212); plot(ecgout);
   title('ecgout');
   linkaxes(ax,'x');
end

end

% SOME NOTES ON THIS FUNCTION:
% == STEP 1: direct BSS on filtered ecg (and not residuals)
    % Adding another step where BSS is performed on the filtered ABD signal was usefull 
    % for the EKF because sometimes the EKF kills the FECG and applying BSS after removing 
    % the MECG does not work really well.
% == STEP 2: ICA on residual in time domain
    % NOTE: jade performed better than PCA and FastICA (default options) when applied to the residual
    % signal but PCA performed better than ICA when applied to the initial
    % prefiltered mixture. Difference in performance was >1% for F1
    % measure but this is due to one specific record that is not separated
    % by ICA but is by ICA. Removing this specific record results in an
    % increase in F1 when using ICA rather than PCA on the ecgf
    % matrix. Having more data in the training set would have been usefull
    % to confirm this observation.
% == STEP 3: ICA on residual in source space    
    % ICA to move to the space domain, performing MECG cancellation in that
    % domain then take ICA of mother-free source signals. 
    % So it is sort of ICA of ICA...   

% == NOTE: in theory ICA would need centering that is subtracting the mean of the
% signal for each channel before applying ICA. Performing this operation
% lowered the F1 measure slightly!
% [~,icaecgr] = jade(bsxfun(@minus,recgs,mean(recgs,2)));