Added robust re-referencing.

Author: chkothe

code/filters/flt_clean_windows.m

@@ -23,6 +23,7 @@
 %                    output) to 0.3 (very lax cleaning of only coarse artifacts). Default: 0.15.
 %
 %
+%
 %   The following are detail parameters that usually do not have to be tuned. If you can't get
 %   the function to do what you want, you might consider adapting these to your data.
 %
@@ -56,6 +57,7 @@
 %                     is quite slow. (default: true)
 %
 %
+%
 %	The following are legacy parameters to enable previous behaviors.
 %
 %   FlaggedQuantile : upper quantile of the per-channel windows that should be flagged for potential
@@ -185,7 +187,7 @@
         % previous method (2011): based on quantiles
         removed_windows = find(sum(1-retain_mask) > min_badchans);
     elseif ~isempty(ignored_chans)
-        % oldest method (2010)
+        % oldest method (2010): attempt at being robust against broken channels
         removed_windows = find(sum(retain_mask) <= ignored_chans);
     else
         error('BCILAB:flt_clean_windows:legacy_options','By setting the flag_quantile parameter you switch to legacy methods; these require that you also set the min_badchans (former default: 0.5) or ignored_chans (former default: 0.3) parameter.');

code/filters/flt_reref.m

@@ -1,6 +1,6 @@
 function signal = flt_reref(varargin)
 % Re-references the data to a new (set of) channel(s) or the average of all channels.
-% Signal = flt_reref(Signal, ReferenceChannels, ExcludeChannels, KeepReference)
+% Signal = flt_reref(Signal, ReferenceChannels, ExcludeChannels, KeepReference, ReferenceType, HuberCutoff, HuberIterations)
 %
 % Re-referencing is a spatial filter in which the (instantaneous) mean signal of some "reference"
 % channels is subtracted from the signal of each channel. This allows to dampen externally-induced
@@ -12,6 +12,9 @@
 % re-referencing, since their spatial filters are adaptively optimized and effectively incorporate
 % re-referencing to the degree that it is necessary.
 %
+% This function also supports robust re-referencing, where artifacts in few channels will not affect
+% the outcome disproportionally.
+%
 % In:
 %   Signal            :  Epoched or continuous data set.
 %
@@ -22,6 +25,19 @@
 %                       (default: [])
 %
 %   KeepReference     : Whether to keep the reference channel (default: false)
+%                       If true and multiple channels were selected, a new channel named CAR is
+%                       appended.
+%
+%   ReferenceType     : Type of referencing to apply, can be one of the following:
+%                       * 'mean': subtract average of reference channels (default)
+%                       * 'median': subtract median of reference channels
+%                       * 'huber': subtract robust mean (under the Huber loss) of reference channels;
+%
+%   HuberCutoff       : cutoff/transition parameter for huber loss; if [], this is set
+%                       to one (robust) standard deviation of the signal
+%
+%   HuberIterations   : number of iterations to compute the Huber fit; a larger number can tolerate
+%                       larger outliers (default: 100)
 %
 % Out:
 %   Signal : Re-referenced data set.
@@ -45,18 +61,52 @@
 %                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
 %                                2010-03-28
 
-% flt_reref_version<1.0> -- for the cache
+% flt_reref_version<1.01> -- for the cache
 
 if ~exp_beginfun('filter') return; end
 
 declare_properties('name',{'Rereferencing','ref'}, 'independent_channels',false,'independent_trials',true);
 
 arg_define(varargin,...
     arg_norep({'signal','Signal'}), ...
-    arg({'chn','ReferenceChannels'}, [], [], 'Cell array of reference channels. The signal data is be referenced to these, defaults to average reference if empty.','type','cellstr','shape','row'),...
-    arg({'exclude','ExcludeChannels'}, [], [], 'Cell array of channels to exclude.','type','cellstr','shape','row'),...
-    arg({'keepref','KeepReference'}, false, [], 'Keep the reference channel.'));
+    arg({'ref_chn','ReferenceChannels','chn'}, [], [], 'Cell array of reference channels. The signal data is be referenced to these, defaults to average reference if empty.','type','cellstr','shape','row'),...
+    arg({'exclude_chn','ExcludeChannels','exclude'}, [], [], 'Cell array of channels to exclude.','type','cellstr','shape','row'),...
+    arg({'keepref','KeepReference'}, false, [], 'Keep the reference channel.'), ...
+    arg({'ref_type','ReferenceType'}, 'mean', {'mean','median','huber'}, 'Type of reference. The traditional average reference operation uses the mean. If this is set to median, the median of the reference channels will be removed, and if set to huber the robust mean under the Huber loss will be removed (slower than median but closer to the mean).'), ...
+    arg({'huber_cut','HuberCutoff'}, [], [], 'Cutoff for huber function. If left empty this is set to one (robust) standard deviation of the signal.','guru',true), ...
+    arg({'huber_iters','HuberIterations'}, 100, [], 'Iterations for huber fitting. A larger number yields tolerance to larger outliers.','guru',true));
+
+if ~isempty(ref_chn)
+    ref_chns = set_chanid(signal,ref_chn);
+else
+    ref_chns = 1:size(signal.data,1);
+end
+    
+if ~isempty(exclude_chn)
+    ref_chns = setdiff(ref_chns,set_chanid(signal,exclude_chn)); end
+
+switch ref_type
+    case 'mean'
+        ref_signal = mean(signal.data(ref_chns,:));
+    case 'median'
+        ref_signal = median(signal.data(ref_chns,:));
+    case 'huber'
+        if isempty(huber_cut)
+            huber_cut = median(mad(signal.data,1,2))*1.4826; end
+        ref_signal = robust_mean(signal.data(ref_chns,:)/huber_cut,1,huber_iters)*huber_cut;
+    otherwise
+        error('Unsupported reference type.');
+end
+    
+signal.data = bsxfun(@minus,signal.data,ref_signal);
 
