function [data] = ft_combineplanar(cfg, data)

% FT_COMBINEPLANAR computes the planar gradient magnitude over both directions

% combining the two gradients at each sensor to a single positive-valued number. This

% can be done for single-trial/averaged planar gradient ERFs or single-trial/averaged

% TFRs.

%

% Use as

% [data] = ft_combineplanar(cfg, data)

% where data contains an averaged planar-gradient ERF or single-trial or

% averaged TFRs.

%

% The configuration can contain

% cfg.method = 'sum', 'svd', 'abssvd', or 'complex' (default = 'sum')

% cfg.updatesens = 'no' or 'yes' (default = 'yes')

% and for timelocked input data (i.e. ERFs), the configuration can also contain

% cfg.demean = 'yes' or 'no' (default = 'no')

% cfg.baselinewindow = [begin end]

%

% To facilitate data-handling and distributed computing you can use

% cfg.inputfile = ...

% cfg.outputfile = ...

% If you specify one of these (or both) the input data will be read from a *.mat

% file on disk and/or the output data will be written to a *.mat file. These mat

% files should contain only a single variable, corresponding with the

% input/output structure.

%

% See also FT_MEGPLANAR

% Undocumented local options:

% cfg.foilim

% cfg.trials

% Copyright (C) 2004, Ole Jensen

% Copyright (C) 2004-2013, Robert Oostenveld

%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org

% for the documentation and details.

%

% FieldTrip 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 3 of the License, or

% (at your option) any later version.

%

% FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.

%

% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts

ft_revision = '$Id$';

ft_nargin = nargin;

ft_nargout = nargout;

% do the general setup of the function

ft_defaults

ft_preamble init

ft_preamble debug

ft_preamble loadvar data

ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init

if ft_abort

return

end

% check if the input data is valid for this function

data = ft_checkdata(data, 'datatype', {'raw', 'freq', 'timelock'}, 'feedback', 'yes', 'senstype', {'ctf64_planar', 'ctf151_planar', 'ctf275_planar', 'neuromag122', 'neuromag306', 'bti148_planar', 'bti248_planar', 'itab153_planar', 'yokogawa64_planar', 'yokogawa160_planar', 'yokogawa208_planar', 'yokogawa440_planar'});

% check if the input cfg is valid for this function

cfg = ft_checkconfig(cfg, 'forbidden', {'combinegrad'});

cfg = ft_checkconfig(cfg, 'deprecated', {'baseline'});

cfg = ft_checkconfig(cfg, 'renamed', {'blc', 'demean'});

cfg = ft_checkconfig(cfg, 'renamed', {'blcwindow', 'baselinewindow'});

cfg = ft_checkconfig(cfg, 'renamed', {'combinemethod', 'method'});

% set the defaults

cfg.demean = ft_getopt(cfg, 'demean', 'no');

cfg.foilim = ft_getopt(cfg, 'foilim', [-inf inf]);

cfg.baselinewindow = ft_getopt(cfg, 'baselinewindow', [-inf inf]);

cfg.trials = ft_getopt(cfg, 'trials', 'all', 1);

cfg.feedback = ft_getopt(cfg, 'feedback', 'none');

cfg.method = ft_getopt(cfg, 'method', 'sum');

cfg.updatesens = ft_getopt(cfg, 'updatesens', 'yes');

if isfield(cfg, 'baseline')

ft_warning('only supporting cfg.baseline for backwards compatibility, please update your cfg');

cfg.demean = 'yes';

cfg.baselinewindow = cfg.baseline;

end

israw = ft_datatype(data, 'raw');

istimelock = ft_datatype(data, 'timelock');

isfreq = ft_datatype(data, 'freq');

if isfield(data, 'dimord')

dimord = data.dimord;

end

% select trials of interest

if ~strcmp(cfg.trials, 'all')

ft_error('trial selection has not been implemented yet') % first fix ft_checkdata (see above)

end

% find the combination of horizontal and vertical channels that should be combined

planar = ft_senslabel(ft_senstype(data), 'output', 'planarcombined');

[sel_pH, sel_dH] = match_str(planar(:,1), data.label); % indices of the horizontal channels

[sel_pV, sel_dV] = match_str(planar(:,2), data.label); % indices of the vertical channels

% identify and remove unnpaired channels

[dum,iH,iV] = intersect(sel_pH,sel_pV);

sel_dH=sel_dH(iH);

sel_dV=sel_dV(iV);

% find the other channels that are present in the data

sel_other = setdiff(1:length(data.label), [sel_dH(:)' sel_dV(:)']);

lab_other = data.label(sel_other);

% define the channel names after combining the planar combinations

