par_autobss.m

function [dataOut] = par_autobss(dataIn, opt)
% autobss() - Performs automatic EOG artifact correction using Blind Source
% Separation (BSS) and identifying the EOG components using fractal analysis
%
% Usage:
%   >>  Y = autobss( X, opt)
%
% Inputs:
%   X              - Input data matrix, d x L
%   opt.wl         - Length of the moving-ICA analysis windows (in samples).
%                    If empty it will be set to the data length L.
%                    default: L
%   opt.ws         - Shift between correlative windows (in samples). If
%                    empty it will be set to the same value as wl.
%                    default: wl
%   opt.wl_dim     - window length (in samples) for the moving window
%                    computation of the fractal dimension
%                    default: .1*opt.wl
%   opt.ws_dim     - window shift (in samples) for the moving window
%                    computation of the fractal dimension
%                    default: opt.wl_dim
%   opt.bss_alg    - BSS algorithm to use
%                    default 'sobi'
%   opt.bss_opt    - Options to pass to the BSS algorithm
%                    default: []
%   opt.crit_alg   - criterion to use for detecting the EOG components
%                    ('fd','svf','joyce'), default: 'fd'
%   opt.crit_opt   - options to pass to the criterion function
%                    default: []
%   opt.tau        - embedding lag for computing the 'svf' criterion
%                    default: 1
%   opt.dim        - embedding dimension for computing the 'svf' criterion
%                    default: 20
%   opt.k          - index of the SVF for the 'svf' criterion
%                    default: 1
%   opt.eogindex   - indexes of the EOG channels. Necessary for computing
%                    'joyce' criterion. 
%                    default: last data channel
%
% Outputs:
%   Y   - Output data matrix (artifact corrected)
%
% Notes:
%   1) BSS will be performed on (possibly overlapping) windows of wl
%      samples.
%   2) The EOG components are detected as described in [3]. Please refer to
%      [3] when using this code in any of your publications.
%   3) Reconciling ovelapping analysis windows is done as in [1]. 
%   4) Available BSS algorithms should have an associated interface
%      function named [cmd]_ifc where [cmd] is the command name used to run
%      the BSS. A sample interface to sobi (named sobi_ifc) is included in
%      folder bss_alg.
%   5) The SVF criterion is described in [4] (not implemented yet!).
%   6) The Joyce criterion is described in [5] (not implemented yet!).
%
% References:
% [1] Wallstrom et al., G.L. International Journal of Psychology 53 (2004)
%     105-119
% [2] Katz, M., Comput. Biol. Med. 18 (1988), 145-156
% [3] Gomez-Herrero, G. , De Clercq, W., Anwar, H., Kara, O., Egiazarian, K. (2006),
%     Proceedings of NORSIG 2006, Reykjavik, Iceland.
% [4] Faul, S., Marnane, L., Lightbody, G., Boylan, G. and Connolly, S. (2005),
%     Proceedings of ICASSP 2005, Philadelphia, USA.
% [5] Joyce, C. A., Gorodnitsky, I.F., Kutas, M. (2004), Psychophysiology,
%     41, 313-325
%
%
% Author: German Gomez-Herrero <german.gomezherrero@gmail.com>
%         Institute of Signal Processing
%         Tampere University of Technology, 2007
%
% See also:
%   POP_AUTOBSSEOG, POP_AUTOBSSEMG, CMERGE_OVERLAP, CMERGE_NOOVERLAP, FD, EEGLAB
%

% Copyright (C) <2007>  German Gomez-Herrero, http://germangh.com


if nargin < 1, help autoica; return; end

if ~exist('opt','var')
    opt = def_autobss;
else
    opt = def_autobss(opt);
end

bss_alg = lower(opt.bss_alg);
ws = opt.ws;
wl = opt.wl;

% Initial and final sample of each analysis epoch
% ----------------------------------------
[d,L] = size(dataIn);

if isempty(wl)
    wl = L; ws=L;
elseif isempty(ws)
    ws = wl;
end


% overlap between correlative samples
% -----------------------------------------
ovlength = wl-ws;

init = 1:ws:L;
final = init+wl-1;
samples = length(init);

% initialize the output (corrected)data
% ----------------------------------------
dataOut = zeros(d,L);

% reference EOG channels data
% ----------------------------------------
if isfield(opt.crit_opt,'eogref')
    eogref = opt.crit_opt.eogref;