-signal = pop_reref(signal,set_chanid(signal,chn),'exclude',fastif(isempty(exclude),[],set_chanid(signal,exclude)),'keepref',fastif(keepref,'on','off')); %#ok<*NODEF>
+if ~keepref && ~isempty(ref_chn)
+    signal = pop_select(signal,'nochannel',ref_chns); end
+if keepref && length(ref_chn)>1
+    signal.data(end+1,:) = ref_signal;
+    signal.nbchan = size(signal.data,1);
+    signal.chanlocs(end+1).labels = 'CAR';
+    signal.chanlocs(end).type = 'REF';
+end
 
-exp_endfun;
+exp_endfun('append_online',{'huber_cut',huber_cut});

code/helpers/hlp_collect_datasets.m

@@ -49,9 +49,6 @@
 %       'conditions',@(EEG) ~isempty(EEG.icawinv) && mean(cellfun('isempty',{EEG.chanlocs.X})) < 0.5, ...
 %       'collect',@(EEG,path){path,EEG.icawinv,EEG.chanlocs,EEG.icachansind})
 %
-% Dependencies: 
-%   env_translatepath.
-%
 %                               Laura Froelich, Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
 %                               2010-09-10
 
@@ -75,7 +72,7 @@
     opts = varargin{1};                                                     % opts is already a struct: fast path
 else
     opts = hlp_varargin2struct(varargin, 'pattern',regexptranslate('wildcard','*.set'), 'nopattern',{}, 'checkset',1, 'nowarnings',1 ,'nodialogs',1, ...
-        'maxsize',5*2^20, 'maxtime',Inf, 'maxnumber',Inf, 'conditions',[], 'collect',@(EEG,path) path, ....
+        'maxsize',5*2^20, 'maxtime',Inf, 'maxnumber',Inf, 'conditions',{}, 'collect',@(EEG,path) path, ....
         'fileconditions',@(path,name,ext) exist([path filesep  name '.fdt'],'file') || exist([path filesep name '.dat'],'file'));
     if ~iscell(opts.conditions)
         opts.conditions = {opts.conditions}; end

code/misc/robust_fit.m

@@ -0,0 +1,37 @@
+function x = robust_fit(A,y,rho,iters)
+% Perform a robust linear regression using the Huber loss function.
+% x = robust_fit(A,y,rho,iters)
+%
+% Input:
+%   A : design matrix
+%   y : target variable
+%   rho : augmented Lagrangian variable (default: 1)
+%   iters : number of iterations to perform (default: 1000)
+%
+% Output:
+%   x : solution for x
+%
+% Notes:
+%   solves the following problem via ADMM for x:
+%     minimize 1/2*sum(huber(A*x - y))
+%
+% Based on the ADMM Matlab codes also found at:
+%   http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
+%
+%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
+%                                2013-03-04
+
+if ~exist('rho','var')
+    rho = 1; end
+if ~exist('iters','var')
+    iters = 1000; end
+
+Aty = A'*y;
+L = sparse(chol(A'*A,'lower')); U = L';
+z = zeros(size(y)); u = z;
+for k = 1:iters
+    x = U \ (L \ (Aty + A'*(z - u)));
+    d = A*x - y + u;
+    z = rho/(1+rho)*d + 1/(1+rho)*max(0,(1-(1+1/rho)./abs(d))).*d;
+    u = d - z;
+end

code/misc/robust_mean.m

@@ -0,0 +1,36 @@
+function X = robust_mean(Y,rho,iters)
+% Perform a robust mean under the Huber loss function.
+% x = robust_mean(Y,rho,iters)
+%
+% Input:
+%   Y : MxN matrix over which to average (columnwise)
+%   rho : augmented Lagrangian variable (default: 1)
+%   iters : number of iterations to perform (default: 1000)
+%
+% Output:
+%   x : 1xN vector that is the roust mean of Y
+%
+% Based on the ADMM Matlab codes also found at:
+%   http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
+%
+%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
+%                                2013-09-26
+
+if ~exist('rho','var')
+    rho = 1; end
+if ~exist('iters','var')
+    iters = 1000; end
+
+m = size(Y,1);
+if m==1
+    X = Y;
+else
+    mu = sum(Y)/m;
+    Z = zeros(size(Y)); U = Z;
+    for k = 1:iters
+        X = mu + sum(Z - U)/m;
+        D = bsxfun(@minus,X,Y - U);
+        Z = (rho/(1+rho) + (1/(1+rho))*max(0,(1-(1+1/rho)./abs(D)))).*D;
+        U = D - Z;
+    end
+end