% they should be sorted according to the order of the planar channels in the data

[dum, sel_planar] = match_str(data.label(sel_dH),planar(:,1));

lab_comb = planar(sel_planar,end);

% perform baseline correction

if strcmp(cfg.demean, 'yes')

if ~(istimelock || israw)

ft_error('baseline correction is only supported for timelocked or raw input data')

end

if ischar(cfg.baselinewindow) && strcmp(cfg.baselinewindow, 'all')

cfg.baselinewindow = [-inf inf];

end

% find the timebins corresponding to the baseline interval

tbeg = nearest(data.time, cfg.baselinewindow(1));

tend = nearest(data.time, cfg.baselinewindow(2));

cfg.baselinewindow(1) = data.time(tbeg);

cfg.baselinewindow(2) = data.time(tend);

data.avg = ft_preproc_baselinecorrect(data.avg, tbeg, tend);

end

if isfreq

switch cfg.method

case 'sum'

if isfield(data, 'powspctrm')

% compute the power of each planar channel, by summing the horizontal and vertical gradients

dimtok = tokenize(dimord, '_');

catdim = strmatch('chan',dimtok);

if catdim==1

combined = data.powspctrm(sel_dH,:,:,:) + data.powspctrm(sel_dV,:,:,:);

other = data.powspctrm(sel_other,:,:,:);

elseif catdim==2

combined = data.powspctrm(:,sel_dH,:,:,:) + data.powspctrm(:,sel_dV,:,:,:);

other = data.powspctrm(:,sel_other,:,:,:);

else

ft_error('unsupported dimension order of frequency data');

end

data.powspctrm = cat(catdim, combined, other);

data.label = cat(1, lab_comb(:), lab_other(:));

else

ft_error('cfg.method = ''%s'' only works for frequency data with powspctrm', cfg.method);

end

case 'svd'

if isfield(data, 'fourierspctrm')

fbin = nearest(data.freq, cfg.foilim(1)):nearest(data.freq, cfg.foilim(2));

Nrpt = size(data.fourierspctrm,1);

Nsgn = length(sel_dH);

Nfrq = length(fbin);

Ntim = size(data.fourierspctrm,4);

fourier = nan(Nrpt,Nsgn,Nfrq,Ntim);

ft_progress('init', cfg.feedback, 'computing the svd');

for j = 1:Nsgn

ft_progress(j/Nsgn, 'computing the svd of signal %d/%d\n', j, Nsgn);

for k = 1:Nfrq

fdat = reshape(data.fourierspctrm(:,[sel_dH(j) sel_dV(j)], fbin(k),:), [Nrpt 2 Ntim]);

fdat = permute(fdat, [2 3 1]); % 2 Ntim Nrpt

fdat = reshape(fdat, [2 Ntim*Nrpt]); % 2 Ntim*Nrpt

timbin = ~isnan(fdat(1,:));

[frot, ut, ori, sin_val] = svdfft(fdat(:,timbin), 2, data.cumtapcnt);

dum = nan(Ntim, Nrpt); % Ntim Nrpt

dum(timbin) = frot(1,:); % Ntim Nrpt, insert the first channel of the rotated data

fourier(:,j,k,:) = transpose(dum); % Nrpt Ntim

data.ori{k} = ori; % to change into a cell

data.eta{k} = sin_val(1)/sum(sin_val(2:end)); % to change into a cell

%for m = 1:Ntim

% dum = data.fourierspctrm(:,[sel_dH(j) sel_dV(j)],fbin(k),m);

% timbin = find(~isnan(dum(:,1)));

% [fourier(timbin,j,k,m)] = svdfft(transpose(dum(timbin,:)),1);

%end

end

end

ft_progress('close');

other = data.fourierspctrm(:,sel_other,fbin,:);

data = rmfield(data, 'fourierspctrm');

data.fourierspctrm = cat(2, fourier, other);

data.label = cat(1, lab_comb(:), lab_other(:));

data.freq = data.freq(fbin);

else

ft_error('cfg.method = ''%s'' only works for frequency data with fourierspctrm', cfg.method);

end

otherwise

ft_error('cfg.method = ''%s'' is not supported for frequency data', cfg.method);

end % switch method

elseif (israw || istimelock)

if istimelock

% convert timelock to raw

data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes');

end

switch cfg.method

case 'sum'

Nrpt = length(data.trial);

for k = 1:Nrpt

combined = sqrt(data.trial{k}(sel_dH,:).^2 + data.trial{k}(sel_dV,:).^2);

other = data.trial{k}(sel_other,:);

data.trial{k} = [combined; other];

end

data.label = cat(1, lab_comb(:), lab_other(:));