else
    eogref = [];
end
finals = cell(1, samples); % Pre-allocate cell array for performance
[finals{:}] = deal(final); % Distribute 'final' to all elements of 'finals'
%finals{samples} = deal(final);

% repositioning of OG code; remove excess samples
for i = 1:samples
    final1 = finals(i);
    final1 = final1{:};
    if final1(i)>L
        if (final1(i)-L)>.5*wl
            samples=i-1;
            %warning('Last analysis window is too short and therefore will not be processed');
            break; 
        else
            final1(i)=L;
            if i < samples
                final1(i+1) = Inf;
            end
        end
    end
    finals{i} = final1;
end


% assigning all variables before loop starts

X{samples} = deal(dataIn);
Xi{samples} = deal([]);
Yis{samples} = deal([]);
Ais{samples} = deal([]);
Wis{samples} = deal([]);
finals{samples} = deal(final);
dataOuts{samples} = deal(dataOut);
index{samples} = deal([]);
opts{samples} = deal(opt);
Ximeans{samples} = deal([]);
thisYs{samples} = deal([]);
thisXs{samples} = deal([]);
EOGs{samples} = deal([]);


if opt.verbose,fprintf('\nRunning BSS filter in sliding windows in PAR');end

%% TEST IF HAS PAR
% RUN IF STATEMENT if par etc.
parfor i = 1:samples

    X = dataIn;
    final = finals{i};
    dataOuti = dataOut;


    Xi{i} = X(:,init(i):final(i));
    % remove data mean
    Ximean = mean(Xi{i},2);
    Xi{i} = Xi{i}-repmat(Ximean,1,size(Xi{i},2));
    if ~isempty(eogref)
        eogrefi = eogref(:,init(i):final(i));
        opts(i).crit_opt.eogref = eogrefi;
    end

    if isfield(opt,'bss_opt') && ~isempty(opt.bss_opt)
        %F = str2fun( strcat( bss_alg ,'_ifc(', DataIni(i),  opt.bss_opt));
        % Convert function name string to a function handle
        bss_func_handle = str2func([bss_alg '_ifc']);
        
        % Run the function using feval
        Wi = feval(bss_func_handle, Xi{i}, opt.bss_opt);
        %[Wi] = feval([bss_alg '_ifc(Xi{i},opt.bss_opt)']);

        Wi = real(Wi);
    else
        % Convert function name string to a function handle
        bss_func_handle = str2func([bss_alg '_ifc']);
        
        % Run the function using feval
        Wi = feval(bss_func_handle, Xi{i});
        %[Wi] = feval([bss_alg '_ifc(Xi{i})']);

        Wi = real(Wi);
    end
    if ~isempty(Wi)
        Ai = pinv(Wi);
        Yi = Wi*Xi{i};
    else
        Ai = [];
        Yi = [];
    end

    % detect the EOG-related components
    switch lower(opt.crit_alg)
        case 'eog_fd'   % fractal dimension
            [index{i}] = eog_fd(Yi,opt.crit_opt);
        case 'eog_corr' % EOG correlation
            [index{i}] = eog_corr(Yi,opt.crit_opt);
        case 'emg_psd'  % PSD ratio EEGband/EMGband
            [index{i}] = emg_psd(Yi,opt.crit_opt);
        case 'eog_svf' % singular value fraction
            error('(par_autobss) criterion eog_svf not working yet. Try fd instead.') 
        case 'eog_joyce' % joyce criterion
            error('(par_autobss) criterion eog_joyce not working yet. Try eog_fd instead.')
    end    

    % correct the EOG (avoiding unnecessary operations and rounding errors)
    if wl>ws
        % if there is overlap between correlative analysis windows
        if isempty(index{i})       % leave the data unchanged
            thisY = X(:,init(i):final(i));        
        elseif length(index{i})<d   % remove just some components
            EOG = real(Ai(:,index{i})*Yi(index{i},:));
            thisY = (Xi{i}-EOG)+repmat(Ximean,1,size(Xi{i},2));
        else                     % remove all data
            thisY = zeros(d,wl);
        end        
        
        if i < 2
            dataOuti(:,1:final(1))=thisY;
        else
            dataOuti(:,1:final(i)) = cmerge_overlap(dataOuti(:,1:final(i-1)),thisY,ovlength);
        end
    else

        % if there is no overlap
        if isempty(index{i})       % leave the data unchanged
            dataOuti(:,init(i):final(i)) = X(:,init(i):final(i));        
        elseif length(index{i})<d   % remove just some components
            EOG = real(Ai(:,index{i})*Yi(index{i},:));
            dataOuti(:,init(i):final(i)) = (Xi{i}-EOG)+repmat(Ximean,1,size(Xi{i},2));
        end                
    end

    % redistribute outputs
    Ais{i} = Ai; Yis{i} = Yi; Wis{i} = Wi; try thisYs{i} = thisY; catch; end;

    dataOuts{i} = dataOuti; % save the output
    Ximeans{i} = Ximean;
    if opt.verbose && ~mod(i,floor(samples/10)), fprintf('.'); end

