% daisy2ref() - Re-reference data acquired with the bipolar daisy chain to 
%               a common-reference or average reference EEG montage. 
%
% Usage:
%   >>  newdata = daisy2ref( data);
%   >>  [newdata chlab] = daisy2ref( data, 'key' , 'val', ... );
%
% Inputs:
%   data     -  2D data matrix ( n.chans. x time )
%
% Optional inputs:
%   'ref'      - ['average'|'common'] convert to average or common reference. 
%                Default is 'common'.
%   'refelec'  - [integer] electrode to use as common reference (as from channel list). 
%                Default is 'P10', or 2.
%   'badchan'  - [integer array] indices of daisy-sequencial bad channels.
%   'chanlocs' - ['string'|'sccn'] channel electrode location file. 'sccn'
%                use the SCCN daisy chain channel location file. Default is [].
% 
% Output:
%   newdata  - Re-referenced data: newdata = M * data
%   chlab    - channel labels for newdata
%
% Author: Arnaud Delorme, CNL / Salk Institute - SCCN, 25 Oct 2002
%
% See also: 
%   averef(), eeglab() 

%  Copyright (C) 2002 Arnaud Delorme, CNL / Salk Institute - SCCN, La Jolla
%
% 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

% $Log: daisy2ref2.m,v $
% Revision 1.4  2003/08/08 22:24:50  arno
% *** empty log message ***
%
% Revision 1.3  2003/04/25 23:51:05  arno
% quit
%
%
%
%
% exit
%
% Revision 1.2  2003/03/03 00:02:41  arno
% default location file
%
% Revision 1.1  2003/03/02 23:50:54  arno
% Initial revision

function [newdata, chanlocs] = daisy2ref( data, varargin );

disp('WARNING: daisy2ref2 is not part of the EEGLAB toolbox and should not be distributed');
disp('         you must contact Arnaud Delorme (arno@salk.edu) for terms of use');

if nargin < 1
	help daisy2ref;
	return;
end;	

g = finputcheck(varargin, { 'ref'      'string'   { 'average' 'common' }    'common'; ...
                            'chanlocs' 'string'   []                        ''; ...
                            'refelec'  'integer'  [1 72]                    2; ...
                            'badchan'  'integer'  [1 72]                    []; ...
                            'badelec'  'integer'  [1 72]                    [] }, 'daisy2ref');

if isstr(g), error(g); end;
if ~isempty(g.badelec),
    g.badchan = g.badelec;
end;
if size( data, 1) > 72
    data = dat(1:72,:);
    disp('daisy2ref: more than 72 channels, removing last channels');
end;
if size( data, 1) < 72
    % artificially add artifactual channels
    tmpdata = zeros(72, size(data,2));
    goodelec = setdiff(1:72, g.badchan);
    tmpdata(goodelec,:) = data;
    data = tmpdata;
end;

A = load('/data/common/projects/eegmri/daisymatrix.txt');

% remove the bad channels rows
% ----------------------------
if any(g.badchan < 0)
    A(-g.badchan,:)    = [];
else
    A(g.badchan,:)    = [];
    data(g.badchan,:) = [];
end;

% scan for channels that can not be estimated
% -------------------------------------------
badelec = find(sum(abs(A),1) == 0);
fprintf('Electrodes [%s] can not be estimated\n', int2str(badelec));
A(:, badelec) = [];
if ~isempty(g.chanlocs)
    if strcmpi(g.chanlocs, 'sccn'), 
        g.chanlocs = '/data/common/matlab/miscfunc/daisy72.locs';
    end;
    chanlocs = readlocs(g.chanlocs);
    chanlocs(:,badelec) = [];
else
    chanlocs = [];
end;

% add some average reference channels
% -----------------------------------
if rank(A) > size(A,2)-1
    fprintf('Low dimension matrix, some electrode projection will be innacurate\n');
end;
if strcmp(g.ref, 'average') 
    fprintf('Adding average reference channel\n');
    A(end+1,:)    = ones(1,size(A,2)); % add an average reference channel
    data(end+1,:) = 0;                 % add a zeros channel
end;

% solve
% -----
newdata = pinv(A)*data; % solve up to some multiplicative constant (common to 
                        % all channels) if the matrix is not invertible 

if ~strcmp(g.ref, 'average')
    if ~isempty ( find(g.refelec == g.badchan) )
        help daisy2ref;
        error('daisy2ref(): the reference channel cannot be a bad electrode.')
    end
    fprintf('Computing common reference with channel %d\n', g.refelec);
    newdata = newdata - repmat(newdata(g.refelec,:), [size(newdata,1) 1]); % re-reference
    newdata(g.refelec,:) = [];
    chanlocs(g.refelec) = [];
else
    fprintf('Recomputing average reference\n');
    newdata = newdata - repmat(mean(newdata(:,:),1), [size(newdata,1) 1]); % re-reference
end;

% see test code at the end
return

% *****************
%        END
% *****************

% OLD VERSION - DOING MORE COMPLEX STUFF BUT PINV (above) SOLVE THEM ALL
% ----------------------------------------------------------------------
% we need to solve the system A*x = data (x = inv(A) * data)
%
%                 electrode dimension
%A     =       1  0  0  0  0  0  0  0 -1 ...
%  channel     0  0  0  0  0  0 -1  0  0 ...
%  dimension   0  0  1  0  0  0  0 -1  0 ...
%              ....

% subtract the reference electrode from every channel
%      [x(2)]                                             [0 -1  0 ...]
% B.(x-[x(2)])  =  A.x   of B.C.x = A.x with C = eye(n) - [0 -1  0 ...] (for ref in second electrode)
%      [x(2)]                                             [0 -1  0 ...]  
%      ...                                                ...
% so B = A.C^-1

C = zeros(size(A));
C(g.refelec,:) = 1;
C = eye(size(A)) - C; 
A = A*pinv(C); % problem: we loose one electrode here
%A(:,g.refelec) = []; % not necessary already removed

% add some average reference channels
% -----------------------------------
if size(A,1) > size(A,2) % overcomplete case
    fprintf('Overcomplete system, removing last channel\n');
    A(end,:) = []; % remove last channel but then have to remove every channel
                   % sequentially and resolve the system
    data(end,:) = [];
elseif size(A,2) > size(A,1)
    fprintf('Non-singular system, adding average reference channels\n');
    for index = 1:size(A,2)-size(A,1)
        A(end+1,:) = [zeros(1,index) ones(1,size(A,2)-index)]; % add an average reference channel
        data(end+1,:) = 0;                                     % add a zeros channel
    end;
end;

% perform the conversion
% ----------------------
invA = pinv(A);
newdata = invA * data;

% put the reference channel to 0
% ------------------------------
newdata(g.refelec,:) = 0;
return;

% testing
% -------
data = rand(72, 1000);
badchans = [8];
a = daisy2ref(data,'refelec', 2, 'badchan',  badchans);
c = daisy2ref4(data,'refelec', 2, 'badchan', badchans);
(c(1:end,1)-a(1:size(c,1),1))'