case 'complex'

Nrpt = length(data.trial);

for k = 1:Nrpt

combined = data.trial{k}(sel_dH,:)*1i + data.trial{k}(sel_dV,:);

other = data.trial{k}(sel_other,:);

data.trial{k} = [combined; other];

end

data.label = cat(1, lab_comb(:), lab_other(:));

case {'svd' 'abssvd'}

Nrpt = length(data.trial);

Nsgn = length(sel_dH);

Nsmp = cellfun('size', data.trial, 2);

Csmp = cumsum([0 Nsmp]);

% do a 'fixed orientation' across all trials approach here

% this is different from the frequency case FIXME

tdat = zeros(2, sum(Nsmp));

for k = 1:Nsgn

for m = 1:Nrpt

tdat(:, (Csmp(m)+1):Csmp(m+1)) = data.trial{m}([sel_dH(k) sel_dV(k)],:);

end

if strcmp(cfg.method, 'abssvd')||strcmp(cfg.method, 'svd')

[rdat, ut, ori, sin_val] = svdfft(tdat, 2);

data.ori{k} = ori; % to change into a cell

data.eta{k} = sin_val(1)/sum(sin_val(2:end)); % to change into a cell

if strcmp(cfg.method, 'abssvd')

rdat = abs(rdat(1,:));

else

rdat = rdat(1,:);

end

end

rdat = mat2cell(rdat, 1, Nsmp);

for m = 1:Nrpt

if k==1, trial{m} = zeros(Nsgn, Nsmp(m)); end

trial{m}(k,:) = rdat{m};

end

end % for each MEG channel

for m = 1:Nrpt

other = data.trial{m}(sel_other,:);

trial{m} = [trial{m}; other];

end

data.trial = trial;

data.label = cat(1, lab_comb(:), lab_other(:));

otherwise

ft_error('cfg.method = ''%s'' is not supported for timelocked or raw data', cfg.method);

end % switch method

if istimelock

% convert raw to timelock

data = ft_checkdata(data, 'datatype', 'timelock', 'feedback', 'yes');

end

else

ft_error('unsupported input data');

end % which ft_datatype

% remove the fields for which the planar gradient could not be combined

data = removefields(data, {'crsspctrm', 'labelcmb'});

if strcmp(cfg.updatesens, 'yes') && isfield(data, 'grad')

% update the grad and only retain the channel related info

[sel_dH, sel_comb] = match_str(data.grad.label, planar(:,1)); % indices of the horizontal channels

[sel_dV ] = match_str(data.grad.label, planar(:,2)); % indices of the vertical channels

% find the other channels that are present in the data

sel_other = setdiff(1:length(data.grad.label), [sel_dH(:)' sel_dV(:)']);

lab_other = data.grad.label(sel_other);

lab_comb = planar(sel_comb,end);

% compute the average position

newpos = [

(data.grad.chanpos(sel_dH,:)+data.grad.chanpos(sel_dV,:))/2

data.grad.chanpos(sel_other,:)

];

% compute the average orientation

newori = [

(data.grad.chanori(sel_dH,:)+data.grad.chanori(sel_dV,:))/2

data.grad.chanori(sel_other,:)

];

newlabel = [

lab_comb

lab_other

];

newtype = [

repmat({'unknown'}, numel(sel_comb), 1) % combined planar

data.grad.chantype(sel_other(:)) % keep the known channel details

];

newunit = [

repmat({'unknown'}, numel(sel_comb), 1) % combined planar

data.grad.chanunit(sel_other(:)) % keep the known channel details

];

newgrad.chanpos = newpos;

newgrad.chanori = newori;

newgrad.label = newlabel;

newgrad.chantype = newtype;

newgrad.chanunit = newunit;

newgrad.unit = data.grad.unit;

newgrad.type = [data.grad.type '_combined'];

% remember the original channel position details

if isfield(data.grad, 'chanposold')

newgrad = copyfields(data.grad, newgrad, {'chanposold', 'chanoriold', 'labelold', 'chantypeold', 'chanunitold'});

else

newgrad.labelold = data.grad.label;

newgrad.chanposold = data.grad.chanpos;

newgrad.chanoriold = data.grad.chanori;

newgrad.chantypeold = data.grad.chantype;

newgrad.chanunitold = data.grad.chanunit;

end

% replace it with the updated gradiometer description

data.grad = newgrad;

end

% convert back to input type if necessary

if istimelock

data = ft_checkdata(data, 'datatype', 'timelock');

end

% do the general cleanup and bookkeeping at the end of the function

ft_postamble debug

ft_postamble previous data

ft_postamble provenance data

ft_postamble history data

ft_postamble savevar data