end % end of parfor loop


% merge all samples


if opt.verbose,fprintf('\nDone.\n');end

% leave the remaining data unchanged 
if final(samples)<L
    dataOut(:,final(samples)+1:end) = dataIn(:,final(samples)+1:end);
end




function [opt] = def_autobss(opt)
% def_autobss() - Sets default analysis parameters for autobss()
%
% Usage:
%   [opt] = def_autobss(opt_in)
%
% Inputs:
%   opt_in          - Input structure containing analysis options
%
% Outputs:
%   opt             - Output structure where all the analysis options not specified
%                     in opt_in have been set to their default values. The possible
%                     analysis options are:
%   opt.wl          - Length (in samples) of the analysis windows. If left
%                     empty it will be set to the data length.
%                     (default: [])
%   opt.ws          - Shift (in samples) between correlative analysis
%                     windows. If empty, the same value as wl will be
%                     assumed.
%                     (default: [])
%   opt.wl_dim      - Length (in samples) of the analysis windows used for
%                     estimating the mean fractal dimension.
%                     (default: [])
%   opt.ws_dim      - Shift (in samples) between correlative windows used for
%                     estimating the mean fractal dimension.
%                     (default: [])
%   opt.bss_alg     - BSS algorithm to use (default: 'sobi')
%   opt.bss_opt     - BSS options (default: [])
%   opt.crit_alg    - criterion to use for identifying the components to be
%                     removed. The following are available:
%                     'eog_fd','eog_corr' -> good for removing EOG
%                     'emg_psd'           -> good for removing EMG
%                     default: 'eog_fd'
%   opt.crit_opt    - Criterion options (default: [])
%   opt.tau         - embedding lag for computing the 'svf' criterion,
%                     default: 1
%   opt.dim         - embedding dimension for computing the 'svf' criterion
%                     default: 20
%   opt.k           - index of the SVF for the 'svf' criterion
%                     default: 1
%   opt.eogindex    - indexes of the EOG channels. Necessary for computing
%                     some of the criteria. 
%                     default: []
%
% Author: German Gomez-Herrero <german.gomezherrero@gmail.com>
%         Institute of Signal Processing
%         Tampere University of Technology, 2006
%
% See also:
%   POP_AUTOBSS, CMERGE_OVERLAP, CMERGE_NOOVERLAP, FD, EEGLAB
%

% Copyright (C) <2006>  <German Gomez-Herrero>
%
% 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.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


if ~exist('opt','var') || ~isfield(opt, 'ws'),
    opt.ws = [];
end
if ~isfield(opt, 'wl'),
    opt.wl = [];
end
if ~isfield(opt, 'wl_dim'),
    opt.wl_dim = [];
end
if ~isfield(opt, 'ws_dim'),
    opt.ws_dim = [];
end
if ~isfield(opt, 'crit_alg') || isempty(opt.crit_alg),
    opt.crit_alg = 'eog_fd';
end
if ~isfield(opt, 'crit_opt') || isempty(opt.crit_opt),
    opt.crit_opt = [];
end
if ~isfield(opt, 'tau') || isempty(opt.tau),
    opt.tau = 1;
end
if ~isfield(opt, 'dim') || isempty(opt.dim),
    opt.dim = 20;
end
if ~isfield(opt,'k') || isempty(opt.k),
    opt.k = 1;
end
if ~isfield(opt, 'verbose') || isempty(opt.verbose),
    opt.verbose = 1;
end
if ~isfield(opt,'bss_alg') || isempty(opt.bss_alg),
    opt.bss_alg = 'sobi';
end
if ~isfield(opt,'bss_opt') || isempty(opt.bss_opt),
    opt.bss_opt = [